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