SSIM index in Mathematica

Question

I was wondering, how one could compute the SSIM index in Mma ? The The structural SIMilarity (SSIM) index, which is a known method for measuring the similarity between two images.

Here is the original paper: http://www.cns.nyu.edu/pub/eero/wang03-reprint.pdf

And this is a nice summary of the essential formula: http://ch.mathworks.com/help/images/ref/ssim.html

The SSIM is an index measuring the structural similarity between two images. It is valued between -1 and 1. When two images are nearly identical, their SSIM is close to 1.

Formula computing the SSIM between two sequences seq1 and seq2 at a given pixel P:

$SSIM(P)= \frac{2*\mu_1(P)*\mu_2(P) + C1}{\mu_1(P)^2 + \mu_2(P)^2 + C1} * \frac{2*cov(P) + C2}{s1(P)^2 + s2(P)^2 + C2}$

With:

mu1(P) and mu2(P): mean value of seq1 and seq2 computed over a small XY window located around P

s1(P) and s2(P): standard deviation of seq1 and seq2 computed over the same window

cov(P): covariance between seq1 and seq2 computed over the same window

C1 = (K1*L)^2: regularization constant (should be as small as possible)

C2 = (K2*L)^2: regularization constant (should be as small as possible)

K1, K2: regularization parameters (must be >0)

L: dynamic range of the pixel values (example: L=255 if the sequence is 8 bit encoded)

Answer 1

Attempt to reconstruct the original author's stuff

I was curious enough to have a look into the Matlab source code ssim.m provided by the original authors Wang et al. 2004. I have tried my best to exactly re-implement this code with Wolfram Language using MMA 9 (have checked it under MMA 11.3 as well and have re-arrangend the listing below for better readability):

