diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs new file mode 100644 index 000000000..eb25583ec --- /dev/null +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/Butteraugli.cs @@ -0,0 +1,2333 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using System.Numerics; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; +using SixLabors.ImageSharp.Formats.Jxl.Memory; +using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; + +namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; + +/// +/// Implementation of Google Butteraugli, an advanced +/// image comparison system. Unlike PSNR or SSIM, +/// Butteraugli compares images the way humans may +/// spot differences instead of pixelwise. +/// +internal static class Butteraugli +{ + // NOTE: The meaning of those constants + // is not perfectly understood. + public const float WMfMalta = 37.0819870399f; + public const float Norm1Mf = 130262059.556f; + public const float WMfMaltaX = 8246.75321353f; + public const float Norm1MfX = 1009002.70582f; + + public const float WHfMalta = 18.7237414387f; + public const float Norm1Hf = 4498534.45232f; + public const float WHfMaltaX = 6923.99476109f; + public const float Norm1HfX = 8051.15833247f; + + public const float WUhfMalta = 1.10039032555f; + public const float Norm1Uhf = 71.7800275169f; + public const float WUhfMaltaX = 173.5f; + public const float Norm1UhfX = 5.0f; + + private const float IntensityTargetNormalizationHack = 0.79079917404f; + + private static readonly float InternalGoodQualityThreshold = + 17.83f * IntensityTargetNormalizationHack; + + private static readonly float GlobalScale = + 1.0f / InternalGoodQualityThreshold; + + public static ReadOnlySpan Wmul => + [ + 400.0, 1.50815703118, 0, + 2150.0, 10.6195433239, 16.2176043152, + 29.2353797994, 0.844626970982, 0.703646627719, + ]; + + public static ReadOnlySpan ComputeKernel(float sigma) + { + const float m = 2.25f; // Accuracy increases when m is increased + float scaler = -1.0f / (2.0f * sigma * sigma); + int diff = Math.Max(1, (int)(m * MathF.Abs(sigma))); + + // Use new because there's only up to 3 elements + float[] kernel = new float[(2 * diff) + 1]; + + for (int i = -diff; i <= diff; i++) + { + kernel[i + diff] = MathF.Exp(scaler * i * i); + } + + return kernel; + } + + public static void ConvolveBorderColumn( + JxlImageF input, + ReadOnlySpan kernel, + int x, + Span rowOut) + { + int offset = kernel.Length / 2; + + int minX = x < offset ? 0 : x - offset; + 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]; + } + + float scale = 1.0f / weight; + + for (int y = 0; y < input.YSize; y++) + { + Span rowIn = input.GetRow(y); + + float sum = 0.0f; + + for (int j = minX; j <= maxX; j++) + { + sum += rowIn[j] * kernel[j - x + offset]; + } + + rowOut[y] = sum * scale; + } + } + + public static bool ConvolutionWithTranspose( + JxlImageF input, + ReadOnlySpan kernel, + JxlImageF output) + { + if (output.XSize != input.YSize) + { + return false; + } + + if (output.YSize != input.XSize) + { + return false; + } + + int len = kernel.Length; + int offset = len / 2; + + float weightNoBorder = 0.0f; + + for (int j = 0; j < len; j++) + { + weightNoBorder += kernel[j]; + } + + float scaleNoBorder = 1.0f / weightNoBorder; + + int border1 = Math.Min(input.XSize, offset); + int border2 = input.XSize > offset ? input.XSize - offset : 0; + + Span scaledKernel = stackalloc float[(len / 2) + 1]; + + for (int i = 0; i <= len / 2; i++) + { + scaledKernel[i] = kernel[i] * scaleNoBorder; + } + + // Middle + switch (len) + { + case 7: + { + float sk0 = scaledKernel[0]; + float sk1 = scaledKernel[1]; + float sk2 = scaledKernel[2]; + float sk3 = scaledKernel[3]; + + for (int y = 0; y < input.YSize; y++) + { + Span rowIn = input.GetRow(y); + + for (int x = border1; x < border2; x++) + { + int i = x - border1; + + float sum0 = (rowIn[i + 0] + rowIn[i + 6]) * sk0; + float sum1 = (rowIn[i + 1] + rowIn[i + 5]) * sk1; + float sum2 = (rowIn[i + 2] + rowIn[i + 4]) * sk2; + float sum = (rowIn[i + 3] * sk3) + sum0 + sum1 + sum2; + + output.GetRow(x)[y] = sum; + } + } + + break; + } + + case 13: + { + for (int y = 0; y < input.YSize; y++) + { + Span rowIn = input.GetRow(y); + + for (int x = border1; x < border2; x++) + { + int i = x - border1; + + float sum0 = (rowIn[i + 0] + rowIn[i + 12]) * scaledKernel[0]; + float sum1 = (rowIn[i + 1] + rowIn[i + 11]) * scaledKernel[1]; + float sum2 = (rowIn[i + 2] + rowIn[i + 10]) * scaledKernel[2]; + float sum3 = (rowIn[i + 3] + rowIn[i + 9]) * scaledKernel[3]; + + sum0 += (rowIn[i + 4] + rowIn[i + 8]) * scaledKernel[4]; + sum1 += (rowIn[i + 5] + rowIn[i + 7]) * scaledKernel[5]; + + float sum = rowIn[i + 6] * scaledKernel[6]; + + output.GetRow(x)[y] = sum + sum0 + sum1 + sum2 + sum3; + } + } + + break; + } + + case 15: + { + for (int y = 0; y < input.YSize; y++) + { + Span rowIn = input.GetRow(y); + + for (int x = border1; x < border2; x++) + { + int i = x - border1; + + float sum0 = (rowIn[i + 0] + rowIn[i + 14]) * scaledKernel[0]; + float sum1 = (rowIn[i + 1] + rowIn[i + 13]) * scaledKernel[1]; + float sum2 = (rowIn[i + 2] + rowIn[i + 12]) * scaledKernel[2]; + float sum3 = (rowIn[i + 3] + rowIn[i + 11]) * scaledKernel[3]; + + sum0 += (rowIn[i + 4] + rowIn[i + 10]) * scaledKernel[4]; + sum1 += (rowIn[i + 5] + rowIn[i + 9]) * scaledKernel[5]; + sum2 += (rowIn[i + 6] + rowIn[i + 8]) * scaledKernel[6]; + + float sum = rowIn[i + 7] * scaledKernel[7]; + + output.GetRow(x)[y] = sum + sum0 + sum1 + sum2 + sum3; + } + } + + break; + } + + case 33: + { + for (int y = 0; y < input.YSize; y++) + { + Span rowIn = input.GetRow(y); + + for (int x = border1; x < border2; x++) + { + int i = x - border1; + + float sum0 = (rowIn[i + 0] + rowIn[i + 32]) * scaledKernel[0]; + float sum1 = (rowIn[i + 1] + rowIn[i + 31]) * scaledKernel[1]; + float sum2 = (rowIn[i + 2] + rowIn[i + 30]) * scaledKernel[2]; + float sum3 = (rowIn[i + 3] + rowIn[i + 29]) * scaledKernel[3]; + + sum0 += (rowIn[i + 4] + rowIn[i + 28]) * scaledKernel[4]; + sum1 += (rowIn[i + 5] + rowIn[i + 27]) * scaledKernel[5]; + sum2 += (rowIn[i + 6] + rowIn[i + 26]) * scaledKernel[6]; + sum3 += (rowIn[i + 7] + rowIn[i + 25]) * scaledKernel[7]; + + sum0 += (rowIn[i + 8] + rowIn[i + 24]) * scaledKernel[8]; + sum1 += (rowIn[i + 9] + rowIn[i + 23]) * scaledKernel[9]; + sum2 += (rowIn[i + 10] + rowIn[i + 22]) * scaledKernel[10]; + sum3 += (rowIn[i + 11] + rowIn[i + 21]) * scaledKernel[11]; + + sum0 += (rowIn[i + 12] + rowIn[i + 20]) * scaledKernel[12]; + sum1 += (rowIn[i + 13] + rowIn[i + 19]) * scaledKernel[13]; + sum2 += (rowIn[i + 14] + rowIn[i + 18]) * scaledKernel[14]; + sum3 += (rowIn[i + 15] + rowIn[i + 17]) * scaledKernel[15]; + + float sum = rowIn[i + 16] * scaledKernel[16]; + + output.GetRow(x)[y] = sum + sum0 + sum1 + sum2 + sum3; + } + } + + break; + } + + default: + throw new NotSupportedException($"Kernel size {len} not implemented."); + } + + // Left border + for (int x = 0; x < border1; x++) + { + ConvolveBorderColumn(input, kernel, x, output.GetRow(x)); + } + + // Right border + for (int x = border2; x < input.XSize; x++) + { + ConvolveBorderColumn(input, kernel, x, output.GetRow(x)); + } + + return true; + } + + private static bool Blur( + JxlImageF input, + float sigma, + in ButteraugliParameters parameters, + ButteraugliBlurTemp temp, + JxlImageF output) + { + ReadOnlySpan kernel = ComputeKernel(sigma); + + // Separable5 does an in-place convolution, so this fast path is not safe + // if input aliases output. + if (kernel.Length == 5 && !ReferenceEquals(input, output)) + { + float sumWeights = 0.0f; + + foreach (float w in kernel) + { + sumWeights += w; + } + + float scale = 1.0f / sumWeights; + + float w0 = kernel[2] * scale; + float w1 = kernel[1] * scale; + float w2 = kernel[0] * scale; + + JxlWeightsSeparable5 weights = default; + FillRep4(ref weights.Horizontal, w0, w1, w2); + FillRep4(ref weights.Vertical, w0, w1, w2); + + if (!Separable5(input, input.GetRectangle(), weights, null, output)) + { + return false; + } + + return true; + } + + if (!temp.GetTransposed(input, out JxlImageF tempT)) + { + return false; + } + + if (!ConvolutionWithTranspose(input, kernel, tempT)) + { + return false; + } + + if (!ConvolutionWithTranspose(tempT, kernel, output)) + { + return false; + } + + return true; + } + + /// + /// Equivalent to HWY_REP4. + /// + private static void FillRep4(ref InlineArray12 values, float a, float b, float c) + { + for (int i = 0; i < 4; i++) + { + values[i] = a; + values[4 + i] = b; + values[8 + i] = c; + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector MaximumClamp(Vector value, float maximum) + { + Vector multiplier = new(0.724216145665f); + Vector maximumValue = new(maximum); + Vector ifPositive = ((value - maximumValue) * multiplier) + maximumValue; + Vector ifNegative = ((value + maximumValue) * multiplier) - maximumValue; + Vector positiveOrValue = Vector.ConditionalSelect(Vector.GreaterThan(value, maximumValue), ifPositive, value); + Vector result = Vector.ConditionalSelect(Vector.LessThan(value, Vector.Negate(maximumValue)), ifNegative, positiveOrValue); + + return result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector RemoveRangeAroundZero(Vector x, float width) + { + Vector w = new(width); + + return Vector.ConditionalSelect( + Vector.GreaterThan(x, w), + x - w, + Vector.ConditionalSelect( + Vector.LessThan(x, -w), + x + w, + Vector.Zero)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector AmplifyRangeAroundZero(Vector x, float width) + { + Vector w = new(width); + + return Vector.ConditionalSelect( + Vector.GreaterThan(x, w), + x + w, + Vector.ConditionalSelect( + Vector.LessThan(x, -w), + x - w, + x + x)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void XybLowFrequencyToValues( + Vector x, + Vector y, + Vector bArg, + out Vector valX, + out Vector valY, + out Vector valB) + { + Vector xMul = new(33.832837186260f); + Vector yMul = new(14.458268100570f); + Vector bMul = new(49.87984651440f); + Vector yToBMul = new(-0.362267051518f); + + Vector b = (yToBMul * y) + bArg; + + valB = b * bMul; + valX = x * xMul; + valY = y * yMul; + } + + public static void XybLowFrequencyToValues(JxlImage3F xybLf) + { + int lanes = Vector.Count; + + for (int y = 0; y < xybLf.YSize; y++) + { + Span rowX = xybLf.PlaneRow(0, y); + Span rowY = xybLf.PlaneRow(1, y); + Span rowB = xybLf.PlaneRow(2, y); + + for (int x = 0; x < xybLf.XSize; x += lanes) + { + Vector valX = new(rowX.Slice(x, lanes)); + Vector valY = new(rowY.Slice(x, lanes)); + Vector valB = new(rowB.Slice(x, lanes)); + + XybLowFrequencyToValues( + valX, + valY, + valB, + out valX, + out valY, + out valB); + + valX.CopyTo(rowX.Slice(x, lanes)); + valY.CopyTo(rowY.Slice(x, lanes)); + valB.CopyTo(rowB.Slice(x, lanes)); + } + } + } + + public static bool SuppressXByY(JxlImageF inY, JxlImageF inOutX) + { + if (!SameSize(inOutX, inY)) + { + return false; + } + + int xSize = inY.XSize; + int ySize = inY.YSize; + int lanes = Vector.Count; + + const float suppress = 46.0f; + const float s = 0.653020556257f; + + Vector sv = new(s); + Vector oneMinusS = new(1.0f - s); + Vector ywv = new(suppress); + + for (int y = 0; y < ySize; y++) + { + ReadOnlySpan rowY = inY.GetRow(y); + Span rowX = inOutX.GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + Vector vx = new(rowX.Slice(x, lanes)); + Vector vy = new(rowY.Slice(x, lanes)); + + Vector scaler = + ((ywv / ((vy * vy) + ywv)) * oneMinusS) + sv; + + (scaler * vx).CopyTo(rowX.Slice(x, lanes)); + } + } + + return true; + } + + public static void Subtract(JxlPlane a, JxlPlane b, JxlPlane c) + { + int lanes = Vector.Count; + + for (int y = 0; y < a.YSize; y++) + { + ReadOnlySpan rowA = a.GetRow(y); + ReadOnlySpan rowB = b.GetRow(y); + Span rowC = c.GetRow(y); + + for (int x = 0; x < a.XSize; x += lanes) + { + Vector va = new(rowA.Slice(x, lanes)); + Vector vb = new(rowB.Slice(x, lanes)); + + (va - vb).CopyTo(rowC.Slice(x, lanes)); + } + } + } + + public static bool SeparateLFAndMF( + in ButteraugliParameters parameters, + JxlImage3F xyb, + JxlImage3F lf, + JxlImage3F mf, + ButteraugliBlurTemp blurTemp) + { + const float sigmaLf = 7.15593339443f; + + for (int i = 0; i < 3; i++) + { + if (!Blur( + xyb.Plane(i), + sigmaLf, + parameters, + blurTemp, + lf.Plane(i))) + { + return false; + } + + Subtract( + xyb.Plane(i), + lf.Plane(i), + mf.Plane(i)); + } + + XybLowFrequencyToValues(lf); + + return true; + } + + public static bool SeparateMfAndHf( + Configuration configuration, + in ButteraugliParameters parameters, + JxlImage3F mf, + JxlImageF[] hf, + BlurTemp 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); + + int lanes = Vector.Count; + + for (int i = 0; i < 3; i++) + { + if (i == 2) + { + if (!Blur(mf.Plane(i), sigmaHf, parameters, blurTemp, mf.Plane(i))) + { + return false; + } + + break; + } + + for (int y = 0; y < ySize; y++) + { + Span rowMf = mf.PlaneRow(i, y); + Span rowHf = hf[i].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + new Vector(rowMf.Slice(x, lanes)) + .CopyTo(rowHf.Slice(x, lanes)); + } + } + + if (!Blur(mf.Plane(i), sigmaHf, parameters, blurTemp, mf.Plane(i))) + { + return false; + } + + if (i == 0) + { + const float removeMfRange = 0.29f; + + for (int y = 0; y < ySize; y++) + { + Span rowMf = mf.PlaneRow(0, y); + Span rowHf = hf[0].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + Vector mfv = new(rowMf.Slice(x, lanes)); + Vector hfv = new Vector(rowHf.Slice(x, lanes)) - mfv; + + mfv = RemoveRangeAroundZero(mfv, removeMfRange); + + mfv.CopyTo(rowMf.Slice(x, lanes)); + hfv.CopyTo(rowHf.Slice(x, lanes)); + } + } + } + else + { + const float addMfRange = 0.1f; + + for (int y = 0; y < ySize; y++) + { + Span rowMf = mf.PlaneRow(1, y); + Span rowHf = hf[1].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + Vector mfv = new(rowMf.Slice(x, lanes)); + Vector hfv = new Vector(rowHf.Slice(x, lanes)) - mfv; + + mfv = AmplifyRangeAroundZero(mfv, addMfRange); + + mfv.CopyTo(rowMf.Slice(x, lanes)); + hfv.CopyTo(rowHf.Slice(x, lanes)); + } + } + } + } + + return SuppressXByY(hf[1], hf[0]); + } + + public static bool SeparateHFAndUHF( + in ButteraugliParameters parameters, + JxlImageF[] hf, + JxlImageF[] uhf, + JxlBlurTemp 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); + + int lanes = Vector.Count; + + for (int i = 0; i < 2; i++) + { + for (int y = 0; y < ySize; y++) + { + Span rowUhf = uhf[i].GetRow(y); + Span rowHf = hf[i].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + new Vector(rowHf.Slice(x, lanes)).CopyTo(rowUhf.Slice(x, lanes)); + } + } + + if (!Blur(hf[i], sigmaUhf, parameters, blurTemp, hf[i])) + { + return false; + } + + if (i == 0) + { + const float removeHfRange = 1.5f; + const float removeUhfRange = 0.04f; + + for (int y = 0; y < ySize; y++) + { + Span rowUhf = uhf[0].GetRow(y); + Span rowHf = hf[0].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + Vector hfv = new(rowHf.Slice(x, lanes)); + Vector uhfv = new Vector(rowUhf.Slice(x, lanes)) - hfv; + + hfv = RemoveRangeAroundZero(hfv, removeHfRange); + uhfv = RemoveRangeAroundZero(uhfv, removeUhfRange); + + hfv.CopyTo(rowHf.Slice(x, lanes)); + uhfv.CopyTo(rowUhf.Slice(x, lanes)); + } + } + } + else + { + const float addHfRange = 0.132f; + const float maxClampHf = 28.4691806922f; + const float maxClampUhf = 5.19175294647f; + const float mulYHf = 2.155f; + const float mulYUhf = 2.69313763794f; + + Vector mulHf = new(mulYHf); + Vector mulUhf = new(mulYUhf); + + for (int y = 0; y < ySize; y++) + { + Span rowUhf = uhf[1].GetRow(y); + Span rowHf = hf[1].GetRow(y); + + for (int x = 0; x < xSize; x += lanes) + { + Vector hfv = new(rowHf.Slice(x, lanes)); + hfv = MaximumClamp(hfv, maxClampHf); + + Vector uhfv = new Vector(rowUhf.Slice(x, lanes)) - hfv; + + uhfv = MaximumClamp(uhfv, maxClampUhf); + uhfv *= mulUhf; + + uhfv.CopyTo(rowUhf.Slice(x, lanes)); + + hfv *= mulHf; + hfv = AmplifyRangeAroundZero(hfv, addHfRange); + + hfv.CopyTo(rowHf.Slice(x, lanes)); + } + } + } + } + + return true; + } + + public static void DeallocateHFAndUHF(JxlImageF[] hf, JxlImageF[] uhf) + { + for (int i = 0; i < 2; i++) + { + hf[i] = new JxlImageF(); + uhf[i] = new JxlImageF(); + } + } + + public static bool SeparateFrequencies( + Configuration configuration, + in ButteraugliParameters parameters, + BlurTemp blurTemp, + JxlImage3F xyb, + PsychoImage ps) + { + ps.Lf = JxlImage3F.Create( + configuration, + xyb.XSize, + xyb.YSize); + + ps.Mf = JxlImage3F.Create( + configuration, + xyb.XSize, + xyb.YSize); + + if (!SeparateLFAndMF( + parameters, + xyb, + ps.Lf, + ps.Mf, + blurTemp)) + { + return false; + } + + if (!SeparateMfAndHf( + parameters, + ps.Mf, + ps.Hf, + blurTemp)) + { + return false; + } + + if (!SeparateHFAndUHF( + parameters, + ps.Hf, + ps.Uhf, + blurTemp)) + { + return false; + } + + return true; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector Sum( + Vector a, + Vector b, + Vector c, + Vector d) + => (a + b) + (c + d); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector Sum( + Vector a, + Vector b, + Vector c, + Vector d, + Vector e) + => Sum(a, b, c, d + e); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector Sum( + Vector a, + Vector b, + Vector c, + Vector d, + Vector e, + Vector f, + Vector g) + => Sum(a, b, c, Sum(d, e, f, g)); + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector Sum( + Vector a, + Vector b, + Vector c, + Vector d, + Vector e, + Vector f, + Vector g, + Vector h, + Vector i) + => (Sum(a, b, c, d) + Sum(e, f, g, h)) + i; + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector MaltaUnitLF( + ReadOnlySpan row, + int index, + int xs) + { + // helper + [MethodImpl(MethodImplOptions.AggressiveInlining)] + static Vector Load(ReadOnlySpan row, int index) + => new(row.Slice(index, Vector.Count)); + + int xs3 = 3 * xs; + + Vector center = Load(row, index); + + Vector sumYConst = Sum( + Load(row, index - 4), + Load(row, index - 2), + center, + Load(row, index + 2), + Load(row, index + 4)); + + Vector retval = sumYConst * sumYConst; + + Vector sum; + + sum = Sum( + Load(row, index - xs3 - xs), + Load(row, index - xs - xs), + center, + Load(row, index + xs + xs), + Load(row, index + xs3 + xs)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - 3), + Load(row, index - xs - xs - 2), + center, + Load(row, index + xs + xs + 2), + Load(row, index + xs3 + 3)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 + 3), + Load(row, index - xs - xs + 2), + center, + Load(row, index + xs + xs - 2), + Load(row, index + xs3 - 3)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - xs + 1), + Load(row, index - xs - xs + 1), + center, + Load(row, index + xs + xs - 1), + Load(row, index + xs3 + xs - 1)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - xs - 1), + Load(row, index - xs - xs - 1), + center, + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + xs + 1)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - 4 - xs), + Load(row, index - 2 - xs), + center, + Load(row, index + 2 + xs), + Load(row, index + 4 + xs)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - 4 + xs), + Load(row, index - 2 + xs), + center, + Load(row, index + 2 - xs), + Load(row, index + 4 - xs)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - 2), + Load(row, index - xs - xs - 1), + center, + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + 2)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 + 2), + Load(row, index - xs - xs + 1), + center, + Load(row, index + xs + xs - 1), + Load(row, index + xs3 - 2)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs - xs - 3), + Load(row, index - xs - 2), + center, + Load(row, index + xs + 2), + Load(row, index + xs + xs + 3)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs - xs + 3), + Load(row, index - xs + 2), + center, + Load(row, index + xs - 2), + Load(row, index + xs + xs - 3)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index + xs + xs - 4), + Load(row, index + xs - 2), + center, + Load(row, index - xs + 2), + Load(row, index - xs - xs + 4)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs - xs - 4), + Load(row, index - xs - 2), + center, + Load(row, index + xs + 2), + Load(row, index + xs + xs + 4)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - xs - 2), + Load(row, index - xs - xs - 1), + center, + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + xs + 2)); + retval = (sum * sum) + retval; + + sum = Sum( + Load(row, index - xs3 - xs + 2), + Load(row, index - xs - xs + 1), + center, + Load(row, index + xs + xs - 1), + Load(row, index + xs3 + xs - 2)); + retval = (sum * sum) + retval; + + return retval; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Vector MaltaUnit( + ReadOnlySpan row, + int index, + int xs) + { + // helper + [MethodImpl(MethodImplOptions.AggressiveInlining)] + static Vector Load(ReadOnlySpan row, int index) + => new(row.Slice(index, Vector.Count)); + + int xs3 = 3 * xs; + + Vector center = Load(row, index); + + Vector sumYConst = Sum( + Load(row, index - 4), + Load(row, index - 3), + Load(row, index - 2), + Load(row, index - 1), + center, + Load(row, index + 1), + Load(row, index + 2), + Load(row, index + 3), + Load(row, index + 4)); + + Vector retval = sumYConst * sumYConst; + + Vector sum; + + sum = Sum( + Load(row, index - xs3 - xs), + Load(row, index - xs3), + Load(row, index - xs - xs), + Load(row, index - xs), + center, + Load(row, index + xs), + Load(row, index + xs + xs), + Load(row, index + xs3), + Load(row, index + xs3 + xs)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - 3), + Load(row, index - xs - xs - 2), + Load(row, index - xs - 1), + center, + Load(row, index + xs + 1), + Load(row, index + xs + xs + 2), + Load(row, index + xs3 + 3)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 + 3), + Load(row, index - xs - xs + 2), + Load(row, index - xs + 1), + center, + Load(row, index + xs - 1), + Load(row, index + xs + xs - 2), + Load(row, index + xs3 - 3)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - xs + 1), + Load(row, index - xs3 + 1), + Load(row, index - xs - xs + 1), + Load(row, index - xs), + center, + Load(row, index + xs), + Load(row, index + xs + xs - 1), + Load(row, index + xs3 - 1), + Load(row, index + xs3 + xs - 1)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - xs - 1), + Load(row, index - xs3 - 1), + Load(row, index - xs - xs - 1), + Load(row, index - xs), + center, + Load(row, index + xs), + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + 1), + Load(row, index + xs3 + xs + 1)); + retval += sum * sum; + + sum = Sum( + Load(row, index - 4 - xs), + Load(row, index - 3 - xs), + Load(row, index - 2 - xs), + Load(row, index - 1), + center, + Load(row, index + 1), + Load(row, index + 2 + xs), + Load(row, index + 3 + xs), + Load(row, index + 4 + xs)); + retval += sum * sum; + + sum = Sum( + Load(row, index - 4 + xs), + Load(row, index - 3 + xs), + Load(row, index - 2 + xs), + Load(row, index - 1), + center, + Load(row, index + 1), + Load(row, index + 2 - xs), + Load(row, index + 3 - xs), + Load(row, index + 4 - xs)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - 2), + Load(row, index - xs - xs - 1), + Load(row, index - xs - 1), + center, + Load(row, index + xs + 1), + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + 2)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 + 2), + Load(row, index - xs - xs + 1), + Load(row, index - xs + 1), + center, + Load(row, index + xs - 1), + Load(row, index + xs + xs - 1), + Load(row, index + xs3 - 2)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs - xs - 3), + Load(row, index - xs - 2), + Load(row, index - xs - 1), + center, + Load(row, index + xs + 1), + Load(row, index + xs + 2), + Load(row, index + xs + xs + 3)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs - xs + 3), + Load(row, index - xs + 2), + Load(row, index - xs + 1), + center, + Load(row, index + xs - 1), + Load(row, index + xs - 2), + Load(row, index + xs + xs - 3)); + retval += sum * sum; + + sum = Sum( + Load(row, index + xs - 4), + Load(row, index + xs - 3), + Load(row, index + xs - 2), + Load(row, index - 1), + center, + Load(row, index + 1), + Load(row, index - xs + 2), + Load(row, index - xs + 3), + Load(row, index - xs + 4)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs - 4), + Load(row, index - xs - 3), + Load(row, index - xs - 2), + Load(row, index - 1), + center, + Load(row, index + 1), + Load(row, index + xs + 2), + Load(row, index + xs + 3), + Load(row, index + xs + 4)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - xs - 1), + Load(row, index - xs3 - 1), + Load(row, index - xs - xs - 1), + Load(row, index - xs), + center, + Load(row, index + xs), + Load(row, index + xs + xs + 1), + Load(row, index + xs3 + 1), + Load(row, index + xs3 + xs + 1)); + retval += sum * sum; + + sum = Sum( + Load(row, index - xs3 - xs + 1), + Load(row, index - xs3 + 1), + Load(row, index - xs - xs + 1), + Load(row, index - xs), + center, + Load(row, index + xs), + Load(row, index + xs + xs - 1), + Load(row, index + xs3 - 1), + Load(row, index + xs3 + xs - 1)); + retval += sum * sum; + + return retval; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static float PaddedMaltaUnit( + JxlImageF diffs, + int x0, + int y0, + bool isLF) + { + if (x0 >= 4 && + y0 >= 4 && + x0 < diffs.XSize - 4 && + y0 < diffs.YSize - 4) + { + return isLF + ? MaltaUnitLF( + diffs.GetRow(y0), + x0, + diffs.PixelsPerRow)[0] + : MaltaUnit( + diffs.GetRow(y0), + x0, + diffs.PixelsPerRow)[0]; + } + + Span borderImage = stackalloc float[12 * 9]; + + for (int dy = 0; dy < 9; dy++) + { + int y = y0 + dy - 4; + + if (y < 0 || y >= diffs.YSize) + { + borderImage.Slice(dy * 12, 12).Clear(); + continue; + } + + ReadOnlySpan rowDiffs = diffs.GetRow(y); + + for (int dx = 0; dx < 9; dx++) + { + int x = x0 + dx - 4; + + borderImage[(dy * 12) + dx] = + x < 0 || x >= diffs.XSize + ? 0.0f + : rowDiffs[x]; + } + + borderImage.Slice((dy * 12) + 9, 3).Clear(); + } + + return isLF + ? MaltaUnitLF( + diffs.GetRow(y0), + x0, + diffs.PixelsPerRow)[0] + : MaltaUnit( + borderImage, + (4 * 12) + 4, + 12)[0]; + } + + public static bool MaltaDiffMap( + bool isLf, + JxlImageF lum0, + JxlImageF lum1, + float w0Gt1, + float w0Lt1, + float norm1, + float len, + float mulli, + JxlImageF diffs, + JxlImageF blockDiffAc) + { + if (!SameSize(lum0, lum1) || !SameSize(lum0, diffs)) + { + return false; + } + + int width = lum0.XSize; + int height = lum0.YSize; + + const float weight0 = 0.5f; + const float weight1 = 0.33f; + + float norm2_0Gt1 = + (float)(mulli * MathF.Sqrt(weight0 * w0Gt1) / ((len * 2) + 1) * norm1); + + float norm2_0Lt1 = + (float)(mulli * MathF.Sqrt(weight1 * w0Lt1) / ((len * 2) + 1) * norm1); + + for (int y = 0; y < height; y++) + { + ReadOnlySpan row0 = lum0.GetRow(y); + ReadOnlySpan row1 = lum1.GetRow(y); + Span rowDiffs = diffs.GetRow(y); + + for (int x = 0; x < width; x++) + { + float absVal = 0.5f * + (MathF.Abs(row0[x]) + MathF.Abs(row1[x])); + + float diff = row0[x] - row1[x]; + + float scaler = norm2_0Gt1 / ((float)norm1 + absVal); + + rowDiffs[x] = scaler * diff; + + float scaler2 = norm2_0Lt1 / ((float)norm1 + absVal); + + float fabs0 = MathF.Abs(row0[x]); + + float tooSmall = 0.55f * fabs0; + float tooBig = 1.05f * fabs0; + + if (row0[x] < 0) + { + if (row1[x] > -tooSmall) + { + rowDiffs[x] -= scaler2 * (row1[x] + tooSmall); + } + else if (row1[x] < -tooBig) + { + rowDiffs[x] += scaler2 * (-row1[x] - tooBig); + } + } + else + { + if (row1[x] < tooSmall) + { + rowDiffs[x] += scaler2 * (tooSmall - row1[x]); + } + else if (row1[x] > tooBig) + { + rowDiffs[x] -= scaler2 * (row1[x] - tooBig); + } + } + } + } + + int y0 = 0; + + for (; y0 < 4; y0++) + { + Span row = blockDiffAc.GetRow(y0); + + for (int x = 0; x < width; x++) + { + row[x] += PaddedMaltaUnit(diffs, x, y0, isLf); + } + } + + int lanes = Vector.Count; + int alignedX = Math.Max(4, lanes); + + int stride = diffs.PixelsPerRow; + + for (; y0 < height - 4; y0++) + { + ReadOnlySpan input = diffs.GetRow(y0); + Span output = blockDiffAc.GetRow(y0); + ref float outputReference = ref MemoryMarshal.GetReference(output); + + int x = 0; + + for (; x < alignedX; x++) + { + output[x] += PaddedMaltaUnit(diffs, x, y0, isLf); + } + + for (; x + lanes + 4 <= width; x += lanes) + { + Vector value = Vector.LoadUnsafe(ref Unsafe.Add(ref outputReference, x)); + + Vector malta = isLf + ? MaltaUnitLF(input, x, stride) + : MaltaUnit(input, x, stride); + + (value + malta).CopyTo(output[x..]); + } + + for (; x < width; x++) + { + output[x] += PaddedMaltaUnit(diffs, x, y0, isLf); + } + } + + for (; y0 < height; y0++) + { + Span row = blockDiffAc.GetRow(y0); + + for (int x = 0; x < width; x++) + { + row[x] += PaddedMaltaUnit(diffs, x, y0, isLf); + } + } + + return true; + } + + public static bool MaltaDiffMap( + JxlImageF lum0, + JxlImageF lum1, + float w0Gt1, + float w0Lt1, + float norm1, + JxlImageF diffs, + JxlImageF blockDiffAc) + { + const float len = 3.75f; + const float mulli = 0.39905817637f; + + return MaltaDiffMap( + false, + lum0, + lum1, + w0Gt1, + w0Lt1, + norm1, + len, + mulli, + diffs, + blockDiffAc); + } + + public static bool MaltaDiffMapLf( + JxlImageF lum0, + JxlImageF lum1, + float w0Gt1, + float w0Lt1, + float norm1, + JxlImageF diffs, + JxlImageF blockDiffAc) + { + const float len = 3.75f; + const float mulli = 0.611612573796f; + + return MaltaDiffMap( + true, + lum0, + lum1, + w0Gt1, + w0Lt1, + norm1, + len, + mulli, + diffs, + blockDiffAc); + } + + public static void CombineChannelsForMasking( + JxlImageF[] hf, + JxlImageF[] uhf, + JxlImageF output) + { + // Only X and Y components are involved in masking. + ReadOnlySpan muls = + [ + 2.5f, + 0.4f, + 0.4f + ]; + + int width = hf[0].XSize; + int height = hf[0].YSize; + + for (int y = 0; y < height; y++) + { + ReadOnlySpan rowYHf = hf[1].GetRow(y); + ReadOnlySpan rowYUhf = uhf[1].GetRow(y); + ReadOnlySpan rowXHf = hf[0].GetRow(y); + ReadOnlySpan rowXUhf = uhf[0].GetRow(y); + + Span row = output.GetRow(y); + + for (int x = 0; x < width; x++) + { + float xDiff = (rowXUhf[x] + rowXHf[x]) * muls[0]; + float yDiff = (rowYUhf[x] * muls[1]) + (rowYHf[x] * muls[2]); + + row[x] = MathF.Sqrt((xDiff * xDiff) + (yDiff * yDiff)); + } + } + } + + public static void DiffPrecompute( + JxlImageF xyb, + float mul, + float biasArg, + JxlImageF output) + { + int width = xyb.XSize; + int height = xyb.YSize; + + float bias = mul * biasArg; + float sqrtBias = MathF.Sqrt(bias); + + for (int y = 0; y < height; y++) + { + ReadOnlySpan rowIn = xyb.GetRow(y); + Span rowOut = output.GetRow(y); + + for (int x = 0; x < width; x++) + { + rowOut[x] = + MathF.Sqrt( + (mul * MathF.Abs(rowIn[x])) + bias) + - sqrtBias; + } + } + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static void StoreMin3( + float value, + ref float min0, + ref float min1, + ref float min2) + { + if (value < min2) + { + if (value < min0) + { + min2 = min1; + min1 = min0; + min0 = value; + } + else if (value < min1) + { + min2 = min1; + min1 = value; + } + else + { + min2 = value; + } + } + } + + public static void FuzzyErosion(JxlImageF from, JxlImageF to) + { + int width = from.XSize; + int height = from.YSize; + + const int step = 3; + + for (int y = 0; y < height; y++) + { + Span output = to.GetRow(y); + + for (int x = 0; x < width; x++) + { + float min0 = from.GetRow(y)[x]; + float min1 = 2 * min0; + float min2 = min1; + + if (x >= step) + { + StoreMin3( + from.GetRow(y)[x - step], + ref min0, + ref min1, + ref min2); + + if (y >= step) + { + StoreMin3( + from.GetRow(y - step)[x - step], + ref min0, + ref min1, + ref min2); + } + + if (y < height - step) + { + StoreMin3( + from.GetRow(y + step)[x - step], + ref min0, + ref min1, + ref min2); + } + } + + if (x < width - step) + { + StoreMin3( + from.GetRow(y)[x + step], + ref min0, + ref min1, + ref min2); + + if (y >= step) + { + StoreMin3( + from.GetRow(y - step)[x + step], + ref min0, + ref min1, + ref min2); + } + + if (y < height - step) + { + StoreMin3( + from.GetRow(y + step)[x + step], + ref min0, + ref min1, + ref min2); + } + } + + if (y >= step) + { + StoreMin3( + from.GetRow(y - step)[x], + ref min0, + ref min1, + ref min2); + } + + if (y < height - step) + { + StoreMin3( + from.GetRow(y + step)[x], + ref min0, + ref min1, + ref min2); + } + + output[x] = + (0.45f * min0) + + (0.3f * min1) + + (0.25f * min2); + } + } + } + + public static bool Mask( + Configuration configuration, + JxlImageF mask0, + JxlImageF mask1, + in ButteraugliParameters parameters, + JxlBlurTemp blurTemp, + JxlImageF? diffAc, + out JxlImageF mask) + { + int width = mask0.XSize; + int height = mask0.YSize; + + mask = new(configuration, width, height); + + const float mul = 6.19424080439f; + const float bias = 12.61050594197f; + const float radius = 2.7f; + + JxlImageF diff0 = new(configuration, width, height); + JxlImageF diff1 = new(configuration, width, height); + JxlImageF blurred0 = new(configuration, width, height); + JxlImageF blurred1 = new(configuration, width, height); + + DiffPrecompute(mask0, mul, bias, diff0); + DiffPrecompute(mask1, mul, bias, diff1); + + if (!Blur(diff0, radius, parameters, blurTemp, blurred0)) + { + return false; + } + + FuzzyErosion(blurred0, diff0); + + if (!Blur(diff1, radius, parameters, blurTemp, blurred1)) + { + return false; + } + + for (int y = 0; y < height; y++) + { + Span maskRow = mask.GetRow(y); + Span diffRow = diffAc is not null ? diffAc.GetRow(y) : []; + + ReadOnlySpan diff0Row = diff0.GetRow(y); + ReadOnlySpan blur0Row = blurred0.GetRow(y); + ReadOnlySpan blur1Row = blurred1.GetRow(y); + + for (int x = 0; x < width; x++) + { + maskRow[x] = diff0Row[x]; + + if (diffRow != null) + { + const float maskToErrorMul = 10.0f; + float diff = blur0Row[x] - blur1Row[x]; + diffRow[x] += maskToErrorMul * diff * diff; + } + } + } + + return true; + } + + public static bool MaskPsychoImage( + Configuration configuration, + ButteraugliPsychoImage pi0, + ButteraugliPsychoImage pi1, + int width, + int height, + in ButteraugliParameters parameters, + BlurTemp blurTemp, + JxlImageF mask, + JxlImageF? diffAc) + { + JxlImageF mask0 = new(configuration, width, height); + JxlImageF mask1 = new(configuration, width, height); + + CombineChannelsForMasking( + pi0.Hf, + pi0.Uhf, + mask0); + + CombineChannelsForMasking( + pi1.Hf, + pi1.Uhf, + mask1); + + return Mask( + mask0, + mask1, + parameters, + blurTemp, + mask, + diffAc); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static float MaskY(float delta) + { + const float offset = 0.829591754942f; + const float scaler = 0.451936922203f; + const float mul = 2.5485944793f; + + float c = mul / ((scaler * delta) + offset); + float result = GlobalScale * (1.0f + c); + + return result * result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static float MaskDcY(float delta) + { + const float offset = 0.20025578522f; + const float scaler = 3.87449418804f; + const float mul = 0.505054525019f; + + float c = mul / ((scaler * delta) + offset); + float result = GlobalScale * (1.0f + c); + + return result * result; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static float MaskColor( + ReadOnlySpan color, + float mask) + => (color[0] * mask) + (color[1] * mask) + (color[2] * mask); + + public static bool CombineChannelsToDiffmap( + JxlImageF mask, + JxlImage3F blockDiffDc, + JxlImage3F blockDiffAc, + float xmul, + JxlImageF result) + { + if (!SameSize(mask, result)) + { + return false; + } + + Span diffDc = stackalloc float[3]; + Span diffAc = stackalloc float[3]; + + int xsize = mask.XSize; + int ysize = mask.YSize; + + for (int y = 0; y < ysize; ++y) + { + Span rowOut = result.GetRow(y); + + for (int x = 0; x < xsize; ++x) + { + float val = mask.GetRow(y)[x]; + + float maskVal = MaskY(val); + float dcMaskVal = MaskDcY(val); + + for (int i = 0; i < 3; ++i) + { + diffDc[i] = blockDiffDc.PlaneRow(i, y)[x]; + diffAc[i] = blockDiffAc.PlaneRow(i, y)[x]; + } + + diffAc[0] *= xmul; + diffDc[0] *= xmul; + + rowOut[x] = MathF.Sqrt( + MaskColor(diffDc, dcMaskVal) + + MaskColor(diffAc, maskVal)); + } + } + + return true; + } + + public static void L2Diff( + JxlImageF i0, + JxlImageF i1, + float w, + JxlImageF diffmap) + { + if (w == 0) + { + return; + } + + for (int y = 0; y < i0.YSize; ++y) + { + ReadOnlySpan row0 = i0.GetRow(y); + ReadOnlySpan row1 = i1.GetRow(y); + Span rowDiff = diffmap.GetRow(y); + + for (int x = 0; x < i0.XSize; ++x) + { + float diff = row0[x] - row1[x]; + rowDiff[x] += diff * diff * w; + } + } + } + + public static void SetL2Diff( + JxlImageF i0, + JxlImageF i1, + float w, + JxlImageF diffmap) + { + if (w == 0) + { + return; + } + + for (int y = 0; y < i0.YSize; ++y) + { + ReadOnlySpan row0 = i0.GetRow(y); + ReadOnlySpan row1 = i1.GetRow(y); + Span rowDiff = diffmap.GetRow(y); + + for (int x = 0; x < i0.XSize; ++x) + { + float diff = row0[x] - row1[x]; + rowDiff[x] = diff * diff * w; + } + } + } + + public static void L2DiffAsymmetric( + JxlImageF i0, + JxlImageF i1, + float w0gt1, + float w0lt1, + JxlImageF diffmap) + { + if (w0gt1 == 0 && w0lt1 == 0) + { + return; + } + + Vector vw0gt1 = new(w0gt1 * 0.8f); + Vector vw0lt1 = new(w0lt1 * 0.8f); + + int lanes = Vector.Count; + + for (int y = 0; y < i0.YSize; ++y) + { + ReadOnlySpan row0 = i0.GetRow(y); + ReadOnlySpan row1 = i1.GetRow(y); + Span rowDiff = diffmap.GetRow(y); + + for (int x = 0; x < i0.XSize; x += lanes) + { + Vector val0 = new(row0[x..]); + Vector val1 = new(row1[x..]); + + // Primary symmetric quadratic objective. + Vector diff = val0 - val1; + + Vector total = + (diff * diff * vw0gt1) + new Vector(rowDiff[x..]); + + Vector fabs0 = Vector.Abs(val0); + + Vector tooSmall = fabs0 * new Vector(0.4f); + + Vector tooBig = fabs0; + + Vector ifNeg = + Vector.ConditionalSelect( + Vector.GreaterThan(val1, -tooSmall), + val1 + tooSmall, + Vector.ConditionalSelect( + Vector.LessThan(val1, -tooBig), + -val1 - tooBig, + Vector.Zero)); + + Vector ifPos = + Vector.ConditionalSelect( + Vector.LessThan(val1, tooSmall), + tooSmall - val1, + Vector.ConditionalSelect( + Vector.GreaterThan(val1, tooBig), + val1 - tooBig, + Vector.Zero)); + + Vector v = + Vector.ConditionalSelect( + Vector.LessThan(val0, Vector.Zero), + ifNeg, + ifPos); + + total += vw0lt1 * v * v; + + total.CopyTo(rowDiff[x..]); + } + } + } + + public static void OpsinAbsorbance( + bool clamp, + Vector in0, + Vector in1, + Vector in2, + out Vector out0, + out Vector out1, + out Vector out2) + { + Vector mix0 = new(0.29956550340058319f); + Vector mix1 = new(0.63373087833825936f); + Vector mix2 = new(0.077705617820981968f); + Vector mix3 = new(1.7557483643287353f); + + Vector mix4 = new(0.22158691104574774f); + Vector mix5 = new(0.69391388044116142f); + Vector mix6 = new(0.0987313588422f); + Vector mix7 = new(1.7557483643287353f); + + Vector mix8 = new(0.02f); + Vector mix9 = new(0.02f); + Vector mix10 = new(0.20480129041026129f); + Vector mix11 = new(12.226454707163354f); + + out0 = (mix0 * in0) + ((mix1 * in1) + ((mix2 * in2) + mix3)); + out1 = (mix4 * in0) + ((mix5 * in1) + ((mix6 * in2) + mix7)); + out2 = (mix8 * in0) + ((mix9 * in1) + ((mix10 * in2) + mix11)); + + if (clamp) + { + out0 = Vector.Max(out0, mix3); + out1 = Vector.Max(out1, mix7); + out2 = Vector.Max(out2, mix11); + } + } + + public static bool OpsinDynamicsImage( + JxlImage3F rgb, + in ButteraugliParameters parameters, + JxlImage3F blurred, + BlurTemp blurTemp, + JxlImage3F xyb) + { + if (blurred == null) + { + return false; + } + + const double sigma = 1.2; + + if (!Blur(rgb.Plane(0), sigma, parameters, blurTemp, blurred.Plane(0))) + { + return false; + } + + if (!Blur(rgb.Plane(1), sigma, parameters, blurTemp, blurred.Plane(1))) + { + return false; + } + + if (!Blur(rgb.Plane(2), sigma, parameters, blurTemp, blurred.Plane(2))) + { + return false; + } + + Vector intensityMultiplier = new((float)parameters.IntensityTarget); + + int lanes = Vector.Count; + + Vector minValue = new(1e-4f); + + for (int y = 0; y < rgb.YSize; ++y) + { + ReadOnlySpan rowR = rgb.PlaneRow(0, y); + ReadOnlySpan rowG = rgb.PlaneRow(1, y); + ReadOnlySpan rowB = rgb.PlaneRow(2, y); + + ReadOnlySpan blurredR = blurred.PlaneRow(0, y); + ReadOnlySpan blurredG = blurred.PlaneRow(1, y); + ReadOnlySpan blurredB = blurred.PlaneRow(2, y); + + Span outX = xyb.PlaneRow(0, y); + Span outY = xyb.PlaneRow(1, y); + Span outB = xyb.PlaneRow(2, y); + + for (int x = 0; x < rgb.XSize; x += lanes) + { + Vector sensitivity0; + Vector sensitivity1; + Vector sensitivity2; + { + OpsinAbsorbance( + true, + new Vector(blurredR[x..]) * intensityMultiplier, + new Vector(blurredG[x..]) * intensityMultiplier, + new Vector(blurredB[x..]) * intensityMultiplier, + out Vector pre0, + out Vector pre1, + out Vector pre2); + + pre0 = Vector.Max(pre0, minValue); + pre1 = Vector.Max(pre1, minValue); + pre2 = Vector.Max(pre2, minValue); + + sensitivity0 = Gamma(pre0) / pre0; + sensitivity1 = Gamma(pre1) / pre1; + sensitivity2 = Gamma(pre2) / pre2; + + sensitivity0 = Vector.Max(sensitivity0, minValue); + sensitivity1 = Vector.Max(sensitivity1, minValue); + sensitivity2 = Vector.Max(sensitivity2, minValue); + } + + OpsinAbsorbance( + false, + new Vector(rowR[x..]) * intensityMultiplier, + new Vector(rowG[x..]) * intensityMultiplier, + new Vector(rowB[x..]) * intensityMultiplier, + out Vector cur0, + out Vector cur1, + out Vector cur2); + + cur0 *= sensitivity0; + cur1 *= sensitivity1; + cur2 *= sensitivity2; + + Vector min01 = new(1.7557483643287353f); + Vector min2 = new(12.226454707163354f); + + cur0 = Vector.Max(cur0, min01); + cur1 = Vector.Max(cur1, min01); + cur2 = Vector.Max(cur2, min2); + + (cur0 - cur1).CopyTo(outX[x..]); + (cur0 + cur1).CopyTo(outY[x..]); + cur2.CopyTo(outB[x..]); + } + } + + return true; + } + + public static bool ButteraugliDiffmapInPlace( + Configuration configuration, + JxlImage3F image0, + JxlImage3F image1, + in ButteraugliParameters parameters, + JxlImageF diffmap) + { + int xSize = image0.XSize; + int ySize = image0.YSize; + + using var blurTemp = new JxlBlurTemp(); + + using (JxlImage3F temp = new(configuration, xSize, ySize)) + { + if (!OpsinDynamicsImage(image0, parameters, temp, blurTemp, image0)) + { + return false; + } + + if (!OpsinDynamicsImage(image1, parameters, temp, blurTemp, image1)) + { + return false; + } + } + + using JxlPlane blockDiffDc = JxlImageF.Create(configuration, xSize, ySize); + blockDiffDc.ZeroFill(); + + // LF/DC + using (JxlImage3F lf0 = new(configuration, xSize, ySize)) + using (JxlImage3F lf1 = new(configuration, xSize, ySize)) + { + if (!SeparateLFAndMF(parameters, image0, lf0, image0, blurTemp)) + { + return false; + } + + if (!SeparateLFAndMF(parameters, image1, lf1, image1, blurTemp)) + { + return false; + } + + for (int c = 0; c < 3; c++) + { + L2Diff( + lf0.Plane(c), + lf1.Plane(c), + Wmul[6 + c], + blockDiffDc); + } + } + + JxlImageF[] hf0 = new JxlImageF[2]; + JxlImageF[] hf1 = new JxlImageF[2]; + + if (!SeparateMfAndHf(parameters, image0, hf0, blurTemp)) + { + return false; + } + + if (!SeparateMfAndHf(parameters, image1, hf1, blurTemp)) + { + return false; + } + + using JxlImageF blockDiffAc = new(configuration, xSize, ySize); + blockDiffAc.ZeroFill(); + + using (JxlImageF diffs = new(configuration, xSize, ySize)) + { + if (!MaltaDiffMap( + true, + image0.Plane(1), + image1.Plane(1), + WMfMalta, + WMfMalta, + Norm1Mf, + diffs, + blockDiffAc)) + { + return false; + } + + if (!MaltaDiffMap( + true, + image0.Plane(0), + image1.Plane(0), + WMfMaltaX, + WMfMaltaX, + Norm1MfX, + diffs, + blockDiffAc)) + { + return false; + } + } + + for (int c = 0; c < 3; c++) + { + L2Diff( + image0.Plane(c), + image1.Plane(c), + Wmul[3 + c], + blockDiffAc); + } + + // Free MF images + image0.Dispose(); + image1.Dispose(); + + JxlImageF[] uhf0 = new JxlImageF[2]; + JxlImageF[] uhf1 = new JxlImageF[2]; + + if (!SeparateHFAndUHF(parameters, hf0, uhf0, blurTemp)) + { + return false; + } + + if (!SeparateHFAndUHF(parameters, hf1, uhf1, blurTemp)) + { + return false; + } + + float hfAsymmetry = parameters.HfAsymmetry; + + using (JxlImageF diffs = new(configuration, xSize, ySize)) + { + MaltaDiffMap( + false, + uhf0[1], + uhf1[1], + WUhfMalta * hfAsymmetry, + WUhfMalta / hfAsymmetry, + Norm1Uhf, + diffs, + blockDiffAc); + + MaltaDiffMap( + false, + uhf0[0], + uhf1[0], + wUhfMaltaX * hfAsymmetry, + wUhfMaltaX / hfAsymmetry, + norm1UhfX, + diffs, + blockDiffAc); + + float sqrtAsym = MathF.Sqrt(hfAsymmetry); + + MaltaDiffMap( + true, + hf0[1], + hf1[1], + WHfMalta * sqrtAsym, + WHfMalta / sqrtAsym, + Norm1Hf, + diffs, + blockDiffAc); + + MaltaDiffMap( + true, + hf0[0], + hf1[0], + WHfMaltaX * sqrtAsym, + WHfMaltaX / sqrtAsym, + Norm1HfX, + diffs, + blockDiffAc); + } + + for (int c = 0; c < 2; c++) + { + L2DiffAsymmetric( + hf0[c], + hf1[c], + Wmul[c] * hfAsymmetry, + Wmul[c] / hfAsymmetry, + blockDiffAc); + } + + // Mask + using JxlImageF mask = new(configuration, xSize, ySize); + using JxlImageF mask0 = new(configuration, xSize, ySize); + using JxlImageF mask1 = new(configuration, xSize, ySize); + + CombineChannelsForMasking(hf0, uhf0, mask0); + CombineChannelsForMasking(hf1, uhf1, mask1); + + DeallocateHFAndUHF(hf0, uhf0); + DeallocateHFAndUHF(hf1, uhf1); + + if (!Mask(mask0, mask1, parameters, blurTemp, mask, blockDiffAc)) + { + return false; + } + + for (int y = 0; y < ySize; y++) + { + ReadOnlySpan dc = blockDiffDc.GetRow(y); + ReadOnlySpan ac = blockDiffAc.GetRow(y); + Span output = diffmap.GetRow(y); + ReadOnlySpan maskRow = mask.GetRow(y); + + for (int x = 0; x < xSize; x++) + { + float m = maskRow[x]; + + output[x] = + MathF.Sqrt( + (dc[x] * (float)MaskDcY(m)) + + (ac[x] * (float)MaskY(m))); + } + } + + return true; + } + + // Calculate a 2x2 subsampled image for purposes of recursive butteraugli at + // multiresolution. + public static JxlImage3F SubSample2x(Configuration configuration, JxlImage3F input) + { + int xs = (input.XSize + 1) / 2; + int ys = (input.YSize + 1) / 2; + + JxlImage3F retval = new(Configuration, xs, ys); + + for (int c = 0; c < 3; ++c) + { + for (int y = 0; y < ys; ++y) + { + for (int x = 0; x < xs; ++x) + { + retval.PlaneRow(c, y)[x] = 0.0f; + } + } + } + + for (int c = 0; c < 3; ++c) + { + for (int y = 0; y < input.YSize; ++y) + { + ReadOnlySpan srcRow = input.PlaneRow(c, y); + + for (int x = 0; x < input.XSize; ++x) + { + retval.PlaneRow(c, y / 2)[x / 2] += + 0.25f * srcRow[x]; + } + } + + if ((input.XSize & 1) != 0) + { + for (int y = 0; y < retval.YSize; ++y) + { + int lastColumn = retval.XSize - 1; + retval.PlaneRow(c, y)[lastColumn] *= 2.0f; + } + } + + if ((input.YSize & 1) != 0) + { + for (int x = 0; x < retval.XSize; ++x) + { + int lastRow = retval.YSize - 1; + retval.PlaneRow(c, lastRow)[x] *= 2.0f; + } + } + } + + return retval; + } + + public static void AddSupersampled2x(JxlImageF src, float w, JxlImageF dest) + { + const float heuristicMixingValue = 0.3f; + + for (int y = 0; y < dest.YSize; ++y) + { + Span destRow = dest.GetRow(y); + + for (int x = 0; x < dest.XSize; ++x) + { + destRow[x] *= 1.0f - (heuristicMixingValue * w); + destRow[x] += w * src.GetRow(y / 2)[x / 2]; + } + } + } +} diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs new file mode 100644 index 000000000..be5837534 --- /dev/null +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliComparator.cs @@ -0,0 +1,46 @@ +// Copyright (c) Six Labors. +// Licensed under the Six Labors Split License. + +using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; + +namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; + +internal sealed class ButteraugliComparator +{ + private int xSize; + private int ySize; + private ButteraugliParameters parameters; + private readonly ButteraugliPsychoImage pi0; + private readonly JxlImage3F temp; + private bool tempInUse; + private readonly ButteraugliBlurTemp blurTemp; + + /// + /// 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) + { + throw new NotImplementedException(); + } + + /// + /// Same as Diffmap but OpsinDynamicsImage() was already applied. + /// + public bool DiffmapOpsinDynamicsImage(JxlImage3F xyb1, JxlImageF result) + { + throw new NotImplementedException(); + } + + /// + /// Same as above but the frequency decomposition was already applied. + /// + public bool DiffmapPsychoImage(ButteraugliPsychoImage pi1, JxlImageF diffmap) + { + throw new NotImplementedException(); + } + + public bool Mask(JxlImageF mask) + { + throw new NotImplementedException(); + } +} diff --git a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliParameters.cs b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliParameters.cs index 9b5fd5761..7af303c8b 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliParameters.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/Butteraugli/ButteraugliParameters.cs @@ -9,18 +9,18 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Butteraugli; internal struct ButteraugliParameters() { /// - /// Multiplier for penalizing new HF artifacts more than + /// Gets or sets the multiplier for penalizing new HF artifacts more than /// blurring away features. Value of 1.0 represents neutral. /// - public float HfAsymmetry = 1f; + public float HfAsymmetry { get; set; } = 1f; /// - /// Multiplier for the psychovisual difference in the X channel. + /// Gets or sets the multiplier for the psychovisual difference in the X channel. /// - public float XMultiplier = 1f; + public float XMultiplier { get; set; } = 1f; /// - /// Number of nits that correspond to 1.0f input values. + /// Gets or sets the number of nits that correspond to 1.0f input values. /// - public float IntensityTarget = 80f; + public float IntensityTarget { get; set; } = 80f; } diff --git a/src/ImageSharp/Formats/Jxl/Processing/JxlQuantizer.cs b/src/ImageSharp/Formats/Jxl/Processing/JxlQuantizer.cs index 5a1405b35..423095286 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/JxlQuantizer.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/JxlQuantizer.cs @@ -4,7 +4,6 @@ using System.Buffers; using System.Runtime.CompilerServices; using SixLabors.ImageSharp.Formats.Jxl.Fields; -using SixLabors.ImageSharp.Formats.Jxl.IO; using SixLabors.ImageSharp.Formats.Jxl.Memory; using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes; diff --git a/src/ImageSharp/Formats/Jxl/Processing/JxlWeightsSeparable5.cs b/src/ImageSharp/Formats/Jxl/Processing/JxlWeightsSeparable5.cs index 249d1ec9a..7fdec43b0 100644 --- a/src/ImageSharp/Formats/Jxl/Processing/JxlWeightsSeparable5.cs +++ b/src/ImageSharp/Formats/Jxl/Processing/JxlWeightsSeparable5.cs @@ -3,9 +3,11 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing; -internal sealed class JxlWeightsSeparable5 +internal struct JxlWeightsSeparable5 { - public InlineArray12 Horizontal { get; set; } + // Don't make these a property so we can ref into them. - public InlineArray12 Vertical { get; set; } + public InlineArray12 Horizontal; + + public InlineArray12 Vertical; }