https://github.com/up42/image-similarity-measures

https://ece.uwaterloo.ca/~z70wang/publications/quality_2c.pdf

https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

<code C++>
const int width = 200;
const int height = 200;

inline int sum (int v)
{
   return (v >> 16) & 255 + (v >> 8) & 255 + v & 255;
}

double meanx,  meany, varx, vary, varxy;

void f (int * imgx, int * imgy)
{
    int n = 0;
    long sumx = 0, sumxx = 0, sumy = 0, sumyy = 0, sumxy = 0;

    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            int valx = imgx [i*width+j];
            int x = (valx >> 16) & 255 + (valx >> 8) & 255 + valx & 255;
            int valy = imgy [i*width+j];
            int y = sum (valy);
            n ++;
            sumx += x;
            sumxx += x*x;
            sumy += y;
            sumyy += y*y;
            sumxy += x*y;
        }

        meanx = sumx/n;
        meany = sumy/n;
        varx  = (sumxx - sumx * sumx / n ) / (n-1);
        vary  = (sumyy - sumy * sumy / n ) / (n-1);
        varxy = (sumxy - sumx * sumy / n ) / n;
}

int main(int argc, char *argv[])
{
    int * imgx = new int [width*height];
    int * imgy = new int [width*height];
    f (imgx, imgy);
}
</code>



<code C++>
class Image;

void f (Image * imgx, Image * imgy)
{
    int n = 0;
    double sumx = 0, sumxx = 0, sumy = 0, sumyy = 0, sumxy = 0;
    
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            x = imgx [i,j];
            y = imgy [i,j];
            n ++;
            sumx += x;
            sumxx += x*x;
            sumy += y;
            sumyy += y*y;
            sumxy += x*y;
        }
                  
                    double meanx = sumx/n;                    
                    double meany = sumy/n;                    
                    double varx  = (sumxx - sumx * sumx / n ) / (n-1);                    
                    double vary  = (sumyy - sumy * sumy / n ) / (n-1);                    
                    double varxy = (sumxy - sumx * sumy / n ) / n / (n-1);
}
</code>