Eliminating for-loops (compact notation?)

I am trying to compare fMRI images, but I seem to have quite an inefficient comparison algorithm.I have heard that efficiency gains can be made by using the slightly more native compact notation, but I have no idea where to start with this. Here is the code, it will end up being run several million times.

for xInd=1:xMax
    for yInd=1:yMax
        for zInd=1:zMax
            if ~isnan(imageArray_1(1,xInd,yInd,zInd))
                statImage(xInd,yInd,zInd)=abs(log(mean(imageArray_1(:,xInd,yInd,zInd)./imageNormalisationDivisor_1))-log(mean(imageArray_2(:,xInd,yInd,zInd)./imageNormalisationDivisor_2)));
            end    
        end
    end
end

Essentially there are 10 three-dimensional arrays per condition, and the difference in log mean value at each voxel across each condition is compared. Each image is normalised by dividing it by the mean value.

It's not a dealbreaker if it can't get much more efficient, I can leave it running for the day or so it would need, but I have a feeling that there are improvements that could be made.

Thanks, Matt

Answer 1

[Thorsten]

Oh, that's easy (just kidding):

 X = squeeze(abs(log(mean(imageArray_1, 1)./mean(imageArray_2,1)*imageNormalisationDivisor_2/imageNormalisationDivisor_1)));

Note that the length of this one-liner is partly due to your long variable names. Compare

 X = squeeze(abs(log(mean(I1, 1)./mean(I2,1)*s2/s1)));

To compare both solutions, use

 max(abs(statImage(:) -X(:)))

That's not zero, but should be a small value in the range of eps, the machine precision.

[Matt Moore]

Thanks Thorsten! That's a good start, although I didn't make it clear in my original post that each image is normalised using its own mean, not the overall mean (i.e. imageNormalisationDivisor is 2D). I think I should be able to put something together from that skeleton though.

Edit: My solution has been to generate a 3D matrix for each scan containing the single normalisation divisor repeated however many times, and then using it thus:

      norm1=nanmean(imageArray_1./imageNormalisationDivisor_1,1);
      norm2=nanmean(imageArray_2./imageNormalisationDivisor_2,1);

      statImage = abs(log(squeeze(norm1./norm2)));

Although it takes a while to generate this normalisation matrix (I would be open to suggestions on how to cut down that time:) ), that does not need to be repeated nearly as often, and the code is overall faster by a factor of 10.