test/pdiff/pdiff.c

   1 /*
   2   Metric
   3   Copyright (C) 2006 Yangli Hector Yee
   4
   5   This program is free software; you can redistribute it and/or modify it under the terms of the
   6   GNU General Public License as published by the Free Software Foundation; either version 2 of the License,
   7   or (at your option) any later version.
   8
   9   This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
  10   without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  11   See the GNU General Public License for more details.
  12
  13   You should have received a copy of the GNU General Public License along with this program;
  14   if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA
  15 */
  16
  17 #define _GNU_SOURCE
  18
  19 #if HAVE_CONFIG_H
  20 #include "config.h"
  21 #endif
  22
  23 #include "lpyramid.h"
  24 #include <math.h>
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27
  28 #if   HAVE_STDINT_H
  29 # include <stdint.h>
  30 #elif HAVE_INTTYPES_H
  31 # include <inttypes.h>
  32 #elif HAVE_SYS_INT_TYPES_H
  33 # include <sys/int_types.h>
  34 #elif defined(_MSC_VER)
  35   typedef __int8 int8_t;
  36   typedef unsigned __int8 uint8_t;
  37   typedef __int16 int16_t;
  38   typedef unsigned __int16 uint16_t;
  39   typedef __int32 int32_t;
  40   typedef unsigned __int32 uint32_t;
  41   typedef __int64 int64_t;
  42   typedef unsigned __int64 uint64_t;
  43 # ifndef HAVE_UINT64_T
  44 #  define HAVE_UINT64_T 1
  45 # endif
  46 # ifndef INT16_MIN
  47 #  define INT16_MIN     (-32767-1)
  48 # endif
  49 # ifndef INT16_MAX
  50 #  define INT16_MAX     (32767)
  51 # endif
  52 # ifndef UINT16_MAX
  53 #  define UINT16_MAX    (65535)
  54 # endif
  55 #else
  56 #error Cannot find definitions for fixed-width integral types (uint8_t, uint32_t, etc.)
  57 #endif
  58
  59 #include "pdiff.h"
  60
  61 #ifndef M_PI
  62 #define M_PI 3.14159265f
  63 #endif
  64
  65 #ifndef __USE_ISOC99
  66 #define expf    exp
  67 #define powf    pow
  68 #define fabsf   fabs
  69 #define sqrtf   sqrt
  70 #define log10f  log10
  71 #endif
  72
  73 /*
  74  * Given the adaptation luminance, this function returns the
  75  * threshold of visibility in cd per m^2
  76  * TVI means Threshold vs Intensity function
  77  * This version comes from Ward Larson Siggraph 1997
  78  */
  79 static float
  80 tvi (float adaptation_luminance)
  81 {
  82     /* returns the threshold luminance given the adaptation luminance
  83        units are candelas per meter squared
  84     */
  85     float log_a, r, result;
  86     log_a = log10f(adaptation_luminance);
  87
  88     if (log_a < -3.94f) {
  89         r = -2.86f;
  90     } else if (log_a < -1.44f) {
  91         r = powf(0.405f * log_a + 1.6f , 2.18f) - 2.86f;
  92     } else if (log_a < -0.0184f) {
  93         r = log_a - 0.395f;
  94     } else if (log_a < 1.9f) {
  95         r = powf(0.249f * log_a + 0.65f, 2.7f) - 0.72f;
  96     } else {
  97         r = log_a - 1.255f;
  98     }
  99
 100     result = powf(10.0f , r);
 101
 102     return result;
 103 }
 104
 105 /* computes the contrast sensitivity function (Barten SPIE 1989)
 106  * given the cycles per degree (cpd) and luminance (lum)
 107  */
 108 static float
 109 csf (float cpd, float lum)
 110 {
 111     float a, b, result;
 112
 113     a = 440.0f * powf((1.0f + 0.7f / lum), -0.2f);
 114     b = 0.3f * powf((1.0f + 100.0f / lum), 0.15f);
 115
 116     result = a * cpd * expf(-b * cpd) * sqrtf(1.0f + 0.06f * expf(b * cpd));
 117
 118     return result;
 119 }
 120
 121 /*
 122  * Visual Masking Function
 123  * from Daly 1993
 124  */
 125 static float
 126 mask (float contrast)
 127 {
 128     float a, b, result;
 129     a = powf(392.498f * contrast,  0.7f);
 130     b = powf(0.0153f * a, 4.0f);
 131     result = powf(1.0f + b, 0.25f);
 132
 133     return result;
 134 }
 135
 136 /* convert Adobe RGB (1998) with reference white D65 to XYZ */
 137 static void
 138 AdobeRGBToXYZ (float r, float g, float b, float *x, float *y, float *z)
 139 {
 140     /* matrix is from http://www.brucelindbloom.com/ */
 141     *x = r * 0.576700f + g * 0.185556f + b * 0.188212f;
 142     *y = r * 0.297361f + g * 0.627355f + b * 0.0752847f;
 143     *z = r * 0.0270328f + g * 0.0706879f + b * 0.991248f;
 144 }
 145
 146 static void
 147 XYZToLAB (float x, float y, float z, float *L, float *A, float *B)
 148 {
 149     static float xw = -1;
 150     static float yw;
 151     static float zw;
 152     const float epsilon  = 216.0f / 24389.0f;
 153     const float kappa = 24389.0f / 27.0f;
 154     float f[3];
 155     float r[3];
 156     int i;
 157
 158     /* reference white */
 159     if (xw < 0) {
 160         AdobeRGBToXYZ(1, 1, 1, &xw, &yw, &zw);
 161     }
 162     r[0] = x / xw;
 163     r[1] = y / yw;
 164     r[2] = z / zw;
 165     for (i = 0; i < 3; i++) {
 166         if (r[i] > epsilon) {
 167             f[i] = powf(r[i], 1.0f / 3.0f);
 168         } else {
 169             f[i] = (kappa * r[i] + 16.0f) / 116.0f;
 170         }
 171     }
 172     *L = 116.0f * f[1] - 16.0f;
 173     *A = 500.0f * (f[0] - f[1]);
 174     *B = 200.0f * (f[1] - f[2]);
 175 }
 176
 177 static uint32_t
 178 _get_pixel (const uint32_t *data, int i)
 179 {
 180     return data[i];
 181 }
 182
 183 static unsigned char
 184 _get_red (const uint32_t *data, int i)
 185 {
 186     uint32_t pixel;
 187     uint8_t alpha;
 188
 189     pixel = _get_pixel (data, i);
 190     alpha = (pixel & 0xff000000) >> 24;
 191     if (alpha == 0)
 192         return 0;
 193     else
 194         return (((pixel & 0x00ff0000) >> 16) * 255 + alpha / 2) / alpha;
 195 }
 196
 197 static unsigned char
 198 _get_green (const uint32_t *data, int i)
 199 {
 200     uint32_t pixel;
 201     uint8_t alpha;
 202
 203     pixel = _get_pixel (data, i);
 204     alpha = (pixel & 0xff000000) >> 24;
 205     if (alpha == 0)
 206         return 0;
 207     else
 208         return (((pixel & 0x0000ff00) >> 8) * 255 + alpha / 2) / alpha;
 209 }
 210
 211 static unsigned char
 212 _get_blue (const uint32_t *data, int i)
 213 {
 214     uint32_t pixel;
 215     uint8_t alpha;
 216
 217     pixel = _get_pixel (data, i);
 218     alpha = (pixel & 0xff000000) >> 24;
 219     if (alpha == 0)
 220         return 0;
 221     else
 222         return (((pixel & 0x000000ff) >> 0) * 255 + alpha / 2) / alpha;
 223 }
 224
 225 static void *
 226 xmalloc (size_t size)
 227 {
 228     void *buf;
 229
 230     buf = malloc (size);
 231     if (buf == NULL) {
 232         fprintf (stderr, "Out of memory.\n");
 233         exit (1);
 234     }
 235
 236     return buf;
 237 }
 238
 239 int
 240 pdiff_compare (cairo_surface_t *surface_a,
 241                cairo_surface_t *surface_b,
 242                double gamma,
 243                double luminance,
 244                double field_of_view)
 245 {
 246     unsigned int dim = (cairo_image_surface_get_width (surface_a)
 247                         * cairo_image_surface_get_height (surface_a));
 248     unsigned int i;
 249
 250     /* assuming colorspaces are in Adobe RGB (1998) convert to XYZ */
 251     float *aX;
 252     float *aY;
 253     float *aZ;
 254     float *bX;
 255     float *bY;
 256     float *bZ;
 257     float *aLum;
 258     float *bLum;
 259
 260     float *aA;
 261     float *bA;
 262     float *aB;
 263     float *bB;
 264
 265     unsigned int x, y, w, h;
 266
 267     lpyramid_t *la, *lb;
 268
 269     float num_one_degree_pixels, pixels_per_degree, num_pixels;
 270     unsigned int adaptation_level;
 271
 272     float cpd[MAX_PYR_LEVELS];
 273     float F_freq[MAX_PYR_LEVELS - 2];
 274     float csf_max;
 275     const uint32_t *data_a, *data_b;
 276
 277     unsigned int pixels_failed;
 278
 279     w = cairo_image_surface_get_width (surface_a);
 280     h = cairo_image_surface_get_height (surface_a);
 281     if (w < 3 || h < 3) /* too small for the Laplacian convolution */
 282         return -1;
 283
 284     aX = xmalloc (dim * sizeof (float));
 285     aY = xmalloc (dim * sizeof (float));
 286     aZ = xmalloc (dim * sizeof (float));
 287     bX = xmalloc (dim * sizeof (float));
 288     bY = xmalloc (dim * sizeof (float));
 289     bZ = xmalloc (dim * sizeof (float));
 290     aLum = xmalloc (dim * sizeof (float));
 291     bLum = xmalloc (dim * sizeof (float));
 292
 293     aA = xmalloc (dim * sizeof (float));
 294     bA = xmalloc (dim * sizeof (float));
 295     aB = xmalloc (dim * sizeof (float));
 296     bB = xmalloc (dim * sizeof (float));
 297
 298     data_a = (uint32_t *) cairo_image_surface_get_data (surface_a);
 299     data_b = (uint32_t *) cairo_image_surface_get_data (surface_b);
 300     for (y = 0; y < h; y++) {
 301         for (x = 0; x < w; x++) {
 302             float r, g, b, l;
 303             i = x + y * w;
 304             r = powf(_get_red (data_a, i) / 255.0f, gamma);
 305             g = powf(_get_green (data_a, i) / 255.0f, gamma);
 306             b = powf(_get_blue (data_a, i) / 255.0f, gamma);
 307
 308             AdobeRGBToXYZ(r,g,b,&aX[i],&aY[i],&aZ[i]);
 309             XYZToLAB(aX[i], aY[i], aZ[i], &l, &aA[i], &aB[i]);
 310             r = powf(_get_red (data_b, i) / 255.0f, gamma);
 311             g = powf(_get_green (data_b, i) / 255.0f, gamma);
 312             b = powf(_get_blue (data_b, i) / 255.0f, gamma);
 313
 314             AdobeRGBToXYZ(r,g,b,&bX[i],&bY[i],&bZ[i]);
 315             XYZToLAB(bX[i], bY[i], bZ[i], &l, &bA[i], &bB[i]);
 316             aLum[i] = aY[i] * luminance;
 317             bLum[i] = bY[i] * luminance;
 318         }
 319     }
 320
 321     la = lpyramid_create (aLum, w, h);
 322     lb = lpyramid_create (bLum, w, h);
 323
 324     num_one_degree_pixels = (float) (2 * tan(field_of_view * 0.5 * M_PI / 180) * 180 / M_PI);
 325     pixels_per_degree = w / num_one_degree_pixels;
 326
 327     num_pixels = 1;
 328     adaptation_level = 0;
 329     for (i = 0; i < MAX_PYR_LEVELS; i++) {
 330         adaptation_level = i;
 331         if (num_pixels > num_one_degree_pixels) break;
 332         num_pixels *= 2;
 333     }
 334
 335     cpd[0] = 0.5f * pixels_per_degree;
 336     for (i = 1; i < MAX_PYR_LEVELS; i++) cpd[i] = 0.5f * cpd[i - 1];
 337     csf_max = csf(3.248f, 100.0f);
 338
 339     for (i = 0; i < MAX_PYR_LEVELS - 2; i++) F_freq[i] = csf_max / csf( cpd[i], 100.0f);
 340
 341     pixels_failed = 0;
 342     for (y = 0; y < h; y++) {
 343         for (x = 0; x < w; x++) {
 344             int index = x + y * w;
 345             float contrast[MAX_PYR_LEVELS - 2];
 346             float F_mask[MAX_PYR_LEVELS - 2];
 347             float factor;
 348             float delta;
 349             float adapt;
 350             bool pass;
 351             float sum_contrast = 0;
 352             for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {
 353                 float n1 = fabsf(lpyramid_get_value (la,x,y,i) - lpyramid_get_value (la,x,y,i + 1));
 354                 float n2 = fabsf(lpyramid_get_value (lb,x,y,i) - lpyramid_get_value (lb,x,y,i + 1));
 355                 float numerator = (n1 > n2) ? n1 : n2;
 356                 float d1 = fabsf(lpyramid_get_value(la,x,y,i+2));
 357                 float d2 = fabsf(lpyramid_get_value(lb,x,y,i+2));
 358                 float denominator = (d1 > d2) ? d1 : d2;
 359                 if (denominator < 1e-5f) denominator = 1e-5f;
 360                 contrast[i] = numerator / denominator;
 361                 sum_contrast += contrast[i];
 362             }
 363             if (sum_contrast < 1e-5) sum_contrast = 1e-5f;
 364             adapt = lpyramid_get_value(la,x,y,adaptation_level) + lpyramid_get_value(lb,x,y,adaptation_level);
 365             adapt *= 0.5f;
 366             if (adapt < 1e-5) adapt = 1e-5f;
 367             for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {
 368                 F_mask[i] = mask(contrast[i] * csf(cpd[i], adapt));
 369             }
 370             factor = 0;
 371             for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {
 372                 factor += contrast[i] * F_freq[i] * F_mask[i] / sum_contrast;
 373             }
 374             if (factor < 1) factor = 1;
 375             if (factor > 10) factor = 10;
 376             delta = fabsf(lpyramid_get_value(la,x,y,0) - lpyramid_get_value(lb,x,y,0));
 377             pass = true;
 378             /* pure luminance test */
 379             if (delta > factor * tvi(adapt)) {
 380                 pass = false;
 381             } else {
 382                 /* CIE delta E test with modifications */
 383                 float color_scale = 1.0f;
 384                 float da = aA[index] - bA[index];
 385                 float db = aB[index] - bB[index];
 386                 float delta_e;
 387                 /* ramp down the color test in scotopic regions */
 388                 if (adapt < 10.0f) {
 389                     color_scale = 1.0f - (10.0f - color_scale) / 10.0f;
 390                     color_scale = color_scale * color_scale;
 391                 }
 392                 da = da * da;
 393                 db = db * db;
 394                 delta_e = (da + db) * color_scale;
 395                 if (delta_e > factor) {
 396                     pass = false;
 397                 }
 398             }
 399             if (!pass)
 400                 pixels_failed++;
 401         }
 402     }
 403
 404     free (aX);
 405     free (aY);
 406     free (aZ);
 407     free (bX);
 408     free (bY);
 409     free (bZ);
 410     free (aLum);
 411     free (bLum);
 412     lpyramid_destroy (la);
 413     lpyramid_destroy (lb);
 414     free (aA);
 415     free (bA);
 416     free (aB);
 417     free (bB);
 418
 419     return pixels_failed;
 420 }