add yee's perceptual metric

experimental/skpdiff/SkPMetric.cpp

Index: experimental/skpdiff/SkPMetric.cpp
diff --git a/experimental/skpdiff/SkPMetric.cpp b/experimental/skpdiff/SkPMetric.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7fbdec87aa4a5e8957f828ed41cba2399176776b
--- /dev/null
+++ b/experimental/skpdiff/SkPMetric.cpp
@@ -0,0 +1,417 @@
+#include <cmath>
+
+#include "SkBitmap.h"
+#include "skpdiff_util.h"
+#include "SkPMetric.h"
+
+struct RGB {
+    float r, g, b;
+};
+
+struct LAB {
+    float l, a, b;
+};
+
+template<class T>
+struct Image2D {
+    int width;
+    int height;
+    T* image;
+
+    Image2D(int w, int h)
+        : width(w),
+          height(h) {
+        SkASSERT(w > 0);
+        SkASSERT(h > 0);
+        image = SkNEW_ARRAY(T, w * h);
+    }
+
+    ~Image2D() {
+        SkDELETE_ARRAY(image);
+    }
+
+    void readPixel(int x, int y, T* pixel) const {
+        SkASSERT(x >= 0);
+        SkASSERT(y >= 0);
+        SkASSERT(x < width);
+        SkASSERT(y < height);
+        *pixel = image[y * width + x];
+    }
+
+    void writePixel(int x, int y, const T& pixel) {
+        SkASSERT(x >= 0);
+        SkASSERT(y >= 0);
+        SkASSERT(x < width);
+        SkASSERT(y < height);
+        image[y * width + x] = pixel;
+    }
+};
+
+typedef Image2D<float> ImageL;
+typedef Image2D<RGB> ImageRGB;
+typedef Image2D<LAB> ImageLAB;
+
+template<class T>

[bsalomon, 2013/06/27 18:31:15]

It seems like the template param should be restricted to be an Image2D.

Maybe make T be one of float, RGB, or LAB and instead of:

T** image;

have

Image2D<T>** image;

But also does this class add much value over an array or Image2Ds? The name
makes it sound like a volume texture, but I don't think that's really what it
is.

[Zach Reizner, 2013/06/27 18:48:49]

I like the idea of moving the template params.

The reason I have this class over an array is for bounds checking. I was having
a hard time debugging my segmentation faults, so I created this and the Image2D.
As an added benefit, I didn't have to pass so array lengths to the functions
below, and my code wasn't littered with complex pointer arithmetic.

