Fold in SIMD optimized PHash implementation

Remove separate `Optimized` file/tests

Author: jaden-young

src/ImageHash/HashAlgorithms/PerceptualHash.cs

@@ -4,8 +4,9 @@
     using System.Collections.Generic;
     using System.Diagnostics;
     using System.Linq;
+    using System.Numerics;
     using System.Runtime.CompilerServices;
-
+    using System.Threading.Tasks;
     using SixLabors.ImageSharp;
     using SixLabors.ImageSharp.PixelFormats;
     using SixLabors.ImageSharp.Processing;
@@ -18,16 +19,18 @@ public class PerceptualHash : IImageHash
         private const int Size = 64;
         private static readonly double Sqrt2DivSize = Math.Sqrt(2D / Size);
         private static readonly double Sqrt2 = 1 / Math.Sqrt(2);
+        private static readonly double[,] _dctCoeffs = GenerateDctCoeffs();
+        private static readonly List<Vector<double>>[] _dctCoeffsSimd = GenerateDctCoeffsSimd();
 
         /// <inheritdoc />
         public ulong Hash(Image<Rgba32> image)
         {
             if (image == null)
                 throw new ArgumentNullException(nameof(image));
 
-            var rows = new double[Size][];
+            var rows = new double[Size, Size];
             var sequence = new double[Size];
-            var matrix = new double[Size][];
+            var matrix = new double[Size, Size];
 
             image.Mutate(ctx => ctx
                                 .Resize(Size, Size)
@@ -40,38 +43,36 @@ public ulong Hash(Image<Rgba32> image)
                 for (var x = 0; x < Size; x++)
                     sequence[x] = image[x, y].R;
 
-                rows[y] = Dct1D(sequence);
+                Dct1D_SIMD(sequence, rows, y);
             }
 
             // Calculate the DCT for each column.
-            for (var x = 0; x < Size; x++)
+            for (var x = 0; x < 8; x++)
             {
                 for (var y = 0; y < Size; y++)
-                    sequence[y] = rows[y][x];
+                    sequence[y] = rows[y, x];
 
-                matrix[x] = Dct1D(sequence);
+                Dct1D_SIMD(sequence, matrix, x, limit: 8);
             }
 
             // Only use the top 8x8 values.
-            var top8X8 = new List<double>(Size);
+            var top8X8 = new double[Size];
             for (var y = 0; y < 8; y++)
             {
                 for (var x = 0; x < 8; x++)
-                    top8X8.Add(matrix[y][x]);
+                    top8X8[(y * 8) + x] = matrix[y, x];
             }
 
-            var topRight = top8X8.ToArray();
-
             // Get Median.
-            var median = CalculateMedian64Values(topRight);
+            var median = CalculateMedian64Values(top8X8);
 
             // Calculate hash.
             var mask = 1UL << (Size - 1);
             var hash = 0UL;
 
             for (var i = 0; i < Size; i++)
             {
-                if (topRight[i] > median)
+                if (top8X8[i] > median)
                     hash |= mask;
 
                 mask = mask >> 1;
@@ -87,28 +88,74 @@ private static double CalculateMedian64Values(IReadOnlyCollection<double> values
             return values.OrderBy(value => value).Skip(31).Take(2).Average();
         }
 
+        private static double[,] GenerateDctCoeffs()
+        {
+            double[,] c = new double[Size, Size];
+            for (var coef = 0; coef < Size; coef++)
+            {
+                for (var i = 0; i < Size; i++)
+                {
+                    c[i, coef] = Math.Cos(((2.0 * i) + 1.0) * coef * Math.PI / (2.0 * Size));
+                }
+            }
+
+            return c;
+        }
+
+        private static List<Vector<double>>[] GenerateDctCoeffsSimd()
+        {
+            List<Vector<double>>[] results = new List<Vector<double>>[Size];
+            for (var coef = 0; coef < Size; coef++)
+            {
+                var singleResultRaw = new double[Size];
+                for (var i = 0; i < Size; i++)
+                {
+                    singleResultRaw[i] = Math.Cos(((2.0 * i) + 1.0) * coef * Math.PI / (2.0 * Size));
+                }
+
+                var singleResultList = new List<Vector<double>>();
+                var stride = Vector<double>.Count;
+                Debug.Assert(Size % stride == 0, "Size must be a multiple of SIMD stride");
+                for (int i = 0; i < Size; i += stride)
+                {
+                    var v = new Vector<double>(singleResultRaw, i);
+                    singleResultList.Add(v);
+                }
+
+                results[coef] = singleResultList;
+            }
+
+            return results;
+        }
+
         /// <summary>
         /// One dimensional Discrete Cosine Transformation.
         /// </summary>
-        /// <param name="values">Should be an array of doubles of length 64.</param>
-        /// <returns>array of doubles of length 64.</returns>
+        /// <param name="valuesRaw">Should be an array of doubles of length 64.</param>
+        /// <param name="coefficients">Coefficients.</param>
+        /// <param name="ci">Coefficients index.</param>
+        /// <param name="limit">Limit.</param>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static double[] Dct1D(IReadOnlyList<double> values)
+        private static void Dct1D_SIMD(double[] valuesRaw, double[,] coefficients, int ci, int limit = Size)
         {
-            Debug.Assert(values.Count == 64, "This DCT method works with 64 doubles.");
-            var coefficients = new double[Size];
+            Debug.Assert(valuesRaw.Length == 64, "This DCT method works with 64 doubles.");
 
-            for (var coef = 0; coef < Size; coef++)
+            var valuesList = new List<Vector<double>>();
+            var stride = Vector<double>.Count;
+            for (int i = 0; i < valuesRaw.Length; i += stride)
             {
-                for (var i = 0; i < Size; i++)
-                    coefficients[coef] += values[i] * Math.Cos(((2.0 * i) + 1.0) * coef * Math.PI / (2.0 * Size));
+                valuesList.Add(new Vector<double>(valuesRaw, i));
+            }
 
-                coefficients[coef] *= Sqrt2DivSize;
+            for (var coef = 0; coef < limit; coef++)
+            {
+                for (int i = 0; i < valuesList.Count; i++)
+                    coefficients[ci, coef] += Vector.Dot(valuesList[i], _dctCoeffsSimd[coef][i]);
+
+                coefficients[ci, coef] *= Sqrt2DivSize;
                 if (coef == 0)
-                    coefficients[coef] *= Sqrt2;
+                    coefficients[ci, coef] *= Sqrt2;
             }
-
-            return coefficients;
         }
     }
 }