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个2068*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

更新:

在下面的评论中(稍微激烈的)辩论之后,我添加了由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(运算符/和\).