+struct Image3D
+{
+    int slices;
+    T** image;
+
+    Image3D(int w, int h, int s)
+        : slices(s) {
+        SkASSERT(s > 0);
+        image = SkNEW_ARRAY(T*, s);
+        for (int sliceIndex = 0; sliceIndex < slices; sliceIndex++) {
+            image[sliceIndex] = SkNEW_ARGS(T, (w, h));
+        }
+    }
+
+    ~Image3D() {
+        for (int sliceIndex = 0; sliceIndex < slices; sliceIndex++) {
+            SkDELETE(image[sliceIndex]);
+        }
+        SkDELETE_ARRAY(image);
+    }
+
+    T* getLayer(int z) const {
+        SkASSERT(z >= 0);
+        SkASSERT(z < slices);
+        return image[z];
+    }
+};
+
+typedef Image3D<ImageL> ImageL3D;
+
+
+#define MAT_ROW_MULT(rc,gc,bc) r*rc + g*gc + b*bc
+
+
+void adobergb_to_cielab(float r, float g, float b, LAB* lab) {

[bsalomon, 2013/06/27 18:31:15]

Does this mean the input is in Adobe RGB?
http://en.wikipedia.org/wiki/Adobe_RGB_color_space

[Zach Reizner, 2013/06/27 18:48:49]

The input is assumed to be Adobe RGB. Adobe RGB has a defined white point that
matches many typical monitors.

+    // Conversion of Adobe RGB to XYZ taken from from "Adobe RGB (1998) ColorImage Encoding"
+    // URL:http://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf
+    // Section: 4.3.5.3
+    // See Also: http://en.wikipedia.org/wiki/Adobe_rgb
+    float x = MAT_ROW_MULT(0.57667f, 0.18556f, 0.18823f);
+    float y = MAT_ROW_MULT(0.29734f, 0.62736f, 0.07529f);
+    float z = MAT_ROW_MULT(0.02703f, 0.07069f, 0.99134f);
+
+    // The following is the white point in XYZ, so it's simply the row wise addition of the above
+    // matrix.
+    const float xw = 0.5767f + 0.185556f + 0.188212f;
+    const float yw = 0.297361f + 0.627355f + 0.0752847f;
+    const float zw = 0.0270328f + 0.0706879f + 0.991248f;
+
+    // This is the XYZ color point relative to the white point
+    float f[3] = { x / xw, y / yw, z / zw };
+
+    // Conversion from XYZ to LAB taken from
+    // http://en.wikipedia.org/wiki/CIELAB#Forward_transformation
+    for (int i = 0; i < 3; i++) {
+        if (f[i] >= 0.008856f) {
+            f[i] = powf(f[i], 1.0f / 3.0f);
+        } else {
+            f[i] = 7.787f * f[i] + 4.0f / 29.0f;
+        }
+    }
+    lab->l = 116.0f * f[1] - 16.0f;
+    lab->a = 500.0f * (f[0] - f[1]);
+    lab->b = 200.0f * (f[1] - f[2]);
+}
+
+/// Converts a 8888 bitmap to LAB color space and puts it into the output
+static void bitmap_to_cielab(const SkBitmap* bitmap, ImageLAB* outImageLAB) {
+    SkASSERT(bitmap->config() == SkBitmap::kARGB_8888_Config);
+
+    int width = bitmap->width();
+    int height = bitmap->height();
+    SkASSERT(outImageLAB->width == width);
+    SkASSERT(outImageLAB->height == height);
+
+    bitmap->lockPixels();
+    RGB rgb;
+    LAB lab;
+    for (int y = 0; y < height; y++) {
+        unsigned char* row = (unsigned char*)bitmap->getAddr(0, y);
+        for (int x = 0; x < width; x++) {
+            // Perform gamma correction which is assumed to be 2.2
+            rgb.r = powf(row[x * 4 + 2] / 255.0f, 2.2f);
+            rgb.g = powf(row[x * 4 + 1] / 255.0f, 2.2f);
+            rgb.b = powf(row[x * 4 + 0] / 255.0f, 2.2f);
+            adobergb_to_cielab(rgb.r, rgb.g, rgb.b, &lab);
+            outImageLAB->writePixel(x, y, lab);
+        }
+    }
+    bitmap->unlockPixels();
+}
+
+// From Barten SPIE 1989
+static float contrast_sensitivity(float cyclesPerDegree, float luminance) {
+    float a = 440.0f * powf(1.0f + 0.7f / luminance, -0.2f);
+    float b = 0.3f * powf(1 + 100.0 / luminance, 0.15f);
+    return a *
+           cyclesPerDegree *
+           expf(-b * cyclesPerDegree) *
+           sqrtf(1.0f + 0.06f * expf(b * cyclesPerDegree));
+}
+
+// From Daly 1993
+static float visual_mask(float contrast) {
+    float x = powf(392.498f * contrast, 0.7f);
+    x = powf(0.0153f * x, 4.0f);
+    return powf(1.0f + x, 0.25f);
+}
+
+// From Ward Larson Siggraph 1997
+static float threshold_vs_intensity(float adaptationLuminance) {
+    float logLum = log10f(adaptationLuminance);
+    float x;
+    if (logLum < -3.94f) {
+        x = -2.86f;
+    } else if (logLum < -1.44f) {
+        x = powf(0.405f * logLum + 1.6f, 2.18) - 2.86f;
+    } else if (logLum < -0.0184f) {
+        x = logLum - 0.395f;
+    } else if (logLum < 1.9f) {
+        x = powf(0.249f * logLum + 0.65f, 2.7f) - 0.72f;
+    } else {
+        x = logLum - 1.255f;
+    }
+    return powf(10.0f, x);
+}
+
+/// Simply takes the L channel from the input and puts it into the output
+static void lab_to_l(const ImageLAB* imageLAB, ImageL* outImageL) {
+    for (int y = 0; y < imageLAB->height; y++) {
+        for (int x = 0; x < imageLAB->width; x++) {
+            LAB lab;
+            imageLAB->readPixel(x, y, &lab);
+            outImageL->writePixel(x, y, lab.l);
+        }
+    }
+}
+
+/// Convolves an image with the given filter in one direction and saves it to the output image
+static void convolve(const ImageL* imageL,
+                     bool vertical, const float* matrix, int radius,
+                     ImageL* outImageL) {
+    SkASSERT(imageL->width == outImageL->width);
+    SkASSERT(imageL->height == outImageL->height);
+    for (int y = 0; y < imageL->height; y++) {
+        for (int x = 0; x < imageL->width; x++) {
+            float lSum = 0.0f;
+            float l;
+            for (int xx = -radius; xx <= radius; xx++) {
+                int nx = x;
+                int ny = y;
+                if (vertical) {
+                    ny += xx;
+                } else {
+                    nx += xx;
+                }
+                if (nx < 0) {

[bsalomon, 2013/06/27 18:31:15]

should these four ifs be done inside the vertical/horizontal if/else?

You mirror at the borders?

+                    nx = -nx;
+                }
+                if (ny < 0) {
+                    ny = -ny;
+                }
+                if (nx >= imageL->width) {
+                    nx = imageL->width + (imageL->width - nx - 1);
+                }
+                if (ny >= imageL->height) {
+                    ny = imageL->height + (imageL->height - ny - 1);
+                }
+                imageL->readPixel(nx, ny, &l);
+                float weight = matrix[xx + radius];
+                lSum += l * weight;
+            }
+            outImageL->writePixel(x, y, lSum);
+        }
+    }
+}
+
+float pmetric(const ImageLAB* baselineLAB, const ImageLAB* testLAB) {
+    int width = baselineLAB->width;
+    int height = baselineLAB->height;
+    int maxLevels = (int)log2(width < height ? width : height);
+
+    const float fov = M_PI / 180.0f * 45.0f;
+    float contrastSensitivityMax = contrast_sensitivity(3.248f, 100.0f);
+    float pixelsPerDegree = width / (2.0f * tanf(fov * 0.5f) * 180.0f / M_PI);
+
+    ImageL3D baselineL(width, height, maxLevels);
+    ImageL3D testL(width, height, maxLevels);
+    ImageL scratchImageL(width, height);
+    float* cyclesPerDegree = SkNEW_ARRAY(float, maxLevels);
+    float* thresholdFactorFrequency = SkNEW_ARRAY(float, maxLevels - 2);
+    float* contrast = SkNEW_ARRAY(float, maxLevels - 2);
+
+    lab_to_l(baselineLAB, baselineL.getLayer(0));
+    lab_to_l(testLAB, testL.getLayer(0));
+
+    // Compute cpd - Cycles per degree on the pyramid
+    cyclesPerDegree[0] = 0.5f * pixelsPerDegree;
+    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
+        cyclesPerDegree[levelIndex] = cyclesPerDegree[levelIndex - 1] * 0.5f;
+    }
+
+    const float filterMatrix[] = { 0.05f, 0.25f, 0.4f, 0.25f, 0.05f };
+    // Compute G - The convolved lum for the baseline
+    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
+        convolve(baselineL.getLayer(levelIndex - 1), false, filterMatrix, 2, &scratchImageL);
+        convolve(&scratchImageL, true, filterMatrix, 2, baselineL.getLayer(levelIndex));
+    }
+    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
+        convolve(testL.getLayer(levelIndex - 1), false, filterMatrix, 2, &scratchImageL);
+        convolve(&scratchImageL, true, filterMatrix, 2, testL.getLayer(levelIndex));
+    }
+
+    // Compute F_freq - The elevation f
+    for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
+        float cpd = cyclesPerDegree[levelIndex];
+        thresholdFactorFrequency[levelIndex] = contrastSensitivityMax /
+                                               contrast_sensitivity(cpd, 100.0f);
+    }
+
+    int failures = 0;
+    // Calculate F
+    for (int y = 0; y < height; y++) {
+        for (int x = 0; x < width; x++) {
+            float lBaseline;
+            float lTest;
+            baselineL.getLayer(0)->readPixel(x, y, &lBaseline);
+            testL.getLayer(0)->readPixel(x, y, &lTest);
+
+            float avgLBaseline;
+            float avgLTest;
+            baselineL.getLayer(maxLevels - 1)->readPixel(x, y, &avgLBaseline);
+            testL.getLayer(maxLevels - 1)->readPixel(x, y, &avgLTest);
+
+            float lAdapt = 0.5f * (avgLBaseline + avgLTest);
+            if (lAdapt < 1e-5) {
+                lAdapt = 1e-5;
+            }
+
+            float contrastSum = 0.0f;
+            for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
+                float baselineL0, baselineL1, baselineL2;
+                float testL0, testL1, testL2;
+                baselineL.getLayer(levelIndex + 0)->readPixel(x, y, &baselineL0);
+                testL.    getLayer(levelIndex + 0)->readPixel(x, y, &testL0);
+                baselineL.getLayer(levelIndex + 1)->readPixel(x, y, &baselineL1);
+                testL.    getLayer(levelIndex + 1)->readPixel(x, y, &testL1);
+                baselineL.getLayer(levelIndex + 2)->readPixel(x, y, &baselineL2);
+                testL.    getLayer(levelIndex + 2)->readPixel(x, y, &testL2);
+
+                float baselineContrast1 = fabsf(baselineL0 - baselineL1);
+                float testContrast1     = fabsf(testL0 - testL1);
+                float numerator = (baselineContrast1 > testContrast1) ?
+                                   baselineContrast1 : testContrast1;
+
+                float baselineContrast2 = fabsf(baselineL2);
+                float testContrast2     = fabsf(testL2);
+                float denominator = (baselineContrast2 > testContrast2) ?
+                                    baselineContrast2 : testContrast2;
+
+                // Avoid divides by close to zero
+                if (denominator < 1e-5) {
+                    denominator = 1e-5;
+                }
+
+                contrast[levelIndex] = numerator / denominator;
+                contrastSum += contrast[levelIndex];
+            }
+
+            if (contrastSum < 1e-5) {
+                contrastSum = 1e-5;
+            }
+
+            float F = 0.0f;
+            for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
+                float mask = visual_mask(contrast[levelIndex] *
+                             contrast_sensitivity(cyclesPerDegree[levelIndex], lAdapt));
+
+                F += contrast[levelIndex] +
+                     thresholdFactorFrequency[levelIndex] * mask / contrastSum;
+            }
+
+            if (F < 1.0f) {
+                F = 1.0f;
+            }
+
+            if (F > 10.0f) {
+                F = 10.0f;
+            }
+
+
+            bool isFailure = false;
+            if (fabsf(lBaseline - lTest) > F * threshold_vs_intensity(lAdapt)) {
+                isFailure = true;
+            } else {
+                LAB baselineColor;
+                LAB testColor;
+                baselineLAB->readPixel(x, y, &baselineColor);
+                testLAB->readPixel(x, y, &testColor);
+                float contrastA = baselineColor.a - testColor.a;
+                float contrastB = baselineColor.b - testColor.b;
+                float colorScale = 1.0f;
+                if (lAdapt < 10.0f) {
+                    colorScale = lAdapt / 10.0f;
+                }
+                colorScale *= colorScale;
+
+                if ((contrastA * contrastA + contrastB * contrastB) * colorScale > F)
+                {
+                    isFailure = true;
+                }
+            }
+
+            if (isFailure) {
+                failures++;
+            }
+        }
+    }
+
+    SkDELETE_ARRAY(cyclesPerDegree);
+    SkDELETE_ARRAY(contrast);
+    SkDELETE_ARRAY(thresholdFactorFrequency);
+    return (double)failures;
+}
+
+const char* SkPMetric::getName() {
+    return "perceptual";
+}
+
+int SkPMetric::queueDiff(SkBitmap* baseline, SkBitmap* test) {
+    int diffID = fQueuedDiffs.count();
+    double startTime = get_seconds();
+    QueuedDiff* diff = fQueuedDiffs.push();
+
+    // Ensure the images are comparable
+    if (baseline->width() != test->width() || baseline->height() != test->height() ||
+                    baseline->width() <= 0 || baseline->height() <= 0) {
+        diff->finished = true;
+        diff->result = 0.0;
+        return diffID;
+    }
+
+    ImageLAB baselineLAB(baseline->width(), baseline->height());
+    ImageLAB testLAB(baseline->width(), baseline->height());
+
+    bitmap_to_cielab(baseline, &baselineLAB);
+    bitmap_to_cielab(test, &testLAB);
+
+    diff->result = pmetric(&baselineLAB, &testLAB);
+
+    SkDebugf("Time: %f\n", (get_seconds() - startTime));
+
+    return diffID;
+}
+
+
+bool SkPMetric::isFinished(int id) {
+    return fQueuedDiffs[id].finished;
+}
+
+double SkPMetric::getResult(int id) {
+    return fQueuedDiffs[id].result;
+}