diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs index 99f1cf2d88..f65a66df92 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs @@ -6,6 +6,7 @@ using System.Runtime.CompilerServices; using System.Runtime.InteropServices; using SixLabors.ImageSharp.Formats.Jxl.Memory; using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; +using SixLabors.ImageSharp.Formats.Jxl.Processing.Image; using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives; namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; @@ -43,6 +44,26 @@ internal static class Butteraugli private const float GlobalScale = 1.0f / InternalGoodQualityThreshold; +#pragma warning disable // Just indentation warnings + /// + /// Underlying data for . + /// + private static readonly float[,] HeatmapData = + { + {0, 0, 0}, {0, 0, 1}, + {0, 1, 1}, {0, 1, 0}, // Good level + {1, 1, 0}, {1, 0, 0}, // Bad level + {1, 0, 1}, {0.5f, 0.5f, 1.0f}, + {1.0f, 0.5f, 0.5f}, // Pastel colors for the very bad quality range. + {1.0f, 1.0f, 0.5f}, {1, 1, 1}, + {1, 1, 1}, // Last color repeated to have a solid range of white. + }; +#pragma warning restore + + private static readonly DenseMatrix Heatmap; + + static Butteraugli() => Heatmap = new(HeatmapData); + public static ReadOnlySpan Wmul => [ 400.0f, 1.50815703118f, 0f, @@ -72,7 +93,7 @@ internal static class Butteraugli } public static void ConvolveBorderColumn( - JxlImageF input, + JxlPlane input, ReadOnlySpan kernel, int x, Span rowOut) @@ -83,6 +104,7 @@ internal static class Butteraugli int maxX = Math.Min(input.XSize - 1, x + offset); float weight = 0.0f; + for (int j = minX; j <= maxX; j++) { weight += kernel[j - x + offset]; @@ -106,9 +128,9 @@ internal static class Butteraugli } public static bool ConvolutionWithTranspose( - JxlImageF input, + JxlPlane input, ReadOnlySpan kernel, - JxlImageF output) + JxlPlane output) { if (output.XSize != input.YSize) { @@ -286,11 +308,11 @@ internal static class Butteraugli } private static bool Blur( - JxlImageF input, + Configuration configuration, + JxlPlane input, float sigma, - in ButteraugliParameters parameters, ButteraugliBlurTemp temp, - JxlImageF output) + JxlPlane output) { ReadOnlySpan kernel = ComputeKernel(sigma); @@ -323,10 +345,7 @@ internal static class Butteraugli return true; } - if (!temp.GetTransposed(input, out JxlImageF tempT)) - { - return false; - } + JxlPlane tempT = temp.GetTransposed(configuration, input); if (!ConvolutionWithTranspose(input, kernel, tempT)) { @@ -447,9 +466,9 @@ internal static class Butteraugli } } - public static bool SuppressXByY(JxlImageF inY, JxlImageF inOutX) + public static bool SuppressXByY(JxlPlane inY, JxlPlane inOutX) { - if (!SameSize(inOutX, inY)) + if (!JxlImageOperations.SameSize(inOutX, inY)) { return false; } @@ -506,7 +525,7 @@ internal static class Butteraugli } public static bool SeparateLFAndMF( - in ButteraugliParameters parameters, + Configuration configuration, JxlImage3F xyb, JxlImage3F lf, JxlImage3F mf, @@ -517,9 +536,9 @@ internal static class Butteraugli for (int i = 0; i < 3; i++) { if (!Blur( + configuration, xyb.Plane(i), sigmaLf, - parameters, blurTemp, lf.Plane(i))) { @@ -539,18 +558,17 @@ internal static class Butteraugli public static bool SeparateMfAndHf( Configuration configuration, - in ButteraugliParameters parameters, JxlImage3F mf, - ref InlineArray2 hf, - BlurTemp blurTemp) + ref InlineArray2> hf, + ButteraugliBlurTemp blurTemp) { const float sigmaHf = 3.22489901262f; int xSize = mf.XSize; int ySize = mf.YSize; - hf[0] = new JxlImageF(configuration, xSize, ySize); - hf[1] = new JxlImageF(configuration, xSize, ySize); + hf[0] = JxlPlane.Create(configuration, xSize, ySize); + hf[1] = JxlPlane.Create(configuration, xSize, ySize); int lanes = Vector.Count; @@ -558,7 +576,7 @@ internal static class Butteraugli { if (i == 2) { - if (!Blur(mf.Plane(i), sigmaHf, parameters, blurTemp, mf.Plane(i))) + if (!Blur(configuration, mf.Plane(i), sigmaHf, blurTemp, mf.Plane(i))) { return false; } @@ -578,7 +596,7 @@ internal static class Butteraugli } } - if (!Blur(mf.Plane(i), sigmaHf, parameters, blurTemp, mf.Plane(i))) + if (!Blur(configuration, mf.Plane(i), sigmaHf, blurTemp, mf.Plane(i))) { return false; } @@ -631,18 +649,18 @@ internal static class Butteraugli } public static bool SeparateHFAndUHF( - in ButteraugliParameters parameters, - JxlImageF[] hf, - JxlImageF[] uhf, - JxlBlurTemp blurTemp) + Configuration configuration, + InlineArray2> hf, + InlineArray2> uhf, + ButteraugliBlurTemp blurTemp) { const float sigmaUhf = 1.56416327805f; int xSize = hf[0].XSize; int ySize = hf[0].YSize; - uhf[0] = new JxlImageF(xSize, ySize); - uhf[1] = new JxlImageF(xSize, ySize); + uhf[0] = JxlPlane.Create(configuration, xSize, ySize); + uhf[1] = JxlPlane.Create(configuration, xSize, ySize); int lanes = Vector.Count; @@ -659,7 +677,7 @@ internal static class Butteraugli } } - if (!Blur(hf[i], sigmaUhf, parameters, blurTemp, hf[i])) + if (!Blur(configuration, hf[i], sigmaUhf, blurTemp, hf[i])) { return false; } @@ -727,34 +745,36 @@ internal static class Butteraugli return true; } - public static void DeallocateHFAndUHF(InlineArray2 hf, InlineArray2 uhf) + public static void DeallocateHFAndUHF(InlineArray2> hf, InlineArray2> uhf) { for (int i = 0; i < 2; i++) { - hf[i] = new JxlImageF(); - uhf[i] = new JxlImageF(); + hf[i].Dispose(); + uhf[i].Dispose(); + + hf[i] = new JxlPlane(); + uhf[i] = new JxlPlane(); } } public static bool SeparateFrequencies( Configuration configuration, - in ButteraugliParameters parameters, - BlurTemp blurTemp, + ButteraugliBlurTemp blurTemp, JxlImage3F xyb, - PsychoImage ps) + ButteraugliPsychoImage ps) { - ps.Lf = JxlImage3F.Create( + ps.Lf = new JxlImage3F( configuration, xyb.XSize, xyb.YSize); - ps.Mf = JxlImage3F.Create( + ps.Mf = new JxlImage3F( configuration, xyb.XSize, xyb.YSize); if (!SeparateLFAndMF( - parameters, + configuration, xyb, ps.Lf, ps.Mf, @@ -764,16 +784,16 @@ internal static class Butteraugli } if (!SeparateMfAndHf( - parameters, + configuration, ps.Mf, - ps.Hf, + ref ps.Hf, blurTemp)) { return false; } if (!SeparateHFAndUHF( - parameters, + configuration, ps.Hf, ps.Uhf, blurTemp)) @@ -1177,7 +1197,7 @@ internal static class Butteraugli [MethodImpl(MethodImplOptions.AggressiveInlining)] private static float PaddedMaltaUnit( - JxlImageF diffs, + JxlPlane diffs, int x0, int y0, bool isLF) @@ -1238,17 +1258,43 @@ internal static class Butteraugli public static bool MaltaDiffMap( bool isLf, - JxlImageF lum0, - JxlImageF lum1, + JxlPlane lum0, + JxlPlane lum1, + float w0Gt1, + float w0Lt1, + float norm1, + JxlPlane diffs, + JxlPlane blockDiffAc) + { + if (isLf) + { + const float len = 3.75f; + const float mulli = 0.39905817637f; + + return MaltaDiffMap(true, lum0, lum1, w0Gt1, w0Lt1, norm1, len, mulli, diffs, blockDiffAc); + } + else + { + const float len = 3.75f; + const float mulli = 0.611612573796f; + + return MaltaDiffMap(false, lum0, lum1, w0Gt1, w0Lt1, norm1, len, mulli, diffs, blockDiffAc); + } + } + + public static bool MaltaDiffMap( + bool isLf, + JxlPlane lum0, + JxlPlane lum1, float w0Gt1, float w0Lt1, float norm1, float len, float mulli, - JxlImageF diffs, - JxlImageF blockDiffAc) + JxlPlane diffs, + JxlPlane blockDiffAc) { - if (!SameSize(lum0, lum1) || !SameSize(lum0, diffs)) + if (!JxlImageOperations.SameSize(lum0, lum1) || !JxlImageOperations.SameSize(lum0, diffs)) { return false; } @@ -1400,13 +1446,13 @@ internal static class Butteraugli } public static bool MaltaDiffMapLf( - JxlImageF lum0, - JxlImageF lum1, + JxlPlane lum0, + JxlPlane lum1, float w0Gt1, float w0Lt1, float norm1, - JxlImageF diffs, - JxlImageF blockDiffAc) + JxlPlane diffs, + JxlPlane blockDiffAc) { const float len = 3.75f; const float mulli = 0.611612573796f; @@ -1425,9 +1471,9 @@ internal static class Butteraugli } public static void CombineChannelsForMasking( - InlineArray2 hf, - InlineArray2 uhf, - JxlImageF output) + InlineArray2> hf, + InlineArray2> uhf, + JxlPlane output) { // Only X and Y components are involved in masking. ReadOnlySpan muls = @@ -1614,8 +1660,7 @@ internal static class Butteraugli Configuration configuration, JxlImageF mask0, JxlImageF mask1, - in ButteraugliParameters parameters, - JxlBlurTemp blurTemp, + ButteraugliBlurTemp blurTemp, JxlImageF? diffAc, out JxlImageF mask) { @@ -1636,14 +1681,14 @@ internal static class Butteraugli DiffPrecompute(mask0, mul, bias, diff0); DiffPrecompute(mask1, mul, bias, diff1); - if (!Blur(diff0, radius, parameters, blurTemp, blurred0)) + if (!Blur(configuration, diff0, radius, blurTemp, blurred0)) { return false; } FuzzyErosion(blurred0, diff0); - if (!Blur(diff1, radius, parameters, blurTemp, blurred1)) + if (!Blur(configuration, diff1, radius, blurTemp, blurred1)) { return false; } @@ -1661,7 +1706,7 @@ internal static class Butteraugli { maskRow[x] = diff0Row[x]; - if (diffRow != null) + if (diffRow.Length > 0) { const float maskToErrorMul = 10.0f; float diff = blur0Row[x] - blur1Row[x]; @@ -1673,16 +1718,15 @@ internal static class Butteraugli return true; } - public static bool MaskPsychoImage( + public static bool MaskButteraugliPsychoImage( Configuration configuration, ButteraugliPsychoImage pi0, ButteraugliPsychoImage pi1, int width, int height, - in ButteraugliParameters parameters, - BlurTemp blurTemp, + ButteraugliBlurTemp blurTemp, JxlImageF mask, - JxlImageF? diffAc) + out JxlImageF? diffAc) { JxlImageF mask0 = new(configuration, width, height); JxlImageF mask1 = new(configuration, width, height); @@ -1698,12 +1742,12 @@ internal static class Butteraugli mask1); return Mask( + configuration, mask0, mask1, - parameters, blurTemp, mask, - diffAc); + out diffAc); } [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -1741,11 +1785,11 @@ internal static class Butteraugli public static bool CombineChannelsToDiffmap( JxlImageF mask, JxlImage3F blockDiffDc, - JxlImage3F blockDiffAc, + JxlImage3 blockDiffAc, float xmul, JxlImageF result) { - if (!SameSize(mask, result)) + if (!JxlImageOperations.SameSize(mask, result)) { return false; } @@ -1786,10 +1830,10 @@ internal static class Butteraugli } public static void L2Diff( - JxlImageF i0, - JxlImageF i1, + JxlPlane i0, + JxlPlane i1, float w, - JxlImageF diffmap) + JxlPlane diffmap) { if (w == 0) { @@ -1811,10 +1855,10 @@ internal static class Butteraugli } public static void SetL2Diff( - JxlImageF i0, - JxlImageF i1, + JxlPlane i0, + JxlPlane i1, float w, - JxlImageF diffmap) + JxlPlane diffmap) { if (w == 0) { @@ -1836,11 +1880,11 @@ internal static class Butteraugli } public static void L2DiffAsymmetric( - JxlImageF i0, - JxlImageF i1, + JxlPlane i0, + JxlPlane i1, float w0gt1, float w0lt1, - JxlImageF diffmap) + JxlPlane diffmap) { if (w0gt1 == 0 && w0lt1 == 0) { @@ -1942,11 +1986,27 @@ internal static class Butteraugli } } + // A simple HDR-compatible gamma function + private static Vector Gamma(Vector v) + { + Vector kRetMul = Vector.Create(19.245013259874995f * JxlMath.InverseLog2E); + Vector kRetAdd = Vector.Create(-23.16046239805755f); + + // if (value < 0) value = 0; + v = Vector.ConditionalSelect(Vector.LessThan(v, Vector.Zero), Vector.Zero, v); + + Vector biased = v + Vector.Create(9.9710635769299145f); + Vector log = Vector.Log2(biased); + + return (kRetMul * log) + kRetAdd; + } + public static bool OpsinDynamicsImage( + Configuration configuration, JxlImage3F rgb, in ButteraugliParameters parameters, JxlImage3F blurred, - BlurTemp blurTemp, + ButteraugliBlurTemp blurTemp, JxlImage3F xyb) { if (blurred == null) @@ -1954,19 +2014,19 @@ internal static class Butteraugli return false; } - const double sigma = 1.2; + const float sigma = 1.2f; - if (!Blur(rgb.Plane(0), sigma, parameters, blurTemp, blurred.Plane(0))) + if (!Blur(configuration, rgb.Plane(0), sigma, blurTemp, blurred.Plane(0))) { return false; } - if (!Blur(rgb.Plane(1), sigma, parameters, blurTemp, blurred.Plane(1))) + if (!Blur(configuration, rgb.Plane(1), sigma, blurTemp, blurred.Plane(1))) { return false; } - if (!Blur(rgb.Plane(2), sigma, parameters, blurTemp, blurred.Plane(2))) + if (!Blur(configuration, rgb.Plane(2), sigma, blurTemp, blurred.Plane(2))) { return false; } @@ -2058,16 +2118,16 @@ internal static class Butteraugli int xSize = image0.XSize; int ySize = image0.YSize; - using var blurTemp = new JxlBlurTemp(); + using ButteraugliBlurTemp blurTemp = new(); using (JxlImage3F temp = new(configuration, xSize, ySize)) { - if (!OpsinDynamicsImage(image0, parameters, temp, blurTemp, image0)) + if (!OpsinDynamicsImage(configuration, image0, parameters, temp, blurTemp, image0)) { return false; } - if (!OpsinDynamicsImage(image1, parameters, temp, blurTemp, image1)) + if (!OpsinDynamicsImage(configuration, image1, parameters, temp, blurTemp, image1)) { return false; } @@ -2080,12 +2140,12 @@ internal static class Butteraugli using (JxlImage3F lf0 = new(configuration, xSize, ySize)) using (JxlImage3F lf1 = new(configuration, xSize, ySize)) { - if (!SeparateLFAndMF(parameters, image0, lf0, image0, blurTemp)) + if (!SeparateLFAndMF(configuration, image0, lf0, image0, blurTemp)) { return false; } - if (!SeparateLFAndMF(parameters, image1, lf1, image1, blurTemp)) + if (!SeparateLFAndMF(configuration, image1, lf1, image1, blurTemp)) { return false; } @@ -2100,15 +2160,15 @@ internal static class Butteraugli } } - InlineArray2 hf0 = default; - InlineArray2 hf1 = default; + InlineArray2> hf0 = default; + InlineArray2> hf1 = default; - if (!SeparateMfAndHf(parameters, image0, hf0, blurTemp)) + if (!SeparateMfAndHf(configuration, image0, ref hf0, blurTemp)) { return false; } - if (!SeparateMfAndHf(parameters, image1, hf1, blurTemp)) + if (!SeparateMfAndHf(configuration, image1, ref hf1, blurTemp)) { return false; } @@ -2158,15 +2218,15 @@ internal static class Butteraugli image0.Dispose(); image1.Dispose(); - InlineArray2 uhf0 = default; - InlineArray2 uhf1 = default; + InlineArray2> uhf0 = default; + InlineArray2> uhf1 = default; - if (!SeparateHFAndUHF(parameters, hf0, uhf0, blurTemp)) + if (!SeparateHFAndUHF(configuration, hf0, uhf0, blurTemp)) { return false; } - if (!SeparateHFAndUHF(parameters, hf1, uhf1, blurTemp)) + if (!SeparateHFAndUHF(configuration, hf1, uhf1, blurTemp)) { return false; } @@ -2175,7 +2235,7 @@ internal static class Butteraugli using (JxlImageF diffs = new(configuration, xSize, ySize)) { - MaltaDiffMap( + _ = MaltaDiffMap( false, uhf0[1], uhf1[1], @@ -2185,19 +2245,19 @@ internal static class Butteraugli diffs, blockDiffAc); - MaltaDiffMap( + _ = MaltaDiffMap( false, uhf0[0], uhf1[0], - wUhfMaltaX * hfAsymmetry, - wUhfMaltaX / hfAsymmetry, - norm1UhfX, + WUhfMaltaX * hfAsymmetry, + WUhfMaltaX / hfAsymmetry, + Norm1UhfX, diffs, blockDiffAc); float sqrtAsym = MathF.Sqrt(hfAsymmetry); - MaltaDiffMap( + _ = MaltaDiffMap( true, hf0[1], hf1[1], @@ -2207,7 +2267,7 @@ internal static class Butteraugli diffs, blockDiffAc); - MaltaDiffMap( + _ = MaltaDiffMap( true, hf0[0], hf1[0], @@ -2239,15 +2299,15 @@ internal static class Butteraugli DeallocateHFAndUHF(hf0, uhf0); DeallocateHFAndUHF(hf1, uhf1); - if (!Mask(mask0, mask1, parameters, blurTemp, mask, blockDiffAc)) + if (!Mask(configuration, mask0, mask1, blurTemp, mask, out JxlImageF newBlockDiffAc)) { return false; } for (int y = 0; y < ySize; y++) { - ReadOnlySpan dc = blockDiffDc.GetRow(y); - ReadOnlySpan ac = blockDiffAc.GetRow(y); + ReadOnlySpan dc = newBlockDiffAc.GetRow(y); + ReadOnlySpan ac = newBlockDiffAc.GetRow(y); Span output = diffmap.GetRow(y); ReadOnlySpan maskRow = mask.GetRow(y); @@ -2255,10 +2315,7 @@ internal static class Butteraugli { float m = maskRow[x]; - output[x] = - MathF.Sqrt( - (dc[x] * (float)MaskDcY(m)) + - (ac[x] * (float)MaskY(m))); + output[x] = MathF.Sqrt((dc[x] * (float)MaskDcY(m)) + (ac[x] * (float)MaskY(m))); } } @@ -2272,7 +2329,7 @@ internal static class Butteraugli int xs = (input.XSize + 1) / 2; int ys = (input.YSize + 1) / 2; - JxlImage3F retval = new(Configuration, xs, ys); + JxlImage3F retval = new(configuration, xs, ys); for (int c = 0; c < 3; ++c) { @@ -2293,8 +2350,7 @@ internal static class Butteraugli for (int x = 0; x < input.XSize; ++x) { - retval.PlaneRow(c, y / 2)[x / 2] += - 0.25f * srcRow[x]; + retval.PlaneRow(c, y / 2)[x / 2] += 0.25f * srcRow[x]; } } @@ -2335,4 +2391,286 @@ internal static class Butteraugli } } } + + private static void ScoreToRgb(float score, float goodThreshold, float badThreshold, ref InlineArray3 rgb) + { + if (score < goodThreshold) + { + score = (score / goodThreshold) * 0.3f; + } + else if (score < badThreshold) + { + score = 0.3f + ((score - goodThreshold) / (badThreshold - goodThreshold) * 0.15f); + } + else + { + score = 0.45f + ((score - badThreshold) / (badThreshold * 12) * 0.5f); + } + + int tableSize = HeatmapData.GetLength(0); + score = Math.Clamp(score * (tableSize - 1), 0f, tableSize - 2); + + int ix = (int)score; + ix = Math.Clamp(ix, 0, tableSize - 2); // handle NaN + float mix = score - ix; + + for (int i = 0; i < 3; ++i) + { + float v = (mix * Heatmap[ix + 1, i]) + ((1 - mix) * Heatmap[ix, i]); + rgb[i] = MathF.Pow(v, 0.5f); + } + } + + public static JxlImage3F CreateHeatMapImage(Configuration configuration, JxlImageF distmap, float goodThreshold, float badThreshold) + { + // Do not dispose (this is the return value; it is caller's responsibility to dispose this + // or keep it) + JxlImage3F heatmap = new(configuration, distmap.XSize, distmap.YSize); + + for (int y = 0; y < distmap.YSize; y++) + { + Span row_distmap = distmap.GetRow(y); + Span row_h0 = heatmap.PlaneRow(0, y); + Span row_h1 = heatmap.PlaneRow(1, y); + Span row_h2 = heatmap.PlaneRow(2, y); + + for (int x = 0; x < distmap.XSize; ++x) + { + float d = row_distmap[x]; + InlineArray3 rgb = default; + ScoreToRgb(d, goodThreshold, badThreshold, ref rgb); + row_h0[x] = rgb[0]; + row_h1[x] = rgb[1]; + row_h2[x] = rgb[2]; + } + } + + return heatmap; + } + + public static float ButteraugliFuzzyInverse(float seek) + { + float pos = 0f; + + for (float range = 1.0f; range >= 1e-10f; range *= 0.5f) + { + float cur = ButteraugliFuzzyClass(pos); + + if (cur < seek) + { + pos -= range; + } + else + { + pos += range; + } + } + + // Normalization (pos) can be printed if seek == 1.0, for example + // when debugging. + return pos; + } + + public static float ButteraugliFuzzyClass(float score) + { + const float fuzzyWidthUp = 4.8f; + const float fuzzyWidthDown = 4.8f; + const float m0 = 2.0f; + const float scaler = 0.7777f; + + float val; + + if (score < 1.0) + { + val = m0 / (1.0f + MathF.Exp((score - 1.0f) * fuzzyWidthDown)); + val -= 1.0f; + val *= 2.0f - scaler; + val += scaler; + } + else + { + val = m0 / (1.0f + MathF.Exp((score - 1.0f) * fuzzyWidthUp)); + val *= scaler; + } + + return val; + } + + public static bool ButteraugliInterfaceInPlace(Configuration configuration, JxlImage3F rgb0, JxlImage3F rgb1, ButteraugliParameters parameters, JxlImageF diffmap, out double diffvalue) + { + diffvalue = 0; + + int xsize = rgb0.XSize; + int ysize = rgb0.YSize; + + if ((xsize & ysize) == 0) // equivalent to xsize == 0 || ysize == 0, minus one branch + { + throw new InvalidOperationException("Zero-sized image"); + } + + if (!JxlImageOperations.SameSize(rgb0, rgb1)) + { + throw new InvalidOperationException("Size mismatch"); + } + + const int max = 8; + + if (xsize < max || ysize < max) + { + bool ok = ButteraugliDiffmapSmall(configuration, max, rgb0, rgb1, parameters, out diffmap); + diffvalue = ButteraugliScoreFromDiffmap(diffmap, parameters); + return ok; + } + + JxlImageF subDiffmap = new(); + + if (xsize >= 15 && ysize >= 15) + { + using JxlImage3F rgb0Sub = SubSample2x(configuration, rgb0); + using JxlImage3F rgb1Sub = SubSample2x(configuration, rgb1); + + if (!ButteraugliDiffmapInPlace(configuration, rgb0Sub, rgb1Sub, parameters, diffmap)) + { + return false; + } + } + + if (!ButteraugliDiffmapInPlace(configuration, rgb0, rgb1, parameters, diffmap)) + { + return false; + } + + if (xsize >= 15 && ysize >= 15) + { + AddSupersampled2x(subDiffmap, 0.5f, diffmap); + } + + diffvalue = ButteraugliScoreFromDiffmap(diffmap, parameters); + return true; + } + + public static bool ButteraugliInterface(Configuration configuration, JxlImage3F rgb0, JxlImage3F rgb1, ButteraugliParameters parameters, JxlImageF diffmap, out double diffValue) + { + diffValue = 0; + + if (!ButteraugliDiffmap(configuration, rgb0, rgb1, parameters, out diffmap)) + { + return false; + } + + diffValue = ButteraugliScoreFromDiffmap(diffmap, parameters); + return true; + } + + public static bool ButteraugliInterface(Configuration configuration, JxlImage3F rgb0, JxlImage3F rgb1, float hfAsymmetry, float xmul, JxlImageF diffmap, out double diffValue) + { + ButteraugliParameters parameters = new() + { + HfAsymmetry = hfAsymmetry, + XMultiplier = xmul + }; + + return ButteraugliInterface(configuration, rgb0, rgb1, parameters, diffmap, out diffValue); + } + + public static bool ButteraugliDiffmap(Configuration configuration, JxlImage3F rgb0, JxlImage3F rgb1, ButteraugliParameters parameters, out JxlImageF diffmap) + { + diffmap = new(); + + int xsize = rgb0.XSize; + int ysize = rgb0.YSize; + + if ((xsize & ysize) == 0) // equivalent to xsize == 0 || ysize == 0, minus one branch + { + throw new InvalidOperationException("Zero-sized image"); + } + + if (!JxlImageOperations.SameSize(rgb0, rgb1)) + { + throw new InvalidOperationException("Size mismatch"); + } + + const int max = 8; + + if (xsize < max || ysize < max) + { + if (!ButteraugliDiffmapSmall(configuration, max, rgb0, rgb1, parameters, out diffmap)) + { + return false; + } + } + + ButteraugliComparator butteraugli = ButteraugliComparator.Make(configuration, rgb0, parameters); + + return butteraugli.Diffmap(configuration, rgb1, diffmap); + } + + public static bool ButteraugliDiffmapSmall(Configuration configuration, int max, JxlImage3F rgb0, JxlImage3F rgb1, ButteraugliParameters parameters, out JxlImageF diffmap) + { + int xsize = rgb0.XSize; + int ysize = rgb0.YSize; + + int xborder = xsize < max ? (max - xsize) / 2 : 0; + int yborder = ysize < max ? (max - ysize) / 2 : 0; + int xscaled = Math.Max(max, xsize); + int yscaled = Math.Max(max, ysize); + + using JxlImage3F scaled0 = new(configuration, xscaled, yscaled); + using JxlImage3F scaled1 = new(configuration, xscaled, yscaled); + + for (int i = 0; i < 3; ++i) + { + for (int y = 0; y < yscaled; y++) + { + for (int x = 0; x < xscaled; x++) + { + int x2 = Math.Min(xsize - 1, x > xborder ? x - xborder : 0); + int y2 = Math.Min(ysize - 1, y > yborder ? y - yborder : 0); + + scaled0.PlaneRow(i, y)[x] = rgb0.PlaneRow(i, y2)[x2]; + scaled1.PlaneRow(i, y)[x] = rgb1.PlaneRow(i, y2)[x2]; + } + } + } + + JxlImageF diffmapScaled = new(); + bool ok = ButteraugliDiffmap(configuration, scaled0, scaled1, parameters, out diffmapScaled); + + diffmap = new(configuration, xsize, ysize); + + for (int y = 0; y < ysize; y++) + { + Span diffmapRow = diffmap.GetRow(y); + Span diffmapScaledRow = diffmapScaled.GetRow(y + yborder); + + for (int x = 0; x < xsize; x++) + { + diffmapRow[x] = diffmapScaledRow[x + xborder]; + } + } + + return ok; + } + + public static float ButteraugliScoreFromDiffmap(JxlImageF diffmap, ButteraugliParameters parameters) + { + float retval = 0.0f; + + for (int y = 0; y < diffmap.YSize; ++y) + { + Span row = diffmap.GetRow(y); + for (int x = 0; x < diffmap.XSize; ++x) + { + retval = Math.Max(retval, row[x]); + } + } + + return retval; + } + + public static bool MaltaDiffMap(JxlPlane lum0, JxlPlane lum1, float w0Gt1, float w0Lt1, float norm1, JxlPlane diffs, JxlImage3 blockDiffAc, int c) + => MaltaDiffMap(isLf: false, lum0, lum1, w0Gt1, w0Lt1, norm1, diffs, blockDiffAc.Plane(c)); + + public static bool MaltaDiffMapLf(JxlPlane lum0, JxlPlane lum1, float w0Gt1, float w0Lt1, float norm1, JxlPlane diffs, JxlImage3 blockDiffAc, int c) + => MaltaDiffMap(isLf: true, lum0, lum1, w0Gt1, w0Lt1, norm1, diffs, blockDiffAc.Plane(c)); } diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliBlurTemp.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliBlurTemp.cs new file mode 100644 index 0000000000..26aa8372ce --- /dev/null +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliBlurTemp.cs @@ -0,0 +1,24 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using SixLabors.ImageSharp.Formats.Jxl.Memory; + +namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; + +internal sealed class ButteraugliBlurTemp : IDisposable +{ + public JxlPlane TransposedTemp { get; set; } = new(); + + public JxlPlane GetTransposed(Configuration configuration, JxlPlane input) + { + if (this.TransposedTemp.XSize == 0) + { + // Yes, YSize and XSize are swapped + this.TransposedTemp = JxlPlane.Create(configuration, input.YSize, input.XSize); + } + + return this.TransposedTemp; + } + + public void Dispose() => this.TransposedTemp.Dispose(); +} diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs index be5837534b..d290268da6 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs @@ -1,46 +1,300 @@ // Copyright (c) Six Labors. // Licensed under the Six Labors Split License. +using SixLabors.ImageSharp.Formats.Jxl.Memory; using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; +using SixLabors.ImageSharp.Formats.Jxl.Processing.Image; namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; -internal sealed class ButteraugliComparator +internal class ButteraugliComparator : IDisposable { - private int xSize; - private int ySize; + private readonly int xSize; + private readonly int ySize; private ButteraugliParameters parameters; - private readonly ButteraugliPsychoImage pi0; - private readonly JxlImage3F temp; - private bool tempInUse; - private readonly ButteraugliBlurTemp blurTemp; + private readonly ButteraugliPsychoImage pi0 = new(); + private JxlImage3F? temp; + private readonly ButteraugliBlurTemp blurTemp = new(); + private ButteraugliComparator? sub; + + public ButteraugliComparator(int xsize, int ysize, ButteraugliParameters parameters) + { + this.xSize = xsize; + this.ySize = ysize; + this.parameters = parameters; + } + + public JxlImage3F Temp => this.temp ??= new(); + + public static ButteraugliComparator Make(Configuration configuration, JxlImage3F rgb0, ButteraugliParameters parameters) + { + int xSize = rgb0.XSize; + int ySize = rgb0.YSize; + + ButteraugliComparator result = new(xSize, ySize, parameters) + { + temp = new JxlImage3F(configuration, xSize, ySize) + }; + + if (xSize < 8 || ySize < 8) + { + return result; + } + + JxlImage3F xyb0 = new(configuration, xSize, ySize); + + if (!Butteraugli.OpsinDynamicsImage(configuration, rgb0, parameters, result.Temp, result.blurTemp, xyb0)) + { + throw new InvalidOperationException("OpsinDynamicsImage failed"); + } + + result.ReleaseTemp(); + + if (!Butteraugli.SeparateFrequencies(configuration, result.blurTemp, xyb0, result.pi0)) + { + throw new InvalidOperationException("Could not separate frequencies"); + } + + JxlImage3F subsampledRgb0 = Butteraugli.SubSample2x(configuration, rgb0); + result.sub = Make(configuration, subsampledRgb0, parameters); + + return result; + } + + public void ReleaseTemp() + { + this.temp?.Dispose(); + this.temp = null; + } /// /// Computes the butteraugli map between the original image given in the constructor and the distorted image given here. /// - public bool Diffmap(JxlImage3F rgb1, JxlImageF result) + public virtual bool Diffmap(Configuration configuration, JxlImage3F rgb1, JxlImageF result) { - throw new NotImplementedException(); + if (this.xSize < 8 || this.ySize < 8) + { + result.Clear(); + return true; + } + + JxlImage3F xyb1 = new(configuration, this.xSize, this.ySize); + + if (!Butteraugli.OpsinDynamicsImage(configuration, rgb1, this.parameters, this.Temp, this.blurTemp, xyb1)) + { + return false; + } + + this.ReleaseTemp(); + + if (!this.DiffmapOpsinDynamicsImage(configuration, xyb1, out result)) + { + return false; + } + + if (this.sub is not null) + { + if (this.sub.xSize < 8 || this.sub.ySize < 8) + { + return true; + } + + JxlImage3F subXyb = new(configuration, this.sub.xSize, this.sub.ySize); + JxlImage3F subsampledRgb1 = Butteraugli.SubSample2x(configuration, rgb1); + + if (!Butteraugli.OpsinDynamicsImage(configuration, subsampledRgb1, this.parameters, this.sub.Temp, this.sub.blurTemp, subXyb)) + { + return false; + } + + this.sub.ReleaseTemp(); + + if (!this.DiffmapOpsinDynamicsImage(configuration, subXyb, out JxlImageF subResult)) + { + return false; + } + + Butteraugli.AddSupersampled2x(subResult, 0.5f, result); + } + + return true; } /// /// Same as Diffmap but OpsinDynamicsImage() was already applied. /// - public bool DiffmapOpsinDynamicsImage(JxlImage3F xyb1, JxlImageF result) + public bool DiffmapOpsinDynamicsImage(Configuration configuration, JxlImage3F xyb1, out JxlImageF result) { - throw new NotImplementedException(); + result = new(); + + if (this.xSize < 8 || this.ySize < 8) + { + result.Clear(); + return true; + } + + ButteraugliPsychoImage pi1 = new(); + + if (!Butteraugli.SeparateFrequencies(configuration, this.blurTemp, xyb1, pi1)) + { + return false; + } + + result = new(configuration, this.xSize, this.ySize); + return this.DiffmapPsychoImage(configuration, pi1, result); } /// /// Same as above but the frequency decomposition was already applied. /// - public bool DiffmapPsychoImage(ButteraugliPsychoImage pi1, JxlImageF diffmap) + public bool DiffmapPsychoImage(Configuration configuration, ButteraugliPsychoImage pi1, JxlImageF diffmap) { - throw new NotImplementedException(); + if (this.xSize < 8 || this.ySize < 8) + { + diffmap.Clear(); + return true; + } + + float hfAsymmetry = this.parameters.HfAsymmetry; + float xmul = this.parameters.XMultiplier; + + JxlImageF diffs = new(configuration, this.xSize, this.ySize); + JxlImage3 blockDiffAc = new(configuration, this.xSize, this.ySize); + + JxlImageOperations.ZeroFillImage(blockDiffAc); + + if (!Butteraugli.MaltaDiffMap( + this.pi0.Uhf[1], + pi1.Uhf[1], + Butteraugli.WUhfMalta * hfAsymmetry, + Butteraugli.WUhfMalta / hfAsymmetry, + Butteraugli.Norm1Uhf, + diffs, + blockDiffAc, + 1)) + { + return false; + } + + if (!Butteraugli.MaltaDiffMap( + this.pi0.Uhf[0], + pi1.Uhf[0], + Butteraugli.WUhfMaltaX * hfAsymmetry, + Butteraugli.WUhfMaltaX / hfAsymmetry, + Butteraugli.Norm1UhfX, + diffs, + blockDiffAc, + 0)) + { + return false; + } + + if (!Butteraugli.MaltaDiffMapLf( + this.pi0.Hf[1], + pi1.Hf[1], + Butteraugli.WHfMalta * MathF.Sqrt(hfAsymmetry), + Butteraugli.WHfMalta / MathF.Sqrt(hfAsymmetry), + Butteraugli.Norm1Hf, + diffs, + blockDiffAc, + 1)) + { + return false; + } + + if (!Butteraugli.MaltaDiffMapLf( + this.pi0.Hf[0], + pi1.Hf[0], + Butteraugli.WHfMaltaX * MathF.Sqrt(hfAsymmetry), + Butteraugli.WHfMaltaX / MathF.Sqrt(hfAsymmetry), + Butteraugli.Norm1HfX, + diffs, + blockDiffAc, + 0)) + { + return false; + } + + if (!Butteraugli.MaltaDiffMapLf( + this.pi0.Mf!.Plane(1), + pi1.Mf!.Plane(1), + Butteraugli.WHfMalta, + Butteraugli.WMfMalta, + Butteraugli.Norm1Mf, + diffs, + blockDiffAc, + 1)) + { + return false; + } + + if (!Butteraugli.MaltaDiffMapLf( + this.pi0.Mf!.Plane(0), + pi1.Mf!.Plane(0), + Butteraugli.WHfMaltaX, + Butteraugli.WMfMaltaX, + Butteraugli.Norm1MfX, + diffs, + blockDiffAc, + 0)) + { + return false; + } + + JxlImage3F blockDiffDc = new(configuration, this.xSize, this.ySize); + + for (int c = 0; c < 3; c++) + { + if (c < 2) + { + Butteraugli.L2DiffAsymmetric( + this.pi0.Hf[c], + pi1.Hf[c], + Butteraugli.Wmul[c] * hfAsymmetry, + Butteraugli.Wmul[c] / hfAsymmetry, + blockDiffAc.Plane(c)); + } + + Butteraugli.L2Diff( + this.pi0.Mf.Plane(c), + pi1.Mf.Plane(c), + Butteraugli.Wmul[3 + c], + blockDiffAc.Plane(c)); + + Butteraugli.SetL2Diff( + this.pi0.Lf!.Plane(c), + pi1.Lf!.Plane(c), + Butteraugli.Wmul[6 + c], + blockDiffDc.Plane(c)); + } + + JxlImageF mask = new(); + + if (!Butteraugli.MaskButteraugliPsychoImage( + configuration, + this.pi0, + pi1, + this.xSize, + this.ySize, + this.blurTemp, + mask, + out JxlImageF? diffAc)) + { + return false; + } + + diffAc?.BytesSpan.CopyTo(blockDiffAc.Plane(1).BytesSpan); + + return Butteraugli.CombineChannelsToDiffmap(mask, blockDiffDc, blockDiffAc, xmul, diffmap); } - public bool Mask(JxlImageF mask) + public virtual bool Mask(Configuration configuration, JxlImageF mask) + => Butteraugli.MaskButteraugliPsychoImage(configuration, this.pi0, this.pi0, this.xSize, this.ySize, this.blurTemp, mask, out _); + + public void Dispose() { - throw new NotImplementedException(); + this.temp?.Dispose(); + this.temp = null; + this.blurTemp.Dispose(); } } diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliPsychoImage.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliPsychoImage.cs new file mode 100644 index 0000000000..de1166cf8e --- /dev/null +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliPsychoImage.cs @@ -0,0 +1,23 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System.Runtime.CompilerServices; +using SixLabors.ImageSharp.Formats.Jxl.Memory; +using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; + +namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; + +// Public fields because InlineArray2 values +// cannot be mutated if it's a property +#pragma warning disable SA1401 // Fields should be private + +internal sealed class ButteraugliPsychoImage +{ + public InlineArray2> Uhf; // XY + + public InlineArray2> Hf; // XY + + public JxlImage3F? Mf { get; set; } // XYB + + public JxlImage3F? Lf { get; set; } // XYB +} diff --git a/src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageOperations.cs b/src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageOperations.cs index 0a4ff8f5d7..b4e98c5ca4 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageOperations.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/Image/JxlImageOperations.cs @@ -23,6 +23,15 @@ internal static class JxlImageOperations /// True if width and height is equal. public static bool SameSize(JxlPlaneBase a, JxlPlaneBase b) => a.XSize == b.XSize && a.YSize == b.YSize; + /// + /// Returns true if first image has same width and height as the second image. + /// + /// First image + /// Second image + /// True if width and height is equal. + public static bool SameSize(JxlImage3 a, JxlImage3 b) + where T : unmanaged => a.XSize == b.XSize && a.YSize == b.YSize; + /// /// Copies everything from one plane to another. /// diff --git a/src/ImageSharp/Formats/Jxl/Processing/JxlToc.cs b/src/ImageSharp/Formats/Jxl/Processing/JxlToc.cs index 21fa5d23b2..b6dfd4e39f 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/JxlToc.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/JxlToc.cs @@ -31,6 +31,7 @@ internal static class JxlToc } private const int BitsPerByte = 8; + private const int MaxTocEntries = 65536; public static bool ReadToc( @@ -98,10 +99,7 @@ internal static class JxlToc } } - if (!reader.JumpToByteBoundary()) - { - return false; - } + reader.JumpToByteBoundary(); if (!CheckBitBudget(tocEntries)) { @@ -113,10 +111,7 @@ internal static class JxlToc sizes[i] = JxlU32Coder.Read(TocDistribution, reader); } - if (!reader.JumpToByteBoundary()) - { - return false; - } + reader.JumpToByteBoundary(); return CheckBitBudget(0); }