From 832436e2e1b8ada1dff89825716139768511c3c8 Mon Sep 17 00:00:00 2001 From: Jyrki Alakuijala Date: Fri, 22 Jun 2018 14:48:50 +0200 Subject: [PATCH] high-frequency asymmetry for improved ringing detection --- butteraugli/butteraugli.cc | 1313 +++++++++++++++++-------------- butteraugli/butteraugli.h | 39 +- butteraugli/butteraugli_main.cc | 33 +- 3 files changed, 791 insertions(+), 594 deletions(-) diff --git a/butteraugli/butteraugli.cc b/butteraugli/butteraugli.cc index 846c1b0..a9973b1 100755 --- a/butteraugli/butteraugli.cc +++ b/butteraugli/butteraugli.cc @@ -58,6 +58,7 @@ namespace butteraugli { void *CacheAligned::Allocate(const size_t bytes) { + PROFILER_FUNC; char *const allocated = static_cast(malloc(bytes + kCacheLineSize)); if (allocated == nullptr) { return nullptr; @@ -73,6 +74,7 @@ void *CacheAligned::Allocate(const size_t bytes) { } void CacheAligned::Free(void *aligned_pointer) { + PROFILER_FUNC; if (aligned_pointer == nullptr) { return; } @@ -100,6 +102,7 @@ static inline bool IsNan(const double x) { } static inline void CheckImage(const ImageF &image, const char *name) { + PROFILER_FUNC; for (size_t y = 0; y < image.ysize(); ++y) { const float * const BUTTERAUGLI_RESTRICT row = image.Row(y); for (size_t x = 0; x < image.xsize(); ++x) { @@ -135,7 +138,7 @@ static inline void CheckImage(const ImageF &image, const char *name) { // Purpose of kInternalGoodQualityThreshold: // Normalize 'ok' image degradation to 1.0 across different versions of // butteraugli. -static const double kInternalGoodQualityThreshold = 11.698467292807441; +static const double kInternalGoodQualityThreshold = 20.35; static const double kGlobalScale = 1.0 / kInternalGoodQualityThreshold; inline float DotProduct(const float u[3], const float v[3]) { @@ -184,6 +187,7 @@ void ConvolveBorderColumn( ImageF Convolution(const ImageF& in, const std::vector& kernel, const float border_ratio) { + PROFILER_FUNC; ImageF out(in.ysize(), in.xsize()); const int len = kernel.size(); const int offset = kernel.size() / 2; @@ -194,26 +198,30 @@ ImageF Convolution(const ImageF& in, float scale_no_border = 1.0f / weight_no_border; const int border1 = in.xsize() <= offset ? in.xsize() : offset; const int border2 = in.xsize() - offset; - int x = 0; + std::vector scaled_kernel = kernel; + for (int i = 0; i < scaled_kernel.size(); ++i) { + scaled_kernel[i] *= scale_no_border; + } // left border - for (; x < border1; ++x) { + for (int x = 0; x < border1; ++x) { ConvolveBorderColumn(in, kernel, weight_no_border, border_ratio, x, out.Row(x)); } // middle - for (; x < border2; ++x) { - float* const BUTTERAUGLI_RESTRICT row_out = out.Row(x); - for (size_t y = 0; y < in.ysize(); ++y) { - const float* const BUTTERAUGLI_RESTRICT row_in = &in.Row(y)[x - offset]; + for (size_t y = 0; y < in.ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_in = in.Row(y); + for (int x = border1; x < border2; ++x) { + const int d = x - offset; + float* const BUTTERAUGLI_RESTRICT row_out = out.Row(x); float sum = 0.0f; for (int j = 0; j < len; ++j) { - sum += row_in[j] * kernel[j]; + sum += row_in[d + j] * scaled_kernel[j]; } - row_out[y] = sum * scale_no_border; + row_out[y] = sum; } } // right border - for (; x < in.xsize(); ++x) { + for (int x = border2; x < in.xsize(); ++x) { ConvolveBorderColumn(in, kernel, weight_no_border, border_ratio, x, out.Row(x)); } @@ -228,27 +236,6 @@ ImageF Blur(const ImageF& in, float sigma, float border_ratio) { kernel, border_ratio); } -// DoGBlur is an approximate of difference of Gaussians. We use it to -// approximate LoG (Laplacian of Gaussians). -// See: https://en.wikipedia.org/wiki/Difference_of_Gaussians -// For motivation see: -// https://en.wikipedia.org/wiki/Pyramid_(image_processing)#Laplacian_pyramid -ImageF DoGBlur(const ImageF& in, float sigma, float border_ratio) { - ImageF blur1 = Blur(in, sigma, border_ratio); - ImageF blur2 = Blur(in, sigma * 2.0f, border_ratio); - static const float mix = 0.5; - ImageF out(in.xsize(), in.ysize()); - for (size_t y = 0; y < in.ysize(); ++y) { - const float* const BUTTERAUGLI_RESTRICT row1 = blur1.Row(y); - const float* const BUTTERAUGLI_RESTRICT row2 = blur2.Row(y); - float* const BUTTERAUGLI_RESTRICT row_out = out.Row(y); - for (size_t x = 0; x < in.xsize(); ++x) { - row_out[x] = (1.0f + mix) * row1[x] - mix * row2[x]; - } - } - return out; -} - // Clamping linear interpolator. inline double InterpolateClampNegative(const double *array, int size, double ix) { @@ -279,46 +266,6 @@ double GammaMaxArg() { return std::max(out0, std::max(out1, out2)); } -// The input images c0 and c1 include the high frequency component only. -// The output scalar images b0 and b1 include the correlation of Y and -// B component at a Gaussian locality around the respective pixel. -ImageF BlurredBlueCorrelation(const std::vector& uhf, - const std::vector& hf) { - const size_t xsize = uhf[0].xsize(); - const size_t ysize = uhf[0].ysize(); - ImageF yb(xsize, ysize); - ImageF yy(xsize, ysize); - for (size_t y = 0; y < ysize; ++y) { - const float* const BUTTERAUGLI_RESTRICT row_uhf_y = uhf[1].Row(y); - const float* const BUTTERAUGLI_RESTRICT row_uhf_b = uhf[2].Row(y); - const float* const BUTTERAUGLI_RESTRICT row_hf_y = hf[1].Row(y); - const float* const BUTTERAUGLI_RESTRICT row_hf_b = hf[2].Row(y); - float* const BUTTERAUGLI_RESTRICT row_yb = yb.Row(y); - float* const BUTTERAUGLI_RESTRICT row_yy = yy.Row(y); - for (size_t x = 0; x < xsize; ++x) { - const float yval = row_hf_y[x] + row_uhf_y[x]; - const float bval = row_hf_b[x] + row_uhf_b[x]; - row_yb[x] = yval * bval; - row_yy[x] = yval * yval; - } - } - const double kSigma = 8.48596332566; - ImageF yy_blurred = Blur(yy, kSigma, 0.0); - ImageF yb_blurred = Blur(yb, kSigma, 0.0); - for (size_t y = 0; y < ysize; ++y) { - const float* const BUTTERAUGLI_RESTRICT row_uhf_y = uhf[1].Row(y); - const float* const BUTTERAUGLI_RESTRICT row_hf_y = hf[1].Row(y); - const float* const BUTTERAUGLI_RESTRICT row_yy = yy_blurred.Row(y); - float* const BUTTERAUGLI_RESTRICT row_yb = yb_blurred.Row(y); - for (size_t x = 0; x < xsize; ++x) { - static const float epsilon = 20.0101389159; - const float yval = row_hf_y[x] + row_uhf_y[x]; - row_yb[x] *= yval / (row_yy[x] + epsilon); - } - } - return yb_blurred; -} - double SimpleGamma(double v) { static const double kGamma = 0.372322653176; static const double limit = 37.8000499603; @@ -382,10 +329,10 @@ std::vector OpsinDynamicsImage(const std::vector& rgb) { PROFILER_FUNC; std::vector xyb(3); std::vector blurred(3); - const double kSigma = 1.44316781537; + const double kSigma = 1.2; for (int i = 0; i < 3; ++i) { xyb[i] = ImageF(rgb[i].xsize(), rgb[i].ysize()); - blurred[i] = Blur(rgb[i], kSigma, 0.0f); + blurred[i] = Blur(rgb[i], kSigma, 0.0); } for (size_t y = 0; y < rgb[0].ysize(); ++y) { const float* const BUTTERAUGLI_RESTRICT row_r = rgb[0].Row(y); @@ -432,33 +379,6 @@ static BUTTERAUGLI_INLINE float AmplifyRangeAroundZero(float w, float x) { return x > w ? x + w : x < -w ? x - w : 2.0f * x; } -std::vector ModifyRangeAroundZero(const double warray[2], - const std::vector& in) { - std::vector out; - for (int k = 0; k < 3; ++k) { - ImageF plane(in[k].xsize(), in[k].ysize()); - for (int y = 0; y < plane.ysize(); ++y) { - auto row_in = in[k].Row(y); - auto row_out = plane.Row(y); - if (k == 2) { - memcpy(row_out, row_in, plane.xsize() * sizeof(row_out[0])); - } else if (warray[k] >= 0) { - const double w = warray[k]; - for (int x = 0; x < plane.xsize(); ++x) { - row_out[x] = RemoveRangeAroundZero(w, row_in[x]); - } - } else { - const double w = -warray[k]; - for (int x = 0; x < plane.xsize(); ++x) { - row_out[x] = AmplifyRangeAroundZero(w, row_in[x]); - } - } - } - out.emplace_back(std::move(plane)); - } - return out; -} - // XybLowFreqToVals converts from low-frequency XYB space to the 'vals' space. // Vals space can be converted to L2-norm space (Euclidean and normalized) // through visual masking. @@ -467,10 +387,10 @@ BUTTERAUGLI_INLINE void XybLowFreqToVals(const V &x, const V &y, const V &b_arg, V *BUTTERAUGLI_RESTRICT valx, V *BUTTERAUGLI_RESTRICT valy, V *BUTTERAUGLI_RESTRICT valb) { - static const double xmuli = 5.55938080599; - static const double ymuli = 4.58944186612; - static const double bmuli = 11.2394147993; - static const double y_to_b_muli = -0.634050875917; + static const double xmuli = 5.57547552483; + static const double ymuli = 1.20828034498; + static const double bmuli = 6.08319517575; + static const double y_to_b_muli = -0.628811683685; const V xmul(xmuli); const V ymul(ymuli); @@ -482,13 +402,12 @@ BUTTERAUGLI_INLINE void XybLowFreqToVals(const V &x, const V &y, const V &b_arg, *valy = y * ymul; } -static ImageF SuppressHfInBrightAreas(size_t xsize, size_t ysize, - const ImageF& hf, - const ImageF& brightness) { +static ImageF SuppressInBrightAreas(size_t xsize, size_t ysize, + double mul, double mul2, double reg, + const ImageF& hf, + const ImageF& brightness) { + PROFILER_FUNC; ImageF inew(xsize, ysize); - static const float mul = 1.12879309857; - static const float mul2 = 2.27308648104; - static const float reg = 2000 * mul2; for (size_t y = 0; y < ysize; ++y) { const float* const rowhf = hf.Row(y); const float* const rowbr = brightness.Row(y); @@ -503,8 +422,35 @@ static ImageF SuppressHfInBrightAreas(size_t xsize, size_t ysize, } +static float SuppressHfInBrightAreas(float hf, float brightness, + float mul, float reg) { + float scaler = mul * reg / (reg + brightness); + return scaler * hf; +} + +static float SuppressUhfInBrightAreas(float hf, float brightness, + float mul, float reg) { + float scaler = mul * reg / (reg + brightness); + return scaler * hf; +} + +static float MaximumClamp(float v, float maxval) { + static const double kMul = 0.688059627878; + if (v >= maxval) { + v -= maxval; + v *= kMul; + v += maxval; + } else if (v < -maxval) { + v += maxval; + v *= kMul; + v -= maxval; + } + return v; +} + static ImageF MaximumClamping(size_t xsize, size_t ysize, const ImageF& ix, double yw) { + static const double kMul = 0.688059627878; ImageF inew(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { const float* const rowx = ix.Row(y); @@ -513,11 +459,11 @@ static ImageF MaximumClamping(size_t xsize, size_t ysize, const ImageF& ix, double v = rowx[x]; if (v >= yw) { v -= yw; - v *= 0.7; + v *= kMul; v += yw; } else if (v < -yw) { v += yw; - v *= 0.7; + v *= kMul; v -= yw; } rownew[x] = v; @@ -526,22 +472,20 @@ static ImageF MaximumClamping(size_t xsize, size_t ysize, const ImageF& ix, return inew; } -double Suppress(double x, double y) { - static const double yw = 16.1797443814; - static const double s = 0.512720106089; - const double scaler = s + (yw * (1.0 - s)) / (yw + y * y); - return scaler * x; -} - static ImageF SuppressXByY(size_t xsize, size_t ysize, - const ImageF& ix, const ImageF& iy, const double w) { + const ImageF& ix, const ImageF& iy, + const double yw) { + static const double s = 0.745954517135; ImageF inew(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { const float* const rowx = ix.Row(y); const float* const rowy = iy.Row(y); float* const rownew = inew.Row(y); for (size_t x = 0; x < xsize; ++x) { - rownew[x] = Suppress(rowx[x], w * rowy[x]); + const double xval = rowx[x]; + const double yval = rowy[x]; + const double scaler = s + (yw * (1.0 - s)) / (yw + yval * yval); + rownew[x] = scaler * xval; } } return inew; @@ -552,14 +496,23 @@ static void SeparateFrequencies( const std::vector& xyb, PsychoImage &ps) { PROFILER_FUNC; - ps.lf.resize(3); - ps.mf.resize(3); - ps.hf.resize(3); - ps.uhf.resize(3); + ps.lf.resize(3); // XYB + ps.mf.resize(3); // XYB + ps.hf.resize(2); // XY + ps.uhf.resize(2); // XY + // Extract lf ... + static const double kSigmaLf = 7.46953768697; + static const double kSigmaHf = 3.734768843485; + static const double kSigmaUhf = 1.8673844217425; + // At borders we move some more of the energy to the high frequency + // parts, because there can be unfortunate continuations in tiling + // background color etc. So we want to represent the borders with + // some more accuracy. + static double border_lf = -0.00457628248637; + static double border_mf = -0.271277366628; + static double border_hf = 0.147068973249; for (int i = 0; i < 3; ++i) { - // Extract lf ... - static const double kSigmaLf = 7.41525493374; - ps.lf[i] = DoGBlur(xyb[i], kSigmaLf, 0.0f); + ps.lf[i] = Blur(xyb[i], kSigmaLf, border_lf); // ... and keep everything else in mf. ps.mf[i] = ImageF(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { @@ -567,52 +520,96 @@ static void SeparateFrequencies( ps.mf[i].Row(y)[x] = xyb[i].Row(y)[x] - ps.lf[i].Row(y)[x]; } } + if (i == 2) { + ps.mf[i] = Blur(ps.mf[i], kSigmaHf, border_mf); + break; + } // Divide mf into mf and hf. - static const double kSigmaHf = 0.5 * kSigmaLf; ps.hf[i] = ImageF(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[i].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[i].Row(y); for (size_t x = 0; x < xsize; ++x) { - ps.hf[i].Row(y)[x] = ps.mf[i].Row(y)[x]; + row_hf[x] = row_mf[x]; } } - ps.mf[i] = DoGBlur(ps.mf[i], kSigmaHf, 0.0f); - for (size_t y = 0; y < ysize; ++y) { - for (size_t x = 0; x < xsize; ++x) { - ps.hf[i].Row(y)[x] -= ps.mf[i].Row(y)[x]; + ps.mf[i] = Blur(ps.mf[i], kSigmaHf, border_mf); + static const double w0 = 0.120079806822; + static const double w1 = 0.03430529365; + if (i == 0) { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[0].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[0].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_hf[x] -= row_mf[x]; + row_mf[x] = RemoveRangeAroundZero(w0, row_mf[x]); + } + } + } else { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[1].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_hf[x] -= row_mf[x]; + row_mf[x] = AmplifyRangeAroundZero(w1, row_mf[x]); + } } } + } + // Suppress red-green by intensity change in the high freq channels. + static const double suppress = 2.96534974403; + ps.hf[0] = SuppressXByY(xsize, ysize, ps.hf[0], ps.hf[1], suppress); + + for (int i = 0; i < 2; ++i) { // Divide hf into hf and uhf. - static const double kSigmaUhf = 0.5 * kSigmaHf; ps.uhf[i] = ImageF(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[i].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[i].Row(y); for (size_t x = 0; x < xsize; ++x) { - ps.uhf[i].Row(y)[x] = ps.hf[i].Row(y)[x]; + row_uhf[x] = row_hf[x]; } } - ps.hf[i] = DoGBlur(ps.hf[i], kSigmaUhf, 0.0f); - for (size_t y = 0; y < ysize; ++y) { - for (size_t x = 0; x < xsize; ++x) { - ps.uhf[i].Row(y)[x] -= ps.hf[i].Row(y)[x]; + ps.hf[i] = Blur(ps.hf[i], kSigmaUhf, border_hf); + static const double kRemoveHfRange = 0.0287615200377; + static const double kMaxclampHf = 78.8223237675; + static const double kMaxclampUhf = 5.8907152736; + static const float kMulSuppressHf = 1.10684769012; + static const float kMulRegHf = 0.478741530298; + static const float kRegHf = 2000 * kMulRegHf; + static const float kMulSuppressUhf = 1.76905001176; + static const float kMulRegUhf = 0.310148420674; + static const float kRegUhf = 2000 * kMulRegUhf; + + if (i == 0) { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[0].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[0].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_uhf[x] -= row_hf[x]; + row_hf[x] = RemoveRangeAroundZero(kRemoveHfRange, row_hf[x]); + } + } + } else { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_lf = ps.lf[1].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_uhf[x] -= row_hf[x]; + row_hf[x] = MaximumClamp(row_hf[x], kMaxclampHf); + row_uhf[x] = MaximumClamp(row_uhf[x], kMaxclampUhf); + row_uhf[x] = SuppressUhfInBrightAreas(row_uhf[x], row_lf[x], + kMulSuppressUhf, kRegUhf); + row_hf[x] = SuppressHfInBrightAreas(row_hf[x], row_lf[x], + kMulSuppressHf, kRegHf); + + } } } } // Modify range around zero code only concerns the high frequency // planes and only the X and Y channels. - static const double uhf_xy_modification[2] = { - -0.0262070567973, - -5.07470663801, - }; - static const double hf_xy_modification[2] = { - 0.0260892336622, - -0.00789413170469, - }; - static const double mf_xy_modification[2] = { - 0.0185433382632, - -0.158111863182, - }; - ps.uhf = ModifyRangeAroundZero(uhf_xy_modification, ps.uhf); - ps.hf = ModifyRangeAroundZero(hf_xy_modification, ps.hf); - ps.mf = ModifyRangeAroundZero(mf_xy_modification, ps.mf); // Convert low freq xyb to vals space so that we can do a simple squared sum // diff on the low frequencies later. for (size_t y = 0; y < ysize; ++y) { @@ -627,143 +624,46 @@ static void SeparateFrequencies( row_b[x] = valb; } } - // Suppress red-green by intensity change. - static const double suppress[3] = { - -0.0636106621652, - 26.8144000514, - }; - ps.uhf[0] = SuppressXByY(xsize, ysize, ps.uhf[0], ps.uhf[1], suppress[0]); - ps.hf[0] = SuppressXByY(xsize, ysize, ps.hf[0], ps.hf[1], suppress[1]); - static const double maxclamp0 = 0.670004157878; - ps.uhf[0] = MaximumClamping(xsize, ysize, ps.uhf[0], maxclamp0); - static const double maxclamp1 = 2.645076392; - ps.hf[0] = MaximumClamping(xsize, ysize, ps.hf[0], maxclamp1); - static const double maxclamp2 = 64.9667578444; - ps.uhf[1] = MaximumClamping(xsize, ysize, ps.uhf[1], maxclamp2); - static const double maxclamp3 = 79.5957602666; - ps.hf[1] = MaximumClamping(xsize, ysize, ps.hf[1], maxclamp3); - - ps.hf[1] = SuppressHfInBrightAreas(xsize, ysize, ps.hf[1], ps.lf[1]); - ps.uhf[1] = SuppressHfInBrightAreas(xsize, ysize, ps.uhf[1], ps.lf[1]); - ps.mf[1] = SuppressHfInBrightAreas(xsize, ysize, ps.mf[1], ps.lf[1]); } -static void SameNoiseLevelsX(const ImageF& i0, const ImageF& i1, - const double kSigma, - const double w, - const double maxclamp, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { - ImageF blurred0 = CopyPixels(i0); - ImageF blurred1 = CopyPixels(i1); - for (size_t y = 0; y < i0.ysize(); ++y) { - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - for (size_t x = i0.xsize() - 1; x != 0; --x) { - row0[x] -= row0[x - 1]; - row1[x] -= row1[x - 1]; - row0[x] = fabs(row0[x]); - row1[x] = fabs(row1[x]); - if (row0[x] > maxclamp) row0[x] = maxclamp; - if (row1[x] > maxclamp) row1[x] = maxclamp; - } - row0[0] = 0.25 * row0[1]; - row1[0] = 0.25 * row0[1]; - } - blurred0 = Blur(blurred0, kSigma, 0.0); - blurred1 = Blur(blurred1, kSigma, 0.0); +static void SameNoiseLevels(const ImageF& i0, const ImageF& i1, + const double kSigma, + const double w, + const double maxclamp, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + ImageF blurred(i0.xsize(), i0.ysize()); for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); + const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); + float* BUTTERAUGLI_RESTRICT const to = blurred.Row(y); for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = row0[x] - row1[x]; - row_diff[x] += w * diff * diff; + double v0 = fabs(row0[x]); + double v1 = fabs(row1[x]); + if (v0 > maxclamp) v0 = maxclamp; + if (v1 > maxclamp) v1 = maxclamp; + to[x] = v0 - v1; } - } -} -static void SameNoiseLevelsY(const ImageF& i0, const ImageF& i1, - const double kSigma, - const double w, - const double maxclamp, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { - ImageF blurred0 = CopyPixels(i0); - ImageF blurred1 = CopyPixels(i1); - for (size_t y = i0.ysize() - 1; y != 0; --y) { - float* BUTTERAUGLI_RESTRICT const row0prev = blurred0.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row1prev = blurred1.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - for (size_t x = 0; x < i0.xsize(); ++x) { - row0[x] -= row0prev[x]; - row1[x] -= row1prev[x]; - row0[x] = fabs(row0[x]); - row1[x] = fabs(row1[x]); - if (row0[x] > maxclamp) row0[x] = maxclamp; - if (row1[x] > maxclamp) row1[x] = maxclamp; - } } - { - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(0); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(0); - float* BUTTERAUGLI_RESTRICT const row0next = blurred0.Row(1); - float* BUTTERAUGLI_RESTRICT const row1next = blurred1.Row(1); - for (size_t x = 0; x < i0.xsize(); ++x) { - row0[x] = 0.25 * row0next[x]; - row1[x] = 0.25 * row1next[x]; - } - } - blurred0 = Blur(blurred0, kSigma, 0.0); - blurred1 = Blur(blurred1, kSigma, 0.0); + blurred = Blur(blurred, kSigma, 0.0); for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); + const float* BUTTERAUGLI_RESTRICT const row = blurred.Row(y); float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = row0[x] - row1[x]; + double diff = row[x]; row_diff[x] += w * diff * diff; } } } -static void SameNoiseLevelsYP1(const ImageF& i0, const ImageF& i1, - const double kSigma, - const double w, - const double maxclamp, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { - ImageF blurred0 = CopyPixels(i0); - ImageF blurred1 = CopyPixels(i1); - for (size_t y = i0.ysize() - 1; y != 0; --y) { - float* BUTTERAUGLI_RESTRICT const row0prev = blurred0.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row1prev = blurred1.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - for (size_t x = 1; x < i0.xsize(); ++x) { - row0[x] -= row0prev[x - 1]; - row1[x] -= row1prev[x - 1]; - row0[x] = fabs(row0[x]); - row1[x] = fabs(row1[x]); - if (row0[x] > maxclamp) row0[x] = maxclamp; - if (row1[x] > maxclamp) row1[x] = maxclamp; - } - row0[0] = 0.25 * row0[1]; - row1[0] = 0.25 * row1[1]; - } - { - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(0); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(0); - float* BUTTERAUGLI_RESTRICT const row0next = blurred0.Row(1); - float* BUTTERAUGLI_RESTRICT const row1next = blurred1.Row(1); - for (size_t x = 0; x < i0.xsize(); ++x) { - row0[x] = 0.25 * row0next[x]; - row1[x] = 0.25 * row1next[x]; - } +static void L2Diff(const ImageF& i0, const ImageF& i1, const double w, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + if (w == 0) { + return; } - blurred0 = Blur(blurred0, kSigma, 0.0); - blurred1 = Blur(blurred1, kSigma, 0.0); for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); + const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); for (size_t x = 0; x < i0.xsize(); ++x) { double diff = row0[x] - row1[x]; @@ -772,97 +672,47 @@ static void SameNoiseLevelsYP1(const ImageF& i0, const ImageF& i1, } } -static void SameNoiseLevelsYM1(const ImageF& i0, const ImageF& i1, - const double kSigma, - const double w, - const double maxclamp, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { - ImageF blurred0 = CopyPixels(i0); - ImageF blurred1 = CopyPixels(i1); - for (size_t y = i0.ysize() - 1; y != 0; --y) { - float* BUTTERAUGLI_RESTRICT const row0prev = blurred0.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row1prev = blurred1.Row(y - 1); - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - for (size_t x = 0; x + 1 < i0.xsize(); ++x) { - row0[x] -= row0prev[x + 1]; - row1[x] -= row1prev[x + 1]; - row0[x] = fabs(row0[x]); - row1[x] = fabs(row1[x]); - if (row0[x] > maxclamp) row0[x] = maxclamp; - if (row1[x] > maxclamp) row1[x] = maxclamp; - } - row0[i0.xsize() - 1] = 0.25 * row0[i0.xsize() - 2]; - row1[i0.xsize() - 1] = 0.25 * row1[i0.xsize() - 2]; - } - { - float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(0); - float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(0); - float* BUTTERAUGLI_RESTRICT const row0next = blurred0.Row(1); - float* BUTTERAUGLI_RESTRICT const row1next = blurred1.Row(1); - for (size_t x = 0; x < i0.xsize(); ++x) { - row0[x] = 0.25 * row0next[x]; - row1[x] = 0.25 * row1next[x]; - } - } - blurred0 = Blur(blurred0, kSigma, 0.0); - blurred1 = Blur(blurred1, kSigma, 0.0); - for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = blurred0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = blurred1.Row(y); - float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); - for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = row0[x] - row1[x]; - row_diff[x] += w * diff * diff; - } +// i0 is the original image. +// i1 is the deformed copy. +static void L2DiffAsymmetric(const ImageF& i0, const ImageF& i1, + double w_0gt1, + double w_0lt1, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + if (w_0gt1 == 0 && w_0lt1 == 0) { + return; } -} - - -static void L2Diff(const ImageF& i0, const ImageF& i1, const double w, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { + w_0gt1 *= 0.8; + w_0lt1 *= 0.8; for (size_t y = 0; y < i0.ysize(); ++y) { const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); for (size_t x = 0; x < i0.xsize(); ++x) { + // Primary symmetric quadratic objective. double diff = row0[x] - row1[x]; - row_diff[x] += w * diff * diff; - } - } -} - -static void LNDiff(const ImageF& i0, const ImageF& i1, const double w, - double n, - ImageF* BUTTERAUGLI_RESTRICT diffmap) { - if (n == 1.0) { - for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); - float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); - for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = fabs(row0[x] - row1[x]); - row_diff[x] += w * diff; - } - } - } else if (n == 2.0) { - for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); - float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); - for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = row0[x] - row1[x]; - row_diff[x] += w * diff * diff; - } - } - } else { - for (size_t y = 0; y < i0.ysize(); ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); - float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); - for (size_t x = 0; x < i0.xsize(); ++x) { - double diff = fabs(row0[x] - row1[x]); - row_diff[x] += w * pow(diff, n); + row_diff[x] += w_0gt1 * diff * diff; + + // Secondary half-open quadratic objectives. + const double fabs0 = fabs(row0[x]); + const double too_small = 0.4 * fabs0; + const double too_big = 1.0 * fabs0; + + if (row0[x] < 0) { + if (row1[x] > -too_small) { + double v = row1[x] + too_small; + row_diff[x] += w_0lt1 * v * v; + } else if (row1[x] < -too_big) { + double v = -row1[x] - too_big; + row_diff[x] += w_0lt1 * v * v; + } + } else { + if (row1[x] < too_small) { + double v = too_small - row1[x]; + row_diff[x] += w_0lt1 * v * v; + } else if (row1[x] > too_big) { + double v = row1[x] - too_big; + row_diff[x] += w_0lt1 * v * v; + } } } } @@ -887,21 +737,6 @@ ImageF CalculateDiffmap(const ImageF& diffmap_in) { : std::sqrt(orig_val)); } } - { - static const double kSigma = 1.72547472444; - static const double mul1 = 0.458794906198; - static const float scale = 1.0f / (1.0f + mul1); - static const double border_ratio = 1.0; // 2.01209066992; - ImageF blurred = Blur(diffmap, kSigma, border_ratio); - for (int y = 0; y < diffmap.ysize(); ++y) { - const float* const BUTTERAUGLI_RESTRICT row_blurred = blurred.Row(y); - float* const BUTTERAUGLI_RESTRICT row = diffmap.Row(y); - for (int x = 0; x < diffmap.xsize(); ++x) { - row[x] += mul1 * row_blurred[x]; - row[x] *= scale; - } - } - } return diffmap; } @@ -913,9 +748,9 @@ void MaskPsychoImage(const PsychoImage& pi0, const PsychoImage& pi1, std::vector mask_xyb1 = CreatePlanes(xsize, ysize, 3); static const double muls[4] = { 0, - 1.75262681671, - 0.962073813832, - 2.587167299, + 1.64178305129, + 0.831081703362, + 3.23680933546, }; for (int i = 0; i < 2; ++i) { double a = muls[2 * i]; @@ -936,10 +771,12 @@ void MaskPsychoImage(const PsychoImage& pi0, const PsychoImage& pi1, Mask(mask_xyb0, mask_xyb1, mask, mask_dc); } -ButteraugliComparator::ButteraugliComparator(const std::vector& rgb0) +ButteraugliComparator::ButteraugliComparator(const std::vector& rgb0, + double hf_asymmetry) : xsize_(rgb0[0].xsize()), ysize_(rgb0[0].ysize()), - num_pixels_(xsize_ * ysize_) { + num_pixels_(xsize_ * ysize_), + hf_asymmetry_(hf_asymmetry) { if (xsize_ < 8 || ysize_ < 8) return; std::vector xyb0 = OpsinDynamicsImage(rgb0); SeparateFrequencies(xsize_, ysize_, xyb0, pi0_); @@ -982,84 +819,77 @@ void ButteraugliComparator::DiffmapPsychoImage(const PsychoImage& pi1, block_diff_ac[c] = ImageF(xsize_, ysize_, 0.0); } - static const double wUhfMalta = 1.23657307981; - static const double norm1Uhf = 466.149933668; - MaltaDiffMap(pi0_.uhf[1], pi1.uhf[1], wUhfMalta, norm1Uhf, - &block_diff_ac[1]); - - static const double wUhfMaltaX = 3.36199686627; - static const double norm1UhfX = norm1Uhf; - MaltaDiffMap(pi0_.uhf[0], pi1.uhf[0], wUhfMaltaX, norm1UhfX, - &block_diff_ac[0]); - - static const double wHfMalta = 15.6469934822; - static const double norm1Hf = norm1Uhf; - MaltaDiffMap(pi0_.hf[1], pi1.hf[1], wHfMalta, norm1Hf, + static const double wUhfMalta = 5.1409625726; + static const double norm1Uhf = 58.5001247061; + MaltaDiffMap(pi0_.uhf[1], pi1.uhf[1], + wUhfMalta * hf_asymmetry_, + wUhfMalta / hf_asymmetry_, + norm1Uhf, &block_diff_ac[1]); - static const double wHfMaltaX = 129.122071602; - static const double norm1HfX = norm1Uhf; - MaltaDiffMap(pi0_.hf[0], pi1.hf[0], wHfMaltaX, norm1HfX, - &block_diff_ac[0]); - - static const double wMfMaltaX = 51.2720081112; - static const double norm1MfX = norm1Uhf; - MaltaDiffMap(pi0_.mf[0], pi1.mf[0], wMfMaltaX, norm1MfX, + static const double wUhfMaltaX = 4.91743441556; + static const double norm1UhfX = 687196.39002; + MaltaDiffMap(pi0_.uhf[0], pi1.uhf[0], + wUhfMaltaX * hf_asymmetry_, + wUhfMaltaX / hf_asymmetry_, + norm1UhfX, &block_diff_ac[0]); - static const double wmul[11] = { + static const double wHfMalta = 153.671655716; + static const double norm1Hf = 83150785.9592; + MaltaDiffMapLF(pi0_.hf[1], pi1.hf[1], + wHfMalta * sqrt(hf_asymmetry_), + wHfMalta / sqrt(hf_asymmetry_), + norm1Hf, + &block_diff_ac[1]); + + static const double wHfMaltaX = 668.358918152; + static const double norm1HfX = 0.882954368025; + MaltaDiffMapLF(pi0_.hf[0], pi1.hf[0], + wHfMaltaX * sqrt(hf_asymmetry_), + wHfMaltaX / sqrt(hf_asymmetry_), + norm1HfX, + &block_diff_ac[0]); + + static const double wMfMalta = 6841.81248144; + static const double norm1Mf = 0.0135134962487; + MaltaDiffMapLF(pi0_.mf[1], pi1.mf[1], wMfMalta, wMfMalta, norm1Mf, + &block_diff_ac[1]); + + static const double wMfMaltaX = 813.901703816; + static const double norm1MfX = 16792.9322251; + MaltaDiffMapLF(pi0_.mf[0], pi1.mf[0], wMfMaltaX, wMfMaltaX, norm1MfX, + &block_diff_ac[0]); + + static const double wmul[9] = { + 0, + 32.4449876135, + 0, + 0, 0, - 2.52211854569, 0, + 1.01370836411, 0, - 7.34229797917, - 1.92307717196, - 0.779146234988, - 4.91012468367, - 1.83755854086, - 0.0, - 234.519844745, + 1.74566011615, }; - - static const double maxclamp = 72.6815019479; - static const double kSigmaHfX = 10.8163829574; - SameNoiseLevelsX(pi0_.hf[1], pi1.hf[1], kSigmaHfX, wmul[10], maxclamp, - &block_diff_ac[1]); - SameNoiseLevelsY(pi0_.hf[1], pi1.hf[1], kSigmaHfX, wmul[10], maxclamp, - &block_diff_ac[1]); - SameNoiseLevelsYP1(pi0_.hf[1], pi1.hf[1], kSigmaHfX, wmul[10], maxclamp, - &block_diff_ac[1]); - SameNoiseLevelsYM1(pi0_.hf[1], pi1.hf[1], kSigmaHfX, wmul[10], maxclamp, - &block_diff_ac[1]); - - - static const double valn[9] = { - 2.0, - 2.0, - 2.0, - 2.0, - 1.0, - 2.0, - 2.0, - 2.0, - 2.0, - }; + static const double maxclamp = 85.7047444518; + static const double kSigmaHfX = 10.6666499623; + static const double w = 884.809801415; + SameNoiseLevels(pi0_.hf[1], pi1.hf[1], kSigmaHfX, w, maxclamp, + &block_diff_ac[1]); for (int c = 0; c < 3; ++c) { - if (wmul[c] != 0) { - LNDiff(pi0_.hf[c], pi1.hf[c], wmul[c], valn[c], &block_diff_ac[c]); + if (c < 2) { + L2DiffAsymmetric(pi0_.hf[c], pi1.hf[c], + wmul[c] * hf_asymmetry_, + wmul[c] / hf_asymmetry_, + &block_diff_ac[c]); } - LNDiff(pi0_.mf[c], pi1.mf[c], wmul[3 + c], valn[3 + c], &block_diff_ac[c]); - LNDiff(pi0_.lf[c], pi1.lf[c], wmul[6 + c], valn[6 + c], &block_diff_dc[c]); + L2Diff(pi0_.mf[c], pi1.mf[c], wmul[3 + c], &block_diff_ac[c]); + L2Diff(pi0_.lf[c], pi1.lf[c], wmul[6 + c], &block_diff_dc[c]); } - static const double wBlueCorr = 0.0122171286852; - ImageF blurred_b_y_correlation0 = BlurredBlueCorrelation(pi0_.uhf, pi0_.hf); - ImageF blurred_b_y_correlation1 = BlurredBlueCorrelation(pi1.uhf, pi1.hf); - L2Diff(blurred_b_y_correlation0, blurred_b_y_correlation1, wBlueCorr, - &block_diff_ac[2]); - std::vector mask_xyb; std::vector mask_xyb_dc; MaskPsychoImage(pi0_, pi1, xsize_, ysize_, &mask_xyb, &mask_xyb_dc); @@ -1068,10 +898,246 @@ void ButteraugliComparator::DiffmapPsychoImage(const PsychoImage& pi1, CombineChannels(mask_xyb, mask_xyb_dc, block_diff_dc, block_diff_ac)); } -static float MaltaUnit(const float *d, const int xs) { +// Allows PaddedMaltaUnit to call either function via overloading. +struct MaltaTagLF {}; +struct MaltaTag {}; + +static float MaltaUnit(MaltaTagLF, const float* BUTTERAUGLI_RESTRICT d, + const int xs) { + const int xs3 = 3 * xs; + float retval = 0; + { + // x grows, y constant + float sum = + d[-4] + + d[-2] + + d[0] + + d[2] + + d[4]; + retval += sum * sum; + } + { + // y grows, x constant + float sum = + d[-xs3 - xs] + + d[-xs - xs] + + d[0] + + d[xs + xs] + + d[xs3 + xs]; + retval += sum * sum; + } + { + // both grow + float sum = + d[-xs3 - 3] + + d[-xs - xs - 2] + + d[0] + + d[xs + xs + 2] + + d[xs3 + 3]; + retval += sum * sum; + } + { + // y grows, x shrinks + float sum = + d[-xs3 + 3] + + d[-xs - xs + 2] + + d[0] + + d[xs + xs - 2] + + d[xs3 - 3]; + retval += sum * sum; + } + { + // y grows -4 to 4, x shrinks 1 -> -1 + float sum = + d[-xs3 - xs + 1] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 + xs - 1]; + retval += sum * sum; + } + { + // y grows -4 to 4, x grows -1 -> 1 + float sum = + d[-xs3 - xs - 1] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + xs + 1]; + retval += sum * sum; + } + { + // x grows -4 to 4, y grows -1 to 1 + float sum = + d[-4 - xs] + + d[-2 - xs] + + d[0] + + d[2 + xs] + + d[4 + xs]; + retval += sum * sum; + } + { + // x grows -4 to 4, y shrinks 1 to -1 + float sum = + d[-4 + xs] + + d[-2 + xs] + + d[0] + + d[2 - xs] + + d[4 - xs]; + retval += sum * sum; + } + { + /* 0_________ + 1__*______ + 2___*_____ + 3_________ + 4____0____ + 5_________ + 6_____*___ + 7______*__ + 8_________ */ + float sum = + d[-xs3 - 2] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + 2]; + retval += sum * sum; + } + { + /* 0_________ + 1______*__ + 2_____*___ + 3_________ + 4____0____ + 5_________ + 6___*_____ + 7__*______ + 8_________ */ + float sum = + d[-xs3 + 2] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 - 2]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_*_______ + 3__*______ + 4____0____ + 5______*__ + 6_______*_ + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 3] + + d[-xs - 2] + + d[0] + + d[xs + 2] + + d[xs + xs + 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_______*_ + 3______*__ + 4____0____ + 5__*______ + 6_*_______ + 7_________ + 8_________ */ + float sum = + d[-xs - xs + 3] + + d[-xs + 2] + + d[0] + + d[xs - 2] + + d[xs + xs - 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2________* + 3______*__ + 4____0____ + 5__*______ + 6*________ + 7_________ + 8_________ */ + + float sum = + d[xs + xs - 4] + + d[xs - 2] + + d[0] + + d[-xs + 2] + + d[-xs - xs + 4]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2*________ + 3__*______ + 4____0____ + 5______*__ + 6________* + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 4] + + d[-xs - 2] + + d[0] + + d[xs + 2] + + d[xs + xs + 4]; + retval += sum * sum; + } + { + /* 0__*______ + 1_________ + 2___*_____ + 3_________ + 4____0____ + 5_________ + 6_____*___ + 7_________ + 8______*__ */ + float sum = + d[-xs3 - xs - 2] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + xs + 2]; + retval += sum * sum; + } + { + /* 0______*__ + 1_________ + 2_____*___ + 3_________ + 4____0____ + 5_________ + 6___*_____ + 7_________ + 8__*______ */ + float sum = + d[-xs3 - xs + 2] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 + xs - 2]; + retval += sum * sum; + } + return retval; +} + +static float MaltaUnit(MaltaTag, const float* BUTTERAUGLI_RESTRICT d, + const int xs) { const int xs3 = 3 * xs; float retval = 0; - static const float kEdgemul = 0.0309255573587; { // x grows, y constant float sum = @@ -1085,18 +1151,6 @@ static float MaltaUnit(const float *d, const int xs) { d[3] + d[4]; retval += sum * sum; - float sum2 = - d[xs - 4] + - d[xs - 3] + - d[xs - 2] + - d[xs - 1] + - d[xs] + - d[xs + 1] + - d[xs + 2] + - d[xs + 3] + - d[xs + 4]; - float edge = sum - sum2; - retval += kEdgemul * edge * edge; } { // y grows, x constant @@ -1111,18 +1165,6 @@ static float MaltaUnit(const float *d, const int xs) { d[xs3] + d[xs3 + xs]; retval += sum * sum; - float sum2 = - d[-xs3 - xs + 1] + - d[-xs3 + 1] + - d[-xs - xs + 1] + - d[-xs + 1] + - d[1] + - d[xs + 1] + - d[xs + xs + 1] + - d[xs3 + 1] + - d[xs3 + xs + 1]; - float edge = sum - sum2; - retval += kEdgemul * edge * edge; } { // both grow @@ -1157,7 +1199,7 @@ static float MaltaUnit(const float *d, const int xs) { d[-xs] + d[0] + d[xs] + - d[xs - 1] + + d[xs + xs - 1] + d[xs3 - 1] + d[xs3 + xs - 1]; retval += sum * sum; @@ -1171,7 +1213,7 @@ static float MaltaUnit(const float *d, const int xs) { d[-xs] + d[0] + d[xs] + - d[xs + 1] + + d[xs + xs + 1] + d[xs3 + 1] + d[xs3 + xs + 1]; retval += sum * sum; @@ -1372,60 +1414,175 @@ static float MaltaUnit(const float *d, const int xs) { return retval; } -void ButteraugliComparator::MaltaDiffMap( - const ImageF& y0, const ImageF& y1, - const double weight, - const double norm1, - ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const { - PROFILER_FUNC; - const double len = 3.75; - static const double mulli = 0.414888221144; - const double w = mulli * sqrt(weight) / (len * 2 + 1); - const double norm2 = w * norm1; +// Returns MaltaUnit. "fastMode" avoids bounds-checks when x0 and y0 are known +// to be far enough from the image borders. +template +static BUTTERAUGLI_INLINE float PaddedMaltaUnit( + float* const BUTTERAUGLI_RESTRICT diffs, const size_t x0, const size_t y0, + const size_t xsize_, const size_t ysize_) { + int ix0 = y0 * xsize_ + x0; + const float* BUTTERAUGLI_RESTRICT d = &diffs[ix0]; + if (fastMode || + (x0 >= 4 && y0 >= 4 && x0 < (xsize_ - 4) && y0 < (ysize_ - 4))) { + return MaltaUnit(Tag(), d, xsize_); + } + + float borderimage[9 * 9]; + for (int dy = 0; dy < 9; ++dy) { + int y = y0 + dy - 4; + if (y < 0 || y >= ysize_) { + for (int dx = 0; dx < 9; ++dx) { + borderimage[dy * 9 + dx] = 0.0f; + } + } else { + for (int dx = 0; dx < 9; ++dx) { + int x = x0 + dx - 4; + if (x < 0 || x >= xsize_) { + borderimage[dy * 9 + dx] = 0.0f; + } else { + borderimage[dy * 9 + dx] = diffs[y * xsize_ + x]; + } + } + } + } + return MaltaUnit(Tag(), &borderimage[4 * 9 + 4], 9); +} + +template +static void MaltaDiffMapImpl(const ImageF& lum0, const ImageF& lum1, + const size_t xsize_, const size_t ysize_, + const double w_0gt1, + const double w_0lt1, + double norm1, + const double len, const double mulli, + ImageF* block_diff_ac) { + const float kWeight0 = 0.5; + const float kWeight1 = 0.33; + + const double w_pre0gt1 = mulli * sqrt(kWeight0 * w_0gt1) / (len * 2 + 1); + const double w_pre0lt1 = mulli * sqrt(kWeight1 * w_0lt1) / (len * 2 + 1); + const float norm2_0gt1 = w_pre0gt1 * norm1; + const float norm2_0lt1 = w_pre0lt1 * norm1; + std::vector diffs(ysize_ * xsize_); - std::vector sums(ysize_ * xsize_); for (size_t y = 0, ix = 0; y < ysize_; ++y) { - const float* BUTTERAUGLI_RESTRICT const row0 = y0.Row(y); - const float* BUTTERAUGLI_RESTRICT const row1 = y1.Row(y); + const float* BUTTERAUGLI_RESTRICT const row0 = lum0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = lum1.Row(y); for (size_t x = 0; x < xsize_; ++x, ++ix) { - double absval = 0.5 * (std::abs(row0[x]) + std::abs(row1[x])); - double diff = row0[x] - row1[x]; - double scaler = norm2 / (norm1 + absval); + const float absval = 0.5 * std::abs(row0[x]) + 0.5 * std::abs(row1[x]); + const float diff = row0[x] - row1[x]; + const float scaler = norm2_0gt1 / (static_cast(norm1) + absval); + + // Primary symmetric quadratic objective. diffs[ix] = scaler * diff; - } - } - float borderimage[9 * 9]; - for (size_t y0 = 0; y0 < ysize_; ++y0) { - float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); - const bool fastModeY = y0 >= 4 && y0 < ysize_ - 4; - for (size_t x0 = 0; x0 < xsize_; ++x0) { - int ix0 = y0 * xsize_ + x0; - const float *d = &diffs[ix0]; - const bool fastModeX = x0 >= 4 && x0 < xsize_ - 4; - if (fastModeY && fastModeX) { - row_diff[x0] += MaltaUnit(d, xsize_); + + const float scaler2 = norm2_0lt1 / (static_cast(norm1) + absval); + const double fabs0 = fabs(row0[x]); + + // Secondary half-open quadratic objectives. + const double too_small = 0.55 * fabs0; + const double too_big = 1.05 * fabs0; + + if (row0[x] < 0) { + if (row1[x] > -too_small) { + double impact = scaler2 * (row1[x] + too_small); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } else if (row1[x] < -too_big) { + double impact = scaler2 * (-row1[x] - too_big); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } } else { - for (int dy = 0; dy < 9; ++dy) { - int y = y0 + dy - 4; - if (y < 0 || y >= ysize_) { - for (int dx = 0; dx < 9; ++dx) { - borderimage[dy * 9 + dx] = 0; - } + if (row1[x] < too_small) { + double impact = scaler2 * (too_small - row1[x]); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } else if (row1[x] > too_big) { + double impact = scaler2 * (row1[x] - too_big); + if (diff < 0) { + diffs[ix] -= impact; } else { - for (int dx = 0; dx < 9; ++dx) { - int x = x0 + dx - 4; - if (x < 0 || x >= xsize_) { - borderimage[dy * 9 + dx] = 0; - } else { - borderimage[dy * 9 + dx] = diffs[y * xsize_ + x]; - } - } + diffs[ix] += impact; } } - row_diff[x0] += MaltaUnit(&borderimage[4 * 9 + 4], 9); } } } + + size_t y0 = 0; + // Top + for (; y0 < 4; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + for (size_t x0 = 0; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } + + // Middle + for (; y0 < ysize_ - 4; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + size_t x0 = 0; + for (; x0 < 4; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + for (; x0 < xsize_ - 4; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + + for (; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } + + // Bottom + for (; y0 < ysize_; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + for (size_t x0 = 0; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } +} + +void ButteraugliComparator::MaltaDiffMap( + const ImageF& lum0, const ImageF& lum1, + const double w_0gt1, + const double w_0lt1, + const double norm1, ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const { + PROFILER_FUNC; + const double len = 3.75; + static const double mulli = 0.354191303559; + MaltaDiffMapImpl(lum0, lum1, xsize_, ysize_, w_0gt1, w_0lt1, + norm1, len, + mulli, block_diff_ac); +} + +void ButteraugliComparator::MaltaDiffMapLF( + const ImageF& lum0, const ImageF& lum1, + const double w_0gt1, + const double w_0lt1, + const double norm1, ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const { + PROFILER_FUNC; + const double len = 3.75; + static const double mulli = 0.405371989604; + MaltaDiffMapImpl(lum0, lum1, xsize_, ysize_, + w_0gt1, w_0lt1, + norm1, len, + mulli, block_diff_ac); } ImageF ButteraugliComparator::CombineChannels( @@ -1487,48 +1644,44 @@ static std::array MakeMask( } double MaskX(double delta) { - PROFILER_FUNC; - static const double extmul = 2.52662693217; - static const double extoff = 2.0577595478; - static const double offset = 0.342502406734; - static const double scaler = 14.4867545374; - static const double mul = 6.03009840821; + static const double extmul = 2.59885507073; + static const double extoff = 3.08805636789; + static const double offset = 0.315424196682; + static const double scaler = 16.2770141832; + static const double mul = 5.62939030582; static const std::array lut = MakeMask(extmul, extoff, mul, offset, scaler); return InterpolateClampNegative(lut.data(), lut.size(), delta); } double MaskY(double delta) { - PROFILER_FUNC; - static const double extmul = 0.965276993931; - static const double extoff = -0.613819681771; - static const double offset = 1.40903146071; - static const double scaler = 1.07806168416; - static const double mul = 7.09705888614; + static const double extmul = 0.9613705131; + static const double extoff = -0.581933100068; + static const double offset = 1.00846207765; + static const double scaler = 2.2342321176; + static const double mul = 6.64307621174; static const std::array lut = MakeMask(extmul, extoff, mul, offset, scaler); return InterpolateClampNegative(lut.data(), lut.size(), delta); } double MaskDcX(double delta) { - PROFILER_FUNC; - static const double extmul = 10.8596436398; - static const double extoff = 1.58374126704; - static const double offset = 0.651968473749; - static const double scaler = 519.45682322; - static const double mul = 4.72871406401; + static const double extmul = 10.0470705878; + static const double extoff = 3.18472654033; + static const double offset = 0.0551512255218; + static const double scaler = 70.0; + static const double mul = 0.373092999662; static const std::array lut = MakeMask(extmul, extoff, mul, offset, scaler); return InterpolateClampNegative(lut.data(), lut.size(), delta); } double MaskDcY(double delta) { - PROFILER_FUNC; - static const double extmul = 0.00538280872633; - static const double extoff = 59.04237604; - static const double offset = 0.0474092064444; - static const double scaler = 5.52679307489; - static const double mul = 22.7326511523; + static const double extmul = 0.0115640939227; + static const double extoff = 45.9483175519; + static const double offset = 0.0142290066313; + static const double scaler = 5.0; + static const double mul = 2.52611324247; static const std::array lut = MakeMask(extmul, extoff, mul, offset, scaler); return InterpolateClampNegative(lut.data(), lut.size(), delta); @@ -1565,9 +1718,9 @@ ImageF DiffPrecompute(const ImageF& xyb0, const ImageF& xyb1) { fabs(row0_in[x] - row0_in2[x])); double sup1 = (fabs(row1_in[x] - row1_in[x2]) + fabs(row1_in[x] - row1_in2[x])); - static const double mul0 = 0.972407512222; + static const double mul0 = 0.918416534734; row_out[x] = mul0 * std::min(sup0, sup1); - static const double cutoff = 123.915065832; + static const double cutoff = 55.0184555849; if (row_out[x] >= cutoff) { row_out[x] = cutoff; } @@ -1585,42 +1738,57 @@ void Mask(const std::vector& xyb0, const size_t ysize = xyb0[0].ysize(); mask->resize(3); *mask_dc = CreatePlanes(xsize, ysize, 3); - double muls[4] = { - 0.05, - 0.144577484346, - 0.231880902493, - 0.529175844348, + double muls[2] = { + 0.207017089891, + 0.267138152891, }; - double normalizer[2] = { + double normalizer = { 1.0 / (muls[0] + muls[1]), - 1.0 / (muls[2] + muls[3]), }; - static const double r0 = 2.32030744494; - static const double r1 = 7.55507439878; - for (int i = 0; i < 2; ++i) { - (*mask)[i] = ImageF(xsize, ysize); - ImageF diff = DiffPrecompute(xyb0[i], xyb1[i]); - ImageF blurred1 = Blur(diff, r0, 0.0f); - ImageF blurred2 = Blur(diff, r1, 0.0f); + static const double r0 = 2.3770330432; + static const double r1 = 9.04353323561; + static const double r2 = 9.24456601467; + static const double border_ratio = -0.0724948220913; + + { + // X component + ImageF diff = DiffPrecompute(xyb0[0], xyb1[0]); + ImageF blurred = Blur(diff, r2, border_ratio); + (*mask)[0] = ImageF(xsize, ysize); for (size_t y = 0; y < ysize; ++y) { for (size_t x = 0; x < xsize; ++x) { - const double val = normalizer[i] * ( - muls[2 * i + 0] * blurred1.Row(y)[x] + - muls[2 * i + 1] * blurred2.Row(y)[x]); - (*mask)[i].Row(y)[x] = val; + (*mask)[0].Row(y)[x] = blurred.Row(y)[x]; } } } + { + // Y component + (*mask)[1] = ImageF(xsize, ysize); + ImageF diff = DiffPrecompute(xyb0[1], xyb1[1]); + ImageF blurred1 = Blur(diff, r0, border_ratio); + ImageF blurred2 = Blur(diff, r1, border_ratio); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + const double val = normalizer * ( + muls[0] * blurred1.Row(y)[x] + + muls[1] * blurred2.Row(y)[x]); + (*mask)[1].Row(y)[x] = val; + } + } + } + // B component (*mask)[2] = ImageF(xsize, ysize); static const double mul[2] = { - 12.5378252408, - 2.31907764902, + 16.6963293877, + 2.1364621982, }; - static const double w00 = 9.27537465315; - static const double w11 = 2.64039747911; - static const double w_ytob_hf = 1.06493691683; - static const double w_ytob_lf = 9.71657276893; - static const double p1_to_p0 = 0.0153146912176; + static const double w00 = 36.4671237619; + static const double w11 = 2.1887170895; + static const double w_ytob_hf = std::max( + 0.086624184478, + 0.0); + static const double w_ytob_lf = 21.6804277046; + static const double p1_to_p0 = 0.0513061271723; for (size_t y = 0; y < ysize; ++y) { for (size_t x = 0; x < xsize; ++x) { @@ -1641,7 +1809,9 @@ void Mask(const std::vector& xyb0, void ButteraugliDiffmap(const std::vector &rgb0_image, const std::vector &rgb1_image, + double hf_asymmetry, ImageF &result_image) { + PROFILER_FUNC; const size_t xsize = rgb0_image[0].xsize(); const size_t ysize = rgb0_image[0].ysize(); static const int kMax = 8; @@ -1667,7 +1837,7 @@ void ButteraugliDiffmap(const std::vector &rgb0_image, } } ImageF diffmap_scaled; - ButteraugliDiffmap(scaled0, scaled1, diffmap_scaled); + ButteraugliDiffmap(scaled0, scaled1, hf_asymmetry, diffmap_scaled); result_image = ImageF(xsize, ysize); for (int y = 0; y < ysize; ++y) { for (int x = 0; x < xsize; ++x) { @@ -1676,12 +1846,13 @@ void ButteraugliDiffmap(const std::vector &rgb0_image, } return; } - ButteraugliComparator butteraugli(rgb0_image); + ButteraugliComparator butteraugli(rgb0_image, hf_asymmetry); butteraugli.Diffmap(rgb1_image, result_image); } bool ButteraugliInterface(const std::vector &rgb0, const std::vector &rgb1, + float hf_asymmetry, ImageF &diffmap, double &diffvalue) { const size_t xsize = rgb0[0].xsize(); @@ -1695,7 +1866,7 @@ bool ButteraugliInterface(const std::vector &rgb0, return false; // Image planes must have same dimensions. } } - ButteraugliDiffmap(rgb0, rgb1, diffmap); + ButteraugliDiffmap(rgb0, rgb1, hf_asymmetry, diffmap); diffvalue = ButteraugliScoreFromDiffmap(diffmap); return true; } @@ -1724,10 +1895,10 @@ bool ButteraugliAdaptiveQuantization(size_t xsize, size_t ysize, } double ButteraugliFuzzyClass(double score) { - static const double fuzzy_width_up = 6.78721575514; - static const double fuzzy_width_down = 5.96507193294; + static const double fuzzy_width_up = 6.07887388532; + static const double fuzzy_width_down = 5.50793514384; static const double m0 = 2.0; - static const double scaler = 0.861077627013; + static const double scaler = 0.840253347958; double val; if (score < 1.0) { // val in [scaler .. 2.0] diff --git a/butteraugli/butteraugli.h b/butteraugli/butteraugli.h index eef5d89..372a4ba 100755 --- a/butteraugli/butteraugli.h +++ b/butteraugli/butteraugli.h @@ -75,6 +75,7 @@ using ImageF = Image; bool ButteraugliInterface(const std::vector &rgb0, const std::vector &rgb1, + float hf_asymmetry, ImageF &diffmap, double &diffvalue); @@ -423,7 +424,7 @@ struct PsychoImage { class ButteraugliComparator { public: - ButteraugliComparator(const std::vector& rgb0); + ButteraugliComparator(const std::vector& rgb0, double hf_asymmetry); // Computes the butteraugli map between the original image given in the // constructor and the distorted image give here. @@ -440,9 +441,17 @@ class ButteraugliComparator { std::vector* BUTTERAUGLI_RESTRICT mask_dc) const; private: + void MaltaDiffMapLF(const ImageF& y0, + const ImageF& y1, + double w_0gt1, + double w_0lt1, + double normalization, + ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const; + void MaltaDiffMap(const ImageF& y0, const ImageF& y1, - double w, + double w_0gt1, + double w_0lt1, double normalization, ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const; @@ -454,11 +463,13 @@ class ButteraugliComparator { const size_t xsize_; const size_t ysize_; const size_t num_pixels_; + float hf_asymmetry_; PsychoImage pi0_; }; void ButteraugliDiffmap(const std::vector &rgb0, const std::vector &rgb1, + double hf_asymmetry, ImageF &diffmap); double ButteraugliScoreFromDiffmap(const ImageF& distmap); @@ -493,18 +504,18 @@ BUTTERAUGLI_INLINE void OpsinAbsorbance(const V &in0, const V &in1, V *BUTTERAUGLI_RESTRICT out1, V *BUTTERAUGLI_RESTRICT out2) { // https://en.wikipedia.org/wiki/Photopsin absorbance modeling. - static const double mixi0 = 0.254034300884; - static const double mixi1 = 0.466691456721; - static const double mixi2 = 0.0927013001396; - static const double mixi3 = 0.971783366726; - static const double mixi4 = 0.231306455359; - static const double mixi5 = 0.546409605324; - static const double mixi6 = mixi2; - static const double mixi7 = mixi3; - static const double mixi8 = 0.423139962818; - static const double mixi9 = 1.24737070912; - static const double mixi10 = 0.61402451409; - static const double mixi11 = 7.63242153651; + static const double mixi0 = 0.254462330846; + static const double mixi1 = 0.488238255095; + static const double mixi2 = 0.0635278003854; + static const double mixi3 = 1.01681026909; + static const double mixi4 = 0.195214015766; + static const double mixi5 = 0.568019861857; + static const double mixi6 = 0.0860755536007; + static const double mixi7 = 1.1510118369; + static const double mixi8 = 0.07374607900105684; + static const double mixi9 = 0.06142425304154509; + static const double mixi10 = 0.24416850520714256; + static const double mixi11 = 1.20481945273; const V mix0(mixi0); const V mix1(mixi1); diff --git a/butteraugli/butteraugli_main.cc b/butteraugli/butteraugli_main.cc index 80c9510..bd8ae09 100755 --- a/butteraugli/butteraugli_main.cc +++ b/butteraugli/butteraugli_main.cc @@ -371,12 +371,26 @@ int Run(int argc, char* argv[]) { std::vector rgb1 = ReadImageOrDie(argv[1]); std::vector rgb2 = ReadImageOrDie(argv[2]); - if (rgb1[0].xsize() != rgb2[0].xsize() || - rgb1[0].ysize() != rgb2[0].ysize()) { - fprintf( - stderr, "The images are not equal in size: (%lu,%lu) vs (%lu,%lu)\n", - rgb1[0].xsize(), rgb2[0].xsize(), rgb1[0].ysize(), rgb2[0].ysize()); - return 1; + if (rgb1.size() == 3 && rgb2.size() == 4) { + // Adding a missing alpha channel to one of the images. + rgb1.push_back(Image8(rgb1[0].xsize(), rgb1[0].ysize(), 255)); + } else if (rgb2.size() == 3 && rgb1.size() == 4) { + // Adding a missing alpha channel to one of the images. + rgb2.push_back(Image8(rgb2[0].xsize(), rgb2[0].ysize(), 255)); + } else if (rgb1.size() != rgb2.size()) { + fprintf(stderr, "Different number of channels: %lu vs %lu\n", rgb1.size(), + rgb2.size()); + exit(1); + } + + for (size_t c = 0; c < rgb1.size(); ++c) { + if (rgb1[c].xsize() != rgb2[c].xsize() || + rgb1[c].ysize() != rgb2[c].ysize()) { + fprintf( + stderr, "The images are not equal in size: (%lu,%lu) vs (%lu,%lu)\n", + rgb1[c].xsize(), rgb2[c].xsize(), rgb1[c].ysize(), rgb2[c].ysize()); + return 1; + } } // TODO: Figure out if it is a good idea to fetch the gamma from the image @@ -387,8 +401,8 @@ int Run(int argc, char* argv[]) { FromSrgbToLinear(rgb2, linear2, 0); ImageF diff_map, diff_map_on_white; double diff_value; - if (!butteraugli::ButteraugliInterface(linear1, linear2, diff_map, - diff_value)) { + if (!butteraugli::ButteraugliInterface(linear1, linear2, 1.0, + diff_map, diff_value)) { fprintf(stderr, "Butteraugli comparison failed\n"); return 1; } @@ -399,7 +413,8 @@ int Run(int argc, char* argv[]) { FromSrgbToLinear(rgb1, linear1, 255); FromSrgbToLinear(rgb2, linear2, 255); double diff_value_on_white; - if (!butteraugli::ButteraugliInterface(linear1, linear2, diff_map_on_white, + if (!butteraugli::ButteraugliInterface(linear1, linear2, 1.0, + diff_map_on_white, diff_value_on_white)) { fprintf(stderr, "Butteraugli comparison failed\n"); return 1;