How to find the mean square error in matlab

Question

Is there a way to find the mean square error in matlab between 2 images A,B(say) in true color of dimension 256*256*3 ? The mathematical formula for a matrix say M1 and M2 is as under

mean sq err=1/n*n { summation (square[M1(i,j)-M2(i,j)])}

where i stands for row and j stands for column

Answer 1

Well, start writing! Eat a programming elephant (even the smallest ones) one byte at a time!

How do we form the difference of two images? First, convert them to doubles in case they are uint8 images, as is common. DO IT! TRY IT! Learn to write matlab code by doing so, and do it in pieces, so you can follow what you did.

First of all, you have not told us if this is to be a MSE over all three channels. Your formula says that we should end up with a different MSE for each of the red, green and blue channels.

double(M1) - double(M2)

Now, how would you form the square of each difference? Use the .^ operator.

(double(M1) - double(M2)).^2

Next, mean squared error implies that we take the mean over all rows and columns. A simple way to do this is with the mean function. This call takes the mean across the rows.

mean((double(M1) - double(M2)).^2,2)

And the next one takes the mean down the columns.

mean(mean((double(M1) - double(M2)).^2,2),1)

The result will be a 1x1x3 vector. Convert that into a 1x3 vector using the reshape function. (The squeeze function would help as well.) Packaging it all into one line, we get this...

MSE = reshape(mean(mean((double(M1) - double(M2)).^2,2),1),[1,3]);

If this seems complex to you, then you are best off splitting it into several lines, with comments that remind you what you did for later.

But the point is, you create an operation in matlab by breaking it down into manageable pieces.

EDIT:

In many cases, people want the RMSE (root-mean-squared-error) which has units the same as your original numbers. It is just the square root of the MSE.

Answer 2

Mean square error for each channel independently:

R1 = M1(:,:,1);
G1 = M1(:,:,2);
B1 = M1(:,:,3);

R2 = M2(:,:,1);
G2 = M2(:,:,2);
B2 = M2(:,:,3);

dR = int32(R1) - int32(R2);
dG = int32(G1) - int32(G2);
dB = int32(B1) - int32(B2);

mseR = mean(dR(:).^2);
mseG = mean(dG(:).^2);
mseB = mean(dB(:).^2);

If this is part of an image registration algorithm, you might want to do away with the squared term:

R1 = M1(:,:,1);
G1 = M1(:,:,2);
B1 = M1(:,:,3);

R2 = M2(:,:,1);
G2 = M2(:,:,2);
B2 = M2(:,:,3);

dR = int32(R1) - int32(R2);
dG = int32(G1) - int32(G2);
dB = int32(B1) - int32(B2);

errR = sum(abs(dR(:))); % 32bits sufficient for sum of 256x256 uint8 img.
errG = sum(abs(dG(:)));
errB = sum(abs(dB(:)));

sumErr = errR + errG + errB;

For additional performance, you might also want to consider converting to a single channel and spatially downsampling, although the effectiveness of this will depend on image content.