MyMSSIM[
   img1raw_Image,
   img2raw_Image,
   window_: GaussianMatrix[{Table[(11 - 1)/2, {2}], 1.5}, 
     Method -> "Gaussian"],
   k1_: 0.01, k2_: 0.03] :=
  Module[{
    f = Max[1, Round[Min[Sequence @@ ImageDimensions[img1raw]]/256]],
    img1in = ColorConvert[Image[img1raw, "Real"], "Grayscale"],
    img2in = ColorConvert[Image[img2raw, "Real"], "Grayscale"],
    mssimmap, mssim},

   If[f > 1,(*if subsampling needs to be done*)
    {img1in, img2in} =
     ImageCorrelate[
        ImagePad[#,(*to mimic the Matlab imfilter edge treatment*)
         NestList[Reverse, {Floor[(f - 1)/2], Ceiling[(f - 1)/2]}, 1],
         "Reversed"],
        ConstantArray[1./Times[f, f], {f, f}],
        Padding -> None] & /@ {img1in, img2in};
    {img1in, img2in} = 
     ImageTake[#, {1, -1, f}, {1, -1, f}] & /@ {img1in, img2in};
    ];

   mssimmap =
    Function[{img1, img2},
      Function[{mu1, mu2},
        Function[{mu1mu1, mu2mu2, mu1mu2},
          ImageApply[#1/#2 &,
           {ImageMultiply[
             ImageAdd[
              ImageMultiply[mu1mu2, 2.],
              k1^2],
             ImageAdd[
              ImageMultiply[
               ImageSubtract[
                ImageCorrelate[
                 ImageMultiply[img1, img2],
                 window, Padding -> None],
                mu1mu2],
               2.],
              k2^2]
             ],
            ImageMultiply[
             ImageAdd[
              ImageAdd[mu1mu1, mu2mu2],
              k1^2],
             ImageAdd[
              ImageAdd[
               ImageSubtract[
                ImageCorrelate[
                 ImageMultiply[img1, img1],
                 window, Padding -> None],
                mu1mu1],
               ImageSubtract[
                ImageCorrelate[
                 ImageMultiply[img2, img2],
                 window, Padding -> None],
                mu2mu2]],
              k2^2]
             ]
            }
           ]
          ][
         ImageMultiply[mu1, mu1],
         ImageMultiply[mu2, mu2],
         ImageMultiply[mu1, mu2]
         ]
        ][
       ImageCorrelate[img1, window, Padding -> None],
       ImageCorrelate[img2, window, Padding -> None]
       ]
      ][img1in, img2in
     ];
   mssim = ImageMeasurements[mssimmap, "Mean"];
   Return@mssim;
   ];

I more or less consequently was using the built-in routines to handle Image data.

I have checked the examples given under http://www.cns.nyu.edu/~lcv/ssim:

urlbase = "http://www.cns.nyu.edu/~lcv/ssim/";
imagelist = {"einstein.gif", "meanshift.gif", "contrast.gif", "impulse.gif", "blur.gif", "jpg.gif"};
Function[{name1, name2}, 
    Column[{Last@StringSplit[name2, "/"], Show@#2, 
        "Mean SSIM: " <> ToString@MyMSSIM[#1, #2]}, Center] &[
     Import@name1, Import@name2]][urlbase <> First@imagelist, 
   urlbase <> #] & /@ imagelist

Please re-check and compare by your own. I also compared a plenty of mean SSIM values for image pairs one-by-one with the output of the Matlab code and found no differences within the precision of the displayed mantissa.

What I skipped in the original Matlab code was any check of minimum image sizes and something like that, and I also did not cover the unlikely case that someone will operate with either k1 or k2 set to 0. Their presets (0.01 and 0.03, resp.) as well as parameters for the Gaussian smoothing kernel window were also taken from that code. As you can see in my code, I also mimic the strange smoothing using a box filter in case of a sub-sampling of images larger than 383 pixels at their short edge. And I generally did exactly the same edge handling using certain Padding settings.

Please look at two basic examples (part of Wang's demo as the other examples given above):

MyMSSIM[Import@"http://www.cns.nyu.edu/~lcv/ssim/einstein.gif", 
 Import@"http://www.cns.nyu.edu/~lcv/ssim/meanshift.gif"]

gives

0.988359

for and

while

MyMSSIM[Import@"http://www.cns.nyu.edu/~lcv/ssim/einstein.gif", 
 Import@"http://www.cns.nyu.edu/~lcv/ssim/jpg.gif"]

gives

0.662363

for and

There is a very small change I made, as the original Matlab code is intended for scalar images and cannot handle RGB images, in case someone inserts RGB images here, I convert them to real valued luminance images. As values now are in the interval [0.,1.], I could cancel the parameter L in the original code, indicating the maximum pixel value. BTW, output range is between -1 and 1, while for typical image pairs the result will be positive.

What I did not consider so far was any multi-scale implementation of the Structural SIMilarity (SSIM), called MS-SSIM msssim.m

Answer 2

I just translate the formula:

$$\begin{aligned} L(X,Y) &= \frac{2u_Xu_Y + C_1}{{u_X}^2+{u_Y}^2+C_1}\\ C(X,Y) &= \frac{2\sigma_X\sigma_Y+C_2}{{\sigma_X}^2+{\sigma_Y}^2 + C_2}\\ S(X,Y) &= \frac{\sigma_{XY} + C_3}{\sigma_X\sigma_Y + C_3} \end{aligned}$$

$$\mathrm{SSIM}= L(X,Y)^\alpha ×C(X,Y)^\beta ×S(X,Y)^\gamma$$

Where:

$$\begin{cases} C_1=(K_1\cdot L)^2, C_2=(K_2\cdot L)^2, C_3=C_2/2\\ K_1=0.01, K_2=0.03, L=255\\ \alpha =\beta =\gamma=1 \end{cases}$$

So we get:

$$\hbox{SSIM}(X,Y)={\frac{(2\mu_x\mu _{y}+C_{1})(2\sigma _{xy}+C_2)}{({\mu_x}^2+{\mu _y}^2+C_1)({\sigma_x}^2+{\sigma_y}^2+C_2)}}$$

MSSIM[img1_Image, img2_Image, n_Integer] := SSIM @@ Transpose[ImagePartition[#, Scaled[1 / n]]& /@ {img1, img2}] / n^2;
SSIM[img1_Image, img2_Image] := Module[
    {i1, i2, mx, my, vx, vy, c1, c2, cov},
    i1 = N@Flatten@ImageData[img1, "Byte"];
    i2 = N@Flatten@ImageData[img2, "Byte"];
    {mx, my} = Mean /@ {i1, i2};
    {vx, vy} = StandardDeviation /@ {i1, i2};
    {c1, c2} = {255 * 0.01, 255 * 0.03}^2;
    cov = Covariance[i1, i2];
    (2mx * my + c1)(2cov + c2) / (mx^2 + my^2 + c1) / (vx^2 + vy^2 + c2)
]

But there's something wrong, it not match the original answer.

http://www.cns.nyu.edu/~lcv/ssim/#test

img1=Import["http://www.cns.nyu.edu/~lcv/ssim/einstein.gif"];
img2=Import["http://www.cns.nyu.edu/~lcv/ssim/meanshift.gif"];
SSIM[img1,img2]

It should be 0.988 but I got 0.997 and I can't figure out what's wrong.

---

Update:

@UDB gave the right answer, so I have to correct my code.

Options[SSIM] = {
    "C1" -> 0.01^2, "C2" -> 0.03^2,
    "Window" -> GaussianMatrix[{{(11 - 1) / 2, (11 - 1) / 2}, 1.5}, Method -> "Gaussian"]
};
SSIM[img1_Image, img2_Image, OptionsPattern[]] := Module[
    {c1, c2, window, mx, my, vx, vy, cov, r},
    {c1, c2, window} = OptionValue[{"C1", "C2", "Window"}];
    mx = ImageCorrelate[img1, window, Padding -> None];
    my = ImageCorrelate[img2, window, Padding -> None];
    vx = ImageCorrelate[img1^2, window, Padding -> None] - mx^2;
    vy = ImageCorrelate[img2^2, window, Padding -> None] - my^2;
    cov = ImageCorrelate[img1 * img2, window, Padding -> None] - mx * my;
    r = (2mx * my + c1) / (mx^2 + my^2 + c1) * (2cov + c2) / (vx + vy + c2);
    Mean@ImageMeasurements[r, "Mean"]
]

The sixth place after the decimal point is a little different.

I don't know if ImageCorrelate is an efficient function, and I am glad to see some more efficient code.