Matlab：如何在一组 2D 向量上向量化嵌套循环

Question

我有一个如下形式的函数：

function Out = DecideIfAPixelIsWithinAnEllipsoidalClass(pixel,means,VarianceCovarianceMatrix)  
   ellipsoid = (pixel-means)'*(VarianceCovarianceMatrix^(-1))*(pixel-means);  
   if ellipsoid <= 1
      Out = 1;
   else
      Out = 0;
   end
end  

我正在使用matlab进行遥感处理，我想对LandSatTM图像进行分类。这张图片有7个波段，大小为2048*2048。所以我将它们存储在3维2048*2048*7矩阵中。在这个函数中意味着是一个 7*1 矩阵，使用名为 ExtractStatisticalParameters 的函数中的类样本计算得出，而 VarianceCovarianceMatrix 实际上是一个 7*7 矩阵，您会看到：

ellipsoid = (pixel-means)'*(VarianceCovarianceMatrix^(-1))*(pixel-means);  

是椭圆体的方程。我的问题是每次你可以将一个像素（它是一个 7*1 向量，其中每一行是分隔带中像素的值）传递给这个函数，所以我需要像这样写一个循环：

for k1=1:2048  
   for k2=1:2048  
      pixel(:,1)=image(k1,k2,:); 
      Out = DecideIfAPixelIsWithinAnEllipsoidalClass(pixel,means,VarianceCovarianceMatrix);  
   end  
end  

您知道这将花费系统的大量时间和精力。您能建议我一种减少施加在系统上的压力的方法吗？

Answer 1

不需要循环！

pMinusMean = bsxfun( @minus, reshape( image, [], 7 ), means' ); %//' subtract means from all pixes
iCv = inv( arianceCovarianceMatrix );
ell = sum( (pMinusMean * iCv ) .* pminusMean, 2 ); % note the .* the second time!
Out = reshape( ell <= 1, size(image(:,:,1)) ); % out is 2048-by-2048 logical image

更新：

在下面的 cmets 中进行了（有些激烈的）辩论之后，我添加了 Rody Oldenhuis 所做的更正：

pMinusMean = bsxfun( @minus, reshape( image, [], 7 ), means' ); %//' subtract means from all pixes
ell = sum( (pMinusMean / varianceCovarianceMatrix ) .* pminusMean, 2 ); % note the .* the second time!
Out = reshape( ell <= 1, size(image(:,:,1)) );

此更改的关键问题是Matlab 的inv() 实现不佳，最好使用mldivide 和mrdivide（运算符/ 和\）代替。