1#include <cmath>
2
3#include "SkBitmap.h"
4#include "skpdiff_util.h"
5#include "SkPMetric.h"
6#include "SkPMetricUtil_generated.h"
7
8struct RGB {
9    float r, g, b;
10};
11
12struct LAB {
13    float l, a, b;
14};
15
16template<class T>
17struct Image2D {
18    int width;
19    int height;
20    T* image;
21
22    Image2D(int w, int h)
23        : width(w),
24          height(h) {
25        SkASSERT(w > 0);
26        SkASSERT(h > 0);
27        image = SkNEW_ARRAY(T, w * h);
28    }
29
30    ~Image2D() {
31        SkDELETE_ARRAY(image);
32    }
33
34    void readPixel(int x, int y, T* pixel) const {
35        SkASSERT(x >= 0);
36        SkASSERT(y >= 0);
37        SkASSERT(x < width);
38        SkASSERT(y < height);
39        *pixel = image[y * width + x];
40    }
41
42    T* getRow(int y) const {
43        return &image[y * width];
44    }
45
46    void writePixel(int x, int y, const T& pixel) {
47        SkASSERT(x >= 0);
48        SkASSERT(y >= 0);
49        SkASSERT(x < width);
50        SkASSERT(y < height);
51        image[y * width + x] = pixel;
52    }
53};
54
55typedef Image2D<float> ImageL;
56typedef Image2D<RGB> ImageRGB;
57typedef Image2D<LAB> ImageLAB;
58
59template<class T>
60struct ImageArray
61{
62    int slices;
63    Image2D<T>** image;
64
65    ImageArray(int w, int h, int s)
66        : slices(s) {
67        SkASSERT(s > 0);
68        image = SkNEW_ARRAY(Image2D<T>*, s);
69        for (int sliceIndex = 0; sliceIndex < slices; sliceIndex++) {
70            image[sliceIndex] = SkNEW_ARGS(Image2D<T>, (w, h));
71        }
72    }
73
74    ~ImageArray() {
75        for (int sliceIndex = 0; sliceIndex < slices; sliceIndex++) {
76            SkDELETE(image[sliceIndex]);
77        }
78        SkDELETE_ARRAY(image);
79    }
80
81    Image2D<T>* getLayer(int z) const {
82        SkASSERT(z >= 0);
83        SkASSERT(z < slices);
84        return image[z];
85    }
86};
87
88typedef ImageArray<float> ImageL3D;
89
90
91#define MAT_ROW_MULT(rc,gc,bc) r*rc + g*gc + b*bc
92
93static void adobergb_to_cielab(float r, float g, float b, LAB* lab) {
94    // Conversion of Adobe RGB to XYZ taken from from "Adobe RGB (1998) ColorImage Encoding"
95    // URL:http://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf
96    // Section: 4.3.5.3
97    // See Also: http://en.wikipedia.org/wiki/Adobe_rgb
98    float x = MAT_ROW_MULT(0.57667f, 0.18556f, 0.18823f);
99    float y = MAT_ROW_MULT(0.29734f, 0.62736f, 0.07529f);
100    float z = MAT_ROW_MULT(0.02703f, 0.07069f, 0.99134f);
101
102    // The following is the white point in XYZ, so it's simply the row wise addition of the above
103    // matrix.
104    const float xw = 0.5767f + 0.185556f + 0.188212f;
105    const float yw = 0.297361f + 0.627355f + 0.0752847f;
106    const float zw = 0.0270328f + 0.0706879f + 0.991248f;
107
108    // This is the XYZ color point relative to the white point
109    float f[3] = { x / xw, y / yw, z / zw };
110
111    // Conversion from XYZ to LAB taken from
112    // http://en.wikipedia.org/wiki/CIELAB#Forward_transformation
113    for (int i = 0; i < 3; i++) {
114        if (f[i] >= 0.008856f) {
115            f[i] = SkPMetricUtil::get_cube_root(f[i]);
116        } else {
117            f[i] = 7.787f * f[i] + 4.0f / 29.0f;
118        }
119    }
120    lab->l = 116.0f * f[1] - 16.0f;
121    lab->a = 500.0f * (f[0] - f[1]);
122    lab->b = 200.0f * (f[1] - f[2]);
123}
124
125/// Converts a 8888 bitmap to LAB color space and puts it into the output
126static void bitmap_to_cielab(const SkBitmap* bitmap, ImageLAB* outImageLAB) {
127    SkASSERT(bitmap->config() == SkBitmap::kARGB_8888_Config);
128
129    int width = bitmap->width();
130    int height = bitmap->height();
131    SkASSERT(outImageLAB->width == width);
132    SkASSERT(outImageLAB->height == height);
133
134    bitmap->lockPixels();
135    RGB rgb;
136    LAB lab;
137    for (int y = 0; y < height; y++) {
138        unsigned char* row = (unsigned char*)bitmap->getAddr(0, y);
139        for (int x = 0; x < width; x++) {
140            // Perform gamma correction which is assumed to be 2.2
141            rgb.r = SkPMetricUtil::get_gamma(row[x * 4 + 2]);
142            rgb.g = SkPMetricUtil::get_gamma(row[x * 4 + 1]);
143            rgb.b = SkPMetricUtil::get_gamma(row[x * 4 + 0]);
144            adobergb_to_cielab(rgb.r, rgb.g, rgb.b, &lab);
145            outImageLAB->writePixel(x, y, lab);
146        }
147    }
148    bitmap->unlockPixels();
149}
150
151// From Barten SPIE 1989
152static float contrast_sensitivity(float cyclesPerDegree, float luminance) {
153    float a = 440.0f * powf(1.0f + 0.7f / luminance, -0.2f);
154    float b = 0.3f * powf(1.0f + 100.0f / luminance, 0.15f);
155    return a *
156           cyclesPerDegree *
157           expf(-b * cyclesPerDegree) *
158           sqrtf(1.0f + 0.06f * expf(b * cyclesPerDegree));
159}
160
161#if 0
162// We're keeping these around for reference and in case the lookup tables are no longer desired.
163// They are no longer called by any code in this file.
164
165// From Daly 1993
166 static float visual_mask(float contrast) {
167    float x = powf(392.498f * contrast, 0.7f);
168    x = powf(0.0153f * x, 4.0f);
169    return powf(1.0f + x, 0.25f);
170}
171
172// From Ward Larson Siggraph 1997
173static float threshold_vs_intensity(float adaptationLuminance) {
174    float logLum = log10f(adaptationLuminance);
175    float x;
176    if (logLum < -3.94f) {
177        x = -2.86f;
178    } else if (logLum < -1.44f) {
179        x = powf(0.405f * logLum + 1.6f, 2.18) - 2.86f;
180    } else if (logLum < -0.0184f) {
181        x = logLum - 0.395f;
182    } else if (logLum < 1.9f) {
183        x = powf(0.249f * logLum + 0.65f, 2.7f) - 0.72f;
184    } else {
185        x = logLum - 1.255f;
186    }
187    return powf(10.0f, x);
188}
189
190#endif
191
192/// Simply takes the L channel from the input and puts it into the output
193static void lab_to_l(const ImageLAB* imageLAB, ImageL* outImageL) {
194    for (int y = 0; y < imageLAB->height; y++) {
195        for (int x = 0; x < imageLAB->width; x++) {
196            LAB lab;
197            imageLAB->readPixel(x, y, &lab);
198            outImageL->writePixel(x, y, lab.l);
199        }
200    }
201}
202
203/// Convolves an image with the given filter in one direction and saves it to the output image
204static void convolve(const ImageL* imageL, bool vertical, ImageL* outImageL) {
205    SkASSERT(imageL->width == outImageL->width);
206    SkASSERT(imageL->height == outImageL->height);
207
208    const float matrix[] = { 0.05f, 0.25f, 0.4f, 0.25f, 0.05f };
209    const int matrixCount = sizeof(matrix) / sizeof(float);
210    const int radius = matrixCount / 2;
211
212    // Keep track of what rows are being operated on for quick access.
213    float* rowPtrs[matrixCount]; // Because matrixCount is constant, this won't create a VLA
214    for (int y = radius; y < matrixCount; y++) {
215        rowPtrs[y] = imageL->getRow(y - radius);
216    }
217    float* writeRow = outImageL->getRow(0);
218
219    for (int y = 0; y < imageL->height; y++) {
220        for (int x = 0; x < imageL->width; x++) {
221            float lSum = 0.0f;
222            for (int xx = -radius; xx <= radius; xx++) {
223                int nx = x;
224                int ny = y;
225
226                // We mirror at edges so that edge pixels that the filter weighting still makes
227                // sense.
228                if (vertical) {
229                    ny += xx;
230                    if (ny < 0) {
231                        ny = -ny;
232                    }
233                    if (ny >= imageL->height) {
234                        ny = imageL->height + (imageL->height - ny - 1);
235                    }
236                } else {
237                    nx += xx;
238                    if (nx < 0) {
239                        nx = -nx;
240                    }
241                    if (nx >= imageL->width) {
242                        nx = imageL->width + (imageL->width - nx - 1);
243                    }
244                }
245
246                float weight = matrix[xx + radius];
247                lSum += rowPtrs[ny - y + radius][nx] * weight;
248            }
249            writeRow[x] = lSum;
250        }
251        // As we move down, scroll the row pointers down with us
252        for (int y = 0; y < matrixCount - 1; y++)
253        {
254            rowPtrs[y] = rowPtrs[y + 1];
255        }
256        rowPtrs[matrixCount - 1] += imageL->width;
257        writeRow += imageL->width;
258    }
259}
260
261static double pmetric(const ImageLAB* baselineLAB, const ImageLAB* testLAB, SkTDArray<SkIPoint>* poi) {
262    int width = baselineLAB->width;
263    int height = baselineLAB->height;
264    int maxLevels = 0;
265
266    // Calculates how many levels to make by how many times the image can be divided in two
267    int smallerDimension = width < height ? width : height;
268    for ( ; smallerDimension > 1; smallerDimension /= 2) {
269        maxLevels++;
270    }
271
272    const float fov = SK_ScalarPI / 180.0f * 45.0f;
273    float contrastSensitivityMax = contrast_sensitivity(3.248f, 100.0f);
274    float pixelsPerDegree = width / (2.0f * tanf(fov * 0.5f) * 180.0f / SK_ScalarPI);
275
276    ImageL3D baselineL(width, height, maxLevels);
277    ImageL3D testL(width, height, maxLevels);
278    ImageL scratchImageL(width, height);
279    float* cyclesPerDegree = SkNEW_ARRAY(float, maxLevels);
280    float* thresholdFactorFrequency = SkNEW_ARRAY(float, maxLevels - 2);
281    float* contrast = SkNEW_ARRAY(float, maxLevels - 2);
282
283    lab_to_l(baselineLAB, baselineL.getLayer(0));
284    lab_to_l(testLAB, testL.getLayer(0));
285
286    // Compute cpd - Cycles per degree on the pyramid
287    cyclesPerDegree[0] = 0.5f * pixelsPerDegree;
288    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
289        cyclesPerDegree[levelIndex] = cyclesPerDegree[levelIndex - 1] * 0.5f;
290    }
291
292    // Contrast sensitivity is based on image dimensions. Therefore it cannot be statically
293    // generated.
294    float* contrastSensitivityTable = SkNEW_ARRAY(float, maxLevels * 1000);
295    for (int levelIndex = 0; levelIndex < maxLevels; levelIndex++) {
296        for (int csLum = 0; csLum < 1000; csLum++) {
297           contrastSensitivityTable[levelIndex * 1000 + csLum] =
298           contrast_sensitivity(cyclesPerDegree[levelIndex], (float)csLum / 10.0f + 1e-5f);
299       }
300    }
301
302    // Compute G - The convolved lum for the baseline
303    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
304        convolve(baselineL.getLayer(levelIndex - 1), false, &scratchImageL);
305        convolve(&scratchImageL, true, baselineL.getLayer(levelIndex));
306    }
307    for (int levelIndex = 1; levelIndex < maxLevels; levelIndex++) {
308        convolve(testL.getLayer(levelIndex - 1), false, &scratchImageL);
309        convolve(&scratchImageL, true, testL.getLayer(levelIndex));
310    }
311
312    // Compute F_freq - The elevation f
313    for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
314        float cpd = cyclesPerDegree[levelIndex];
315        thresholdFactorFrequency[levelIndex] = contrastSensitivityMax /
316                                               contrast_sensitivity(cpd, 100.0f);
317    }
318
319    int failures = 0;
320    // Calculate F
321    for (int y = 0; y < height; y++) {
322        for (int x = 0; x < width; x++) {
323            float lBaseline;
324            float lTest;
325            baselineL.getLayer(0)->readPixel(x, y, &lBaseline);
326            testL.getLayer(0)->readPixel(x, y, &lTest);
327
328            float avgLBaseline;
329            float avgLTest;
330            baselineL.getLayer(maxLevels - 1)->readPixel(x, y, &avgLBaseline);
331            testL.getLayer(maxLevels - 1)->readPixel(x, y, &avgLTest);
332
333            float lAdapt = 0.5f * (avgLBaseline + avgLTest);
334            if (lAdapt < 1e-5f) {
335                lAdapt = 1e-5f;
336            }
337
338            float contrastSum = 0.0f;
339            for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
340                float baselineL0, baselineL1, baselineL2;
341                float testL0, testL1, testL2;
342                baselineL.getLayer(levelIndex + 0)->readPixel(x, y, &baselineL0);
343                testL.    getLayer(levelIndex + 0)->readPixel(x, y, &testL0);
344                baselineL.getLayer(levelIndex + 1)->readPixel(x, y, &baselineL1);
345                testL.    getLayer(levelIndex + 1)->readPixel(x, y, &testL1);
346                baselineL.getLayer(levelIndex + 2)->readPixel(x, y, &baselineL2);
347                testL.    getLayer(levelIndex + 2)->readPixel(x, y, &testL2);
348
349                float baselineContrast1 = fabsf(baselineL0 - baselineL1);
350                float testContrast1     = fabsf(testL0 - testL1);
351                float numerator = (baselineContrast1 > testContrast1) ?
352                                   baselineContrast1 : testContrast1;
353
354                float baselineContrast2 = fabsf(baselineL2);
355                float testContrast2     = fabsf(testL2);
356                float denominator = (baselineContrast2 > testContrast2) ?
357                                    baselineContrast2 : testContrast2;
358
359                // Avoid divides by close to zero
360                if (denominator < 1e-5f) {
361                    denominator = 1e-5f;
362                }
363                contrast[levelIndex] = numerator / denominator;
364                contrastSum += contrast[levelIndex];
365            }
366
367            if (contrastSum < 1e-5f) {
368                contrastSum = 1e-5f;
369            }
370
371            float F = 0.0f;
372            for (int levelIndex = 0; levelIndex < maxLevels - 2; levelIndex++) {
373                float contrastSensitivity = contrastSensitivityTable[levelIndex * 1000 +
374                                                                     (int)(lAdapt * 10.0)];
375                float mask = SkPMetricUtil::get_visual_mask(contrast[levelIndex] *
376                                                            contrastSensitivity);
377
378                F += contrast[levelIndex] +
379                     thresholdFactorFrequency[levelIndex] * mask / contrastSum;
380            }
381
382            if (F < 1.0f) {
383                F = 1.0f;
384            }
385
386            if (F > 10.0f) {
387                F = 10.0f;
388            }
389
390
391            bool isFailure = false;
392            if (fabsf(lBaseline - lTest) > F * SkPMetricUtil::get_threshold_vs_intensity(lAdapt)) {
393                isFailure = true;
394            } else {
395                LAB baselineColor;
396                LAB testColor;
397                baselineLAB->readPixel(x, y, &baselineColor);
398                testLAB->readPixel(x, y, &testColor);
399                float contrastA = baselineColor.a - testColor.a;
400                float contrastB = baselineColor.b - testColor.b;
401                float colorScale = 1.0f;
402                if (lAdapt < 10.0f) {
403                    colorScale = lAdapt / 10.0f;
404                }
405                colorScale *= colorScale;
406
407                if ((contrastA * contrastA + contrastB * contrastB) * colorScale > F)
408                {
409                    isFailure = true;
410                }
411            }
412
413            if (isFailure) {
414                failures++;
415                poi->push()->set(x, y);
416            }
417        }
418    }
419
420    SkDELETE_ARRAY(cyclesPerDegree);
421    SkDELETE_ARRAY(contrast);
422    SkDELETE_ARRAY(thresholdFactorFrequency);
423    SkDELETE_ARRAY(contrastSensitivityTable);
424    return 1.0 - (double)failures / (width * height);
425}
426
427const char* SkPMetric::getName() {
428    return "perceptual";
429}
430
431int SkPMetric::queueDiff(SkBitmap* baseline, SkBitmap* test) {
432    double startTime = get_seconds();
433    int diffID = fQueuedDiffs.count();
434    QueuedDiff& diff = fQueuedDiffs.push_back();
435    diff.result = 0.0;
436
437    // Ensure the images are comparable
438    if (baseline->width() != test->width() || baseline->height() != test->height() ||
439                    baseline->width() <= 0 || baseline->height() <= 0) {
440        diff.finished = true;
441        return diffID;
442    }
443
444    ImageLAB baselineLAB(baseline->width(), baseline->height());
445    ImageLAB testLAB(baseline->width(), baseline->height());
446
447    bitmap_to_cielab(baseline, &baselineLAB);
448    bitmap_to_cielab(test, &testLAB);
449
450    diff.result = pmetric(&baselineLAB, &testLAB, &diff.poi);
451
452    SkDebugf("Time: %f\n", (get_seconds() - startTime));
453
454    return diffID;
455}
456
457
458void SkPMetric::deleteDiff(int id) {
459
460}
461
462bool SkPMetric::isFinished(int id) {
463    return fQueuedDiffs[id].finished;
464}
465
466double SkPMetric::getResult(int id) {
467    return fQueuedDiffs[id].result;
468}
469
470int SkPMetric::getPointsOfInterestCount(int id) {
471    return fQueuedDiffs[id].poi.count();
472}
473
474SkIPoint* SkPMetric::getPointsOfInterest(int id) {
475    return fQueuedDiffs[id].poi.begin();
476}
477