add table pregeneration script for pmetric

experimental/skpdiff/SkPMetric.cpp

 #include <cmath>
 
 #include "SkBitmap.h"
 #include "skpdiff_util.h"
 #include "SkPMetric.h"
+#include "SkPMetricUtil_generated.h"
 
 struct RGB {
     float r, g, b;
 };
 
 struct LAB {
     float l, a, b;
 };
 
 template<class T>
@@ skipping 89 matching lines @@
     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);
+            f[i] = SkPMetricUtil::get_cube_root(f[i]);
         } 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);
+            rgb.r = SkPMetricUtil::get_gamma(row[x * 4 + 2]);
+            rgb.g = SkPMetricUtil::get_gamma(row[x * 4 + 1]);
+            rgb.b = SkPMetricUtil::get_gamma(row[x * 4 + 0]);
             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 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 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;
@@ skipping 23 matching lines @@
     const int matrixCount = sizeof(matrix) / sizeof(float);
     const int radius = matrixCount / 2;
 
     // Keep track of what rows are being operated on for quick access.
     float* rowPtrs[matrixCount]; // Because matrixCount is constant, this won't create a VLA
     for (int y = radius; y < matrixCount; y++) {
         rowPtrs[y] = imageL->getRow(y - radius);
     }
     float* writeRow = outImageL->getRow(0);
 
-
     for (int y = 0; y < imageL->height; y++) {
         for (int x = 0; x < imageL->width; x++) {
             float lSum = 0.0f;
             for (int xx = -radius; xx <= radius; xx++) {
                 int nx = x;
                 int ny = y;
 
                 // We mirror at edges so that edge pixels that the filter weighting still makes
                 // sense.
                 if (vertical) {
@@ skipping 47 matching lines @@
 
     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;
     }
 
+    // Contrast sensitivity is based on image dimensions. Therefore it cannot be statically
+    // generated.
+    float * contrastSensitivityTable = SkNEW_ARRAY(float, maxLevels * 1000);

[djsollen, 2013/07/16 13:55:11]

float*

+    for (int levelIndex = 0; levelIndex < maxLevels; levelIndex++) {
+        for (int csLum = 0; csLum < 1000; csLum++) {
+           contrastSensitivityTable[levelIndex * 1000 + csLum] =
+           contrast_sensitivity(cyclesPerDegree[levelIndex], (float)csLum / 10.0f + 1e-5f);
+       }
+    }
+
     // Compute G - The convolved lum for the baseline
     for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
         convolve(baselineL.getLayer(levelIndex - 1), false, &scratchImageL);
         convolve(&scratchImageL, true, baselineL.getLayer(levelIndex));
     }
     for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
         convolve(testL.getLayer(levelIndex - 1), false, &scratchImageL);
         convolve(&scratchImageL, true, testL.getLayer(levelIndex));
     }
 
@@ skipping 41 matching lines @@
 
                 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));
+                float mask = SkPMetricUtil::get_visual_mask(contrast[levelIndex] *
+                                              contrastSensitivityTable[levelIndex * 1000 + (int)(lAdapt * 10.0)]);

[djsollen, 2013/07/16 13:55:11]

100 char limit

 
                 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)) {
+            if (fabsf(lBaseline - lTest) > F * SkPMetricUtil::get_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) {
@@ skipping 10 matching lines @@
             if (isFailure) {
                 failures++;
                 poi->push()->set(x, y);
             }
         }
     }
 
     SkDELETE_ARRAY(cyclesPerDegree);
     SkDELETE_ARRAY(contrast);
     SkDELETE_ARRAY(thresholdFactorFrequency);
+    SkDELETE_ARRAY(contrastSensitivityTable);
     return 1.0 - (double)failures / (width * height);
 }
 
 const char* SkPMetric::getName() {
     return "perceptual";
 }
 
 int SkPMetric::queueDiff(SkBitmap* baseline, SkBitmap* test) {
+    double startTime = get_seconds();
     int diffID = fQueuedDiffs.count();
-    double startTime = get_seconds();
     QueuedDiff& diff = fQueuedDiffs.push_back();
     diff.result = 0.0;
 
     // Ensure the images are comparable
     if (baseline->width() != test->width() || baseline->height() != test->height() ||
                     baseline->width() <= 0 || baseline->height() <= 0) {
         diff.finished = true;
         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, &diff.poi);
 
     SkDebugf("Time: %f\n", (get_seconds() - startTime));
 
     return diffID;
 }
 
 
 void SkPMetric::deleteDiff(int id) {
-   fQueuedDiffs[id].poi.reset();
+   //fQueuedDiffs[id].poi.reset();
 }
 
 bool SkPMetric::isFinished(int id) {
     return fQueuedDiffs[id].finished;
 }
 
 double SkPMetric::getResult(int id) {
     return fQueuedDiffs[id].result;
 }
 
 int SkPMetric::getPointsOfInterestCount(int id) {
     return fQueuedDiffs[id].poi.count();
 }
 
 SkIPoint* SkPMetric::getPointsOfInterest(int id) {
     return fQueuedDiffs[id].poi.begin();
 }