1 /* 2 Metric 3 Copyright (C) 2006-2011 Yangli Hector Yee 4 Copyright (C) 2011-2016 Steven Myint, Jeff Terrace 5 6 This program is free software; you can redistribute it and/or modify it under 7 the terms of the GNU General Public License as published by the Free Software 8 Foundation; either version 2 of the License, or (at your option) any later 9 version. 10 11 This program is distributed in the hope that it will be useful, but WITHOUT ANY 12 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 13 PARTICULAR PURPOSE. See the GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License along with 16 this program; if not, write to the Free Software Foundation, Inc., 59 Temple 17 Place, Suite 330, Boston, MA 02111-1307 USA 18 */ 19 20 #include "metric.h" 21 22 #include "lpyramid.h" 23 #include "rgba_image.h" 24 25 #include <ciso646> 26 #include <cmath> 27 #include <cstddef> 28 #include <iostream> 29 #include <vector> 30 #include <algorithm> 31 32 33 namespace pdiff 34 { 35 #if _MSC_VER <= 1800 36 static const auto pi = 3.14159265f; 37 38 to_radians(const float degrees)39 static float to_radians(const float degrees) // LCOV_EXCL_LINE 40 { 41 return degrees * pi / 180.f; // LCOV_EXCL_LINE 42 } 43 44 to_degrees(const float radians)45 static float to_degrees(const float radians) // LCOV_EXCL_LINE 46 { 47 return radians * 180.f / pi; // LCOV_EXCL_LINE 48 } 49 #else 50 constexpr auto pi = 3.14159265f; 51 52 53 constexpr float to_radians(const float degrees) // LCOV_EXCL_LINE 54 { 55 return degrees * pi / 180.f; // LCOV_EXCL_LINE 56 } 57 58 59 constexpr float to_degrees(const float radians) // LCOV_EXCL_LINE 60 { 61 return radians * 180.f / pi; // LCOV_EXCL_LINE 62 } 63 #endif 64 65 // Given the adaptation luminance, this function returns the 66 // threshold of visibility in cd per m^2. 67 // 68 // TVI means Threshold vs Intensity function. 69 // This version comes from Ward Larson Siggraph 1997. 70 // 71 // Returns the threshold luminance given the adaptation luminance. 72 // Units are candelas per meter squared. tvi(const float adaptation_luminance)73 static float tvi(const float adaptation_luminance) 74 { 75 const auto log_a = log10f(adaptation_luminance); 76 77 float r; 78 if (log_a < -3.94f) 79 { 80 r = -2.86f; 81 } 82 else if (log_a < -1.44f) 83 { 84 r = powf(0.405f * log_a + 1.6f, 2.18f) - 2.86f; 85 } 86 else if (log_a < -0.0184f) 87 { 88 r = log_a - 0.395f; 89 } 90 else if (log_a < 1.9f) 91 { 92 r = powf(0.249f * log_a + 0.65f, 2.7f) - 0.72f; 93 } 94 else 95 { 96 r = log_a - 1.255f; 97 } 98 99 return powf(10.0f, r); 100 } 101 102 103 // computes the contrast sensitivity function (Barten SPIE 1989) 104 // given the cycles per degree (cpd) and luminance (lum) csf(const float cpd,const float lum)105 static float csf(const float cpd, const float lum) 106 { 107 const auto a = 440.f * powf((1.f + 0.7f / lum), -0.2f); 108 const auto b = 0.3f * powf((1.0f + 100.0f / lum), 0.15f); 109 110 return a * cpd * expf(-b * cpd) * sqrtf(1.0f + 0.06f * expf(b * cpd)); 111 } 112 113 114 /* 115 * Visual Masking Function 116 * from Daly 1993 117 */ mask(const float contrast)118 static float mask(const float contrast) 119 { 120 const auto a = powf(392.498f * contrast, 0.7f); 121 const auto b = powf(0.0153f * a, 4.f); 122 return powf(1.0f + b, 0.25f); 123 } 124 125 126 // convert Adobe RGB (1998) with reference white D65 to XYZ adobe_rgb_to_xyz(const float r,const float g,const float b,float & x,float & y,float & z)127 static void adobe_rgb_to_xyz(const float r, const float g, const float b, 128 float &x, float &y, float &z) 129 { 130 // matrix is from http://www.brucelindbloom.com/ 131 x = r * 0.576700f + g * 0.185556f + b * 0.188212f; 132 y = r * 0.297361f + g * 0.627355f + b * 0.0752847f; 133 z = r * 0.0270328f + g * 0.0706879f + b * 0.991248f; 134 } 135 136 137 struct White 138 { Whitepdiff::White139 White() 140 { 141 adobe_rgb_to_xyz(1.f, 1.f, 1.f, x, y, z); 142 } 143 144 float x; 145 float y; 146 float z; 147 }; 148 149 150 static const White global_white; 151 152 xyz_to_lab(const float x,const float y,const float z,float & l,float & a,float & b)153 static void xyz_to_lab(const float x, const float y, const float z, 154 float &l, float &a, float &b) 155 { 156 const float epsilon = 216.0f / 24389.0f; 157 const float kappa = 24389.0f / 27.0f; 158 const float r[] = { 159 x / global_white.x, 160 y / global_white.y, 161 z / global_white.z 162 }; 163 float f[3]; 164 for (auto i = 0u; i < 3; i++) 165 { 166 if (r[i] > epsilon) 167 { 168 f[i] = powf(r[i], 1.0f / 3.0f); 169 } 170 else 171 { 172 f[i] = (kappa * r[i] + 16.0f) / 116.0f; 173 } 174 } 175 l = 116.0f * f[1] - 16.0f; 176 a = 500.0f * (f[0] - f[1]); 177 b = 200.0f * (f[1] - f[2]); 178 } 179 180 adaptation(const float num_one_degree_pixels)181 static unsigned int adaptation(const float num_one_degree_pixels) 182 { 183 auto num_pixels = 1.f; 184 auto adaptation_level = 0u; 185 for (auto i = 0u; i < MAX_PYR_LEVELS; i++) 186 { 187 adaptation_level = i; 188 if (num_pixels > num_one_degree_pixels) 189 { 190 break; 191 } 192 num_pixels *= 2; 193 } 194 return adaptation_level; // LCOV_EXCL_LINE 195 } 196 197 PerceptualDiffParameters()198 PerceptualDiffParameters::PerceptualDiffParameters() 199 { 200 luminance_only = false; 201 field_of_view = 45.0f; 202 gamma = 2.2f; 203 threshold_pixels = 100; 204 luminance = 100.0f; 205 color_factor = 1.0f; 206 down_sample = 0; 207 } 208 209 yee_compare(const PerceptualDiffParameters & args,const RGBAImage & image_a,const RGBAImage & image_b,size_t * output_num_pixels_failed,float * output_error_sum,std::string * output_reason,RGBAImage * output_image_difference,std::ostream * output_verbose)210 bool yee_compare(const PerceptualDiffParameters &args, 211 const RGBAImage &image_a, 212 const RGBAImage &image_b, 213 size_t *output_num_pixels_failed, 214 float *output_error_sum, 215 std::string *output_reason, 216 RGBAImage *output_image_difference, 217 std::ostream *output_verbose) 218 { 219 if ((image_a.get_width() != image_b.get_width()) or 220 (image_a.get_height() != image_b.get_height())) 221 { 222 if (output_reason) 223 { 224 *output_reason = "Image dimensions do not match\n"; 225 } 226 return false; 227 } 228 229 const auto w = image_a.get_width(); 230 const auto h = image_a.get_height(); 231 const auto dim = w * h; 232 233 auto identical = true; 234 for (auto i = 0u; i < dim; i++) 235 { 236 if (image_a.get(i) != image_b.get(i)) 237 { 238 identical = false; 239 break; 240 } 241 } 242 if (identical) 243 { 244 if (output_reason) 245 { 246 *output_reason = "Images are binary identical\n"; 247 } 248 return true; 249 } 250 251 // Assuming colorspaces are in Adobe RGB (1998) convert to XYZ. 252 std::vector<float> a_lum(dim); 253 std::vector<float> b_lum(dim); 254 255 std::vector<float> a_a(dim); 256 std::vector<float> b_a(dim); 257 std::vector<float> a_b(dim); 258 std::vector<float> b_b(dim); 259 260 if (output_verbose) 261 { 262 *output_verbose << "Converting RGB to XYZ\n"; 263 } 264 265 const auto gamma = args.gamma; 266 const auto luminance = args.luminance; 267 268 #pragma omp parallel for shared(args, a_lum, b_lum, a_a, a_b, b_a, b_b) 269 for (auto y = 0; y < static_cast<ptrdiff_t>(h); y++) 270 { 271 for (auto x = 0u; x < w; x++) 272 { 273 const auto i = x + y * w; 274 275 // perceptualdiff used to use premultiplied alphas when loading 276 // the image. This is no longer the case since the switch to 277 // FreeImage. We need to do the multiplication here now. As was 278 // the case with premultiplied alphas, differences in alphas 279 // won't be detected where the color is black. 280 281 const auto a_alpha = image_a.get_alpha(i) / 255.f; 282 283 const auto a_color_r = powf( 284 image_a.get_red(i) / 255.f * a_alpha, 285 gamma); 286 const auto a_color_g = powf( 287 image_a.get_green(i) / 255.f * a_alpha, 288 gamma); 289 const auto a_color_b = powf( 290 image_a.get_blue(i) / 255.f * a_alpha, 291 gamma); 292 293 float a_x; 294 float a_y; 295 float a_z; 296 adobe_rgb_to_xyz(a_color_r, a_color_g, a_color_b, 297 a_x, a_y, a_z); 298 float l; 299 xyz_to_lab(a_x, a_y, a_z, l, a_a[i], a_b[i]); 300 301 const auto b_alpha = image_b.get_alpha(i) / 255.f; 302 303 const auto b_color_r = powf( 304 image_b.get_red(i) / 255.f * b_alpha, 305 gamma); 306 const auto b_color_g = powf( 307 image_b.get_green(i) / 255.f * b_alpha, 308 gamma); 309 const auto b_color_b = powf( 310 image_b.get_blue(i) / 255.f * b_alpha, 311 gamma); 312 313 float b_x; 314 float b_y; 315 float b_z; 316 adobe_rgb_to_xyz(b_color_r, b_color_g, b_color_b, 317 b_x, b_y, b_z); 318 xyz_to_lab(b_x, b_y, b_z, l, b_a[i], b_b[i]); 319 320 a_lum[i] = a_y * luminance; 321 b_lum[i] = b_y * luminance; 322 } 323 } 324 325 if (output_verbose) 326 { 327 *output_verbose << "Constructing Laplacian Pyramids\n"; 328 } 329 330 const LPyramid la(a_lum, w, h); 331 const LPyramid lb(b_lum, w, h); 332 333 const auto num_one_degree_pixels = 334 to_degrees(2 * 335 std::tan(args.field_of_view * to_radians(.5f))); 336 const auto pixels_per_degree = w / num_one_degree_pixels; 337 338 if (output_verbose) 339 { 340 *output_verbose << "Performing test\n"; 341 } 342 343 const auto adaptation_level = adaptation(num_one_degree_pixels); 344 345 float cpd[MAX_PYR_LEVELS]; 346 cpd[0] = 0.5f * pixels_per_degree; 347 for (auto i = 1u; i < MAX_PYR_LEVELS; i++) 348 { 349 cpd[i] = 0.5f * cpd[i - 1]; 350 } 351 const auto csf_max = csf(3.248f, 100.0f); 352 353 static_assert(MAX_PYR_LEVELS > 2, 354 "MAX_PYR_LEVELS must be greater than 2"); 355 356 float f_freq[MAX_PYR_LEVELS - 2]; 357 for (auto i = 0u; i < MAX_PYR_LEVELS - 2; i++) 358 { 359 f_freq[i] = csf_max / csf(cpd[i], 100.0f); 360 } 361 362 auto pixels_failed = 0u; 363 auto error_sum = 0.; 364 365 #pragma omp parallel for reduction(+ : pixels_failed, error_sum) \ 366 shared(args, a_a, a_b, b_a, b_b, cpd, f_freq) 367 for (auto y = 0; y < static_cast<ptrdiff_t>(h); y++) 368 { 369 for (auto x = 0u; x < w; x++) 370 { 371 const auto index = y * w + x; 372 373 const auto adapt = std::max( 374 (la.get_value(x, y, adaptation_level) + 375 lb.get_value(x, y, adaptation_level)) * 0.5f, 376 1e-5f); 377 378 auto sum_contrast = 0.f; 379 auto factor = 0.f; 380 381 for (auto i = 0u; i < MAX_PYR_LEVELS - 2; i++) 382 { 383 const auto n1 = 384 std::abs(la.get_value(x, y, i) - 385 la.get_value(x, y, i + 1)); 386 387 const auto n2 = 388 std::abs(lb.get_value(x, y, i) - 389 lb.get_value(x, y, i + 1)); 390 391 const auto numerator = std::max(n1, n2); 392 const auto d1 = std::abs(la.get_value(x, y, i + 2)); 393 const auto d2 = std::abs(lb.get_value(x, y, i + 2)); 394 const auto denominator = std::max(std::max(d1, d2), 1e-5f); 395 const auto contrast = numerator / denominator; 396 const auto f_mask = mask(contrast * csf(cpd[i], adapt)); 397 factor += contrast * f_freq[i] * f_mask; 398 sum_contrast += contrast; 399 } 400 sum_contrast = std::max(sum_contrast, 1e-5f); 401 factor /= sum_contrast; 402 factor = std::min(std::max(factor, 1.f), 10.f); 403 const auto delta = 404 std::abs(la.get_value(x, y, 0) - lb.get_value(x, y, 0)); 405 error_sum += delta; 406 auto pass = true; 407 408 // Pure luminance test. 409 if (delta > factor * tvi(adapt)) 410 { 411 pass = false; 412 } 413 414 if (not args.luminance_only) 415 { 416 // CIE delta E test with modifications. 417 auto color_scale = args.color_factor; 418 419 // Ramp down the color test in scotopic regions. 420 if (adapt < 10.0f) 421 { 422 // Don't do color test at all. 423 color_scale = 0.0; 424 } 425 426 const auto da = a_a[index] - b_a[index]; 427 const auto db = a_b[index] - b_b[index]; 428 const auto delta_e = (da * da + db * db) * color_scale; 429 error_sum += delta_e; 430 if (delta_e > factor) 431 { 432 pass = false; 433 } 434 } 435 436 if (not pass) 437 { 438 pixels_failed++; 439 if (output_image_difference) 440 { 441 output_image_difference->set(255, 0, 0, 255, index); 442 } 443 } 444 else 445 { 446 if (output_image_difference) 447 { 448 output_image_difference->set(0, 0, 0, 255, index); 449 } 450 } 451 } 452 } 453 454 const auto different = 455 std::to_string(pixels_failed) + " pixels are different\n"; 456 457 const auto passed = pixels_failed < args.threshold_pixels; 458 459 if (output_reason) 460 { 461 if (passed) 462 { 463 *output_reason = 464 "Images are perceptually indistinguishable\n" + different; 465 } 466 else 467 { 468 *output_reason = "Images are visibly different\n" + different; 469 } 470 } 471 472 if (output_num_pixels_failed) 473 { 474 *output_num_pixels_failed = pixels_failed; 475 } 476 477 if (output_error_sum) 478 { 479 *output_error_sum = error_sum; 480 } 481 482 return passed; 483 } 484 } 485