Correct adaptive threshold windows for narrow regions and borders

src/ImageSharp/Processing/Processors/Binarization/AdaptiveThresholdProcessor.cs

@@ -9,8 +9,12 @@ namespace SixLabors.ImageSharp.Processing.Processors.Binarization;
 /// Performs Bradley Adaptive Threshold filter against an image.
 /// </summary>
 /// <remarks>
+/// Each pixel is compared with the mean of a centered window, clipped to the processing rectangle.
+/// The window includes the current pixel, including when the rectangle is only one pixel wide or tall.
+/// <para>
 /// Implements "Adaptive Thresholding Using the Integral Image",
 /// see paper: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.420.7883&amp;rep=rep1&amp;type=pdf
+/// </para>
 /// </remarks>
 public class AdaptiveThresholdProcessor : IImageProcessor
 {

src/ImageSharp/Processing/Processors/Binarization/AdaptiveThresholdProcessor{TPixel}.cs

@@ -38,18 +38,18 @@ protected override void OnFrameApply(ImageFrame<TPixel> source)
         TPixel lower = this.definition.Lower.ToPixel<TPixel>();
         float thresholdLimit = this.definition.ThresholdLimit;
 
-        // ClusterSize defines the size of cluster to used to check for average.
-        // Tweaked to support up to 4k wide pixels and not more. 4096 / 16 is 256 thus the '-1'
-        byte clusterSize = (byte)Math.Clamp(interest.Width / 16F, 0, 255);
+        // Use a centered window with at least one neighboring pixel on each side, where the bounds allow it.
+        // Retain the radius cap of 255 so the window does not keep growing for images wider than 4096 pixels.
+        byte windowRadius = (byte)Math.Clamp(interest.Width / 16, 1, 255);
 
         using Buffer2D<ulong> intImage = source.CalculateIntegralImage(interest);
-        RowOperation operation = new(configuration, interest, source.PixelBuffer, intImage, upper, lower, thresholdLimit, clusterSize);
-        ParallelRowIterator.IterateRows<RowOperation, L8>(
-            configuration,
-            interest,
-            in operation);
+        RowOperation operation = new(configuration, interest, source.PixelBuffer, intImage, upper, lower, thresholdLimit, windowRadius);
+        ParallelRowIterator.IterateRows<RowOperation, L8>(configuration, interest, in operation);
     }
 
+    /// <summary>
+    /// Applies the local threshold using the inclusive sums stored in the integral image.
+    /// </summary>
     private readonly struct RowOperation : IRowOperation<L8>
     {
         private readonly Configuration configuration;
@@ -61,8 +61,19 @@ protected override void OnFrameApply(ImageFrame<TPixel> source)
         private readonly float thresholdLimit;
         private readonly int startX;
         private readonly int startY;
-        private readonly byte clusterSize;
+        private readonly byte windowRadius;
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="RowOperation"/> struct.
+        /// </summary>
+        /// <param name="configuration">The configuration used for bulk luminance conversion.</param>
+        /// <param name="bounds">The image rectangle to process.</param>
+        /// <param name="source">The pixel buffer to update.</param>
+        /// <param name="intImage">The inclusive integral image, with coordinates relative to <paramref name="bounds"/>.</param>
+        /// <param name="upper">The output for samples above the local threshold.</param>
+        /// <param name="lower">The output for samples at or below the local threshold.</param>
+        /// <param name="thresholdLimit">The fraction of the local mean used as the threshold.</param>
+        /// <param name="windowRadius">The number of neighboring pixels to include on each side.</param>
         [MethodImpl(InliningOptions.ShortMethod)]
         public RowOperation(
             Configuration configuration,
@@ -72,7 +83,7 @@ public RowOperation(
             TPixel upper,
             TPixel lower,
             float thresholdLimit,
-            byte clusterSize)
+            byte windowRadius)
         {
             this.configuration = configuration;
             this.bounds = bounds;
@@ -83,7 +94,7 @@ public RowOperation(
             this.upper = upper;
             this.lower = lower;
             this.thresholdLimit = thresholdLimit;
-            this.clusterSize = clusterSize;
+            this.windowRadius = windowRadius;
         }
 
         /// <inheritdoc/>
@@ -98,22 +109,37 @@ public void Invoke(int y, Span<L8> span)
             Span<TPixel> rowSpan = this.source.DangerousGetRowSpan(y).Slice(this.startX, span.Length);
             PixelOperations<TPixel>.Instance.ToL8(this.configuration, rowSpan, span);
 
-            int startY = this.startY;
             int maxX = this.bounds.Width - 1;
-            int maxY = this.bounds.Height - 1;
-            int clusterSize = this.clusterSize;
+            int windowRadius = this.windowRadius;
             float thresholdLimit = this.thresholdLimit;
-            Buffer2D<ulong> image = this.intImage;
+
+            // Integral coordinates are relative to the processing rectangle. Include the current pixel
+            // and clip each side of its centered window to that rectangle, including one-pixel dimensions.
+            int relativeY = y - this.startY;
+            int top = Math.Max(relativeY - windowRadius, 0);
+            int bottom = Math.Min(relativeY + windowRadius, this.bounds.Height - 1);
+            uint windowHeight = (uint)(bottom - top + 1);
+
+            // I(x,y) includes row y and column x. Subtract the prefix ending just before the window.
+            // An empty prefix span represents the image's top boundary. Acquire both spans once per row.
+            ReadOnlySpan<ulong> bottomRow = this.intImage.DangerousGetRowSpan(bottom);
+            ReadOnlySpan<ulong> precedingRow = top > 0 ? this.intImage.DangerousGetRowSpan(top - 1) : default;
+
             for (int x = 0; x < rowSpan.Length; x++)
             {
-                int x1 = Math.Clamp(x - clusterSize + 1, 0, maxX);
-                int x2 = Math.Min(x + clusterSize + 1, maxX);
-                int y1 = Math.Clamp(y - startY - clusterSize + 1, 0, maxY);
-                int y2 = Math.Min(y - startY + clusterSize + 1, maxY);
+                int left = Math.Max(x - windowRadius, 0);
+                int right = Math.Min(x + windowRadius, maxX);
+                uint count = (uint)(right - left + 1) * windowHeight;
 
-                uint count = (uint)((x2 - x1) * (y2 - y1));
-                ulong sum = Math.Min(image[x2, y2] - image[x1, y2] - image[x2, y1] + image[x1, y1], ulong.MaxValue);
+                ulong sum = bottomRow[right] - (precedingRow.IsEmpty ? 0 : precedingRow[right]);
+                if (left > 0)
+                {
+                    // Remove the columns before the window from the same vertical band. Subtracting
+                    // these nonnegative band sums avoids unsigned underflow between integral corners.
+                    sum -= bottomRow[left - 1] - (precedingRow.IsEmpty ? 0 : precedingRow[left - 1]);
+                }
 
+                // The inclusive window always contains the current pixel, so count cannot be zero.
                 if (span[x].PackedValue * count <= sum * thresholdLimit)
                 {
                     rowSpan[x] = this.lower;

tests/ImageSharp.Tests/Processing/Binarization/AdaptiveThresholdTests.cs

@@ -113,6 +113,76 @@ public void AdaptiveThreshold_Works<TPixel>(TestImageProvider<TPixel> provider)
         }
     }
 
+    /// <summary>
+    /// Checks every pixel, including the last row and column, in a narrow image.
+    /// </summary>
+    /// <param name="width">The image width.</param>
+    /// <param name="height">The image height.</param>
+    [Theory]
+    [InlineData(3, 3)]
+    [InlineData(8, 8)]
+    [InlineData(15, 15)]
+    [InlineData(15, 3)]
+    public void AdaptiveThreshold_NarrowImage_SeparatesDarkPixel(int width, int height)
+    {
+        using Image<Rgba32> image = new(width, height, new Rgba32(200, 200, 200));
+        Point dark = new(width / 2, height / 2);
+        image[dark.X, dark.Y] = new Rgba32(10, 10, 10);
+
+        image.Mutate(img => img.AdaptiveThreshold());
+
+        AssertOnlyDarkPixelIsLower(image, new Rectangle(0, 0, width, height), dark);
+    }
+
+    /// <summary>
+    /// Checks a narrow processing rectangle that does not start at the image origin.
+    /// </summary>
+    [Fact]
+    public void AdaptiveThreshold_NarrowRectangle_SeparatesDarkPixel()
+    {
+        using Image<Rgba32> image = new(100, 100, new Rgba32(200, 200, 200));
+        Rectangle rectangle = new(10, 10, 8, 50);
+        Point dark = new(14, 30);
+        image[dark.X, dark.Y] = new Rgba32(10, 10, 10);
+
+        image.Mutate(img => img.AdaptiveThreshold(Color.White, Color.Black, rectangle));
+
+        AssertOnlyDarkPixelIsLower(image, rectangle, dark);
+    }
+
+    /// <summary>
+    /// Checks centered windows against direct sums, including single-pixel dimensions and rectangle borders.
+    /// </summary>
+    /// <param name="width">The processing rectangle width.</param>
+    /// <param name="height">The processing rectangle height.</param>
+    /// <param name="radius">The expected window radius.</param>
+    [Theory]
+    [InlineData(1, 1, 1)]
+    [InlineData(1, 8, 1)]
+    [InlineData(8, 1, 1)]
+    [InlineData(15, 7, 1)]
+    [InlineData(16, 9, 1)]
+    [InlineData(31, 11, 1)]
+    [InlineData(32, 17, 2)]
+    public void AdaptiveThreshold_Rectangle_MatchesLocalMean(int width, int height, int radius)
+    {
+        Rectangle rectangle = new(3, 4, width, height);
+        using Image<L8> source = new(width + 6, height + 8);
+        for (int y = 0; y < source.Height; y++)
+        {
+            for (int x = 0; x < source.Width; x++)
+            {
+                // An asymmetric pattern detects windows shifted away from the current pixel.
+                source[x, y] = new L8((byte)((17 + (37 * x) + (61 * y)) % 256));
+            }
+        }
+
+        const float thresholdLimit = .85F;
+        using Image<L8> actual = source.Clone(context => context.AdaptiveThreshold(Color.White, Color.Black, thresholdLimit, rectangle));
+
+        AssertMatchesLocalMean(source, actual, rectangle, radius, thresholdLimit);
+    }
+
     [Theory]
     [WithFile(TestImages.Png.Bradley02, PixelTypes.Rgba32)]
     public void AdaptiveThreshold_WithRectangle_Works<TPixel>(TestImageProvider<TPixel> provider)
@@ -125,4 +195,62 @@ public void AdaptiveThreshold_WithRectangle_Works<TPixel>(TestImageProvider<TPix
             image.CompareToReferenceOutput(ImageComparer.Exact, provider);
         }
     }
+
+    /// <summary>
+    /// Asserts that a thresholded region of uniform light pixels with a single dark pixel
+    /// maps the dark pixel to black and every other pixel to white.
+    /// </summary>
+    /// <param name="image">The thresholded image.</param>
+    /// <param name="region">The region that was thresholded.</param>
+    /// <param name="dark">The location of the dark pixel.</param>
+    private static void AssertOnlyDarkPixelIsLower(Image<Rgba32> image, Rectangle region, Point dark)
+    {
+        Rgba32 black = Color.Black.ToPixel<Rgba32>();
+        Rgba32 white = Color.White.ToPixel<Rgba32>();
+
+        for (int y = region.Top; y < region.Bottom; y++)
+        {
+            for (int x = region.Left; x < region.Right; x++)
+            {
+                Assert.Equal(x == dark.X && y == dark.Y ? black : white, image[x, y]);
+            }
+        }
+    }
+
+    /// <summary>
+    /// Checks the local mean without using an integral image, and checks that pixels outside the rectangle are unchanged.
+    /// </summary>
+    /// <param name="source">The unchanged grayscale input.</param>
+    /// <param name="actual">The processed image.</param>
+    /// <param name="rectangle">The processing rectangle.</param>
+    /// <param name="radius">The expected window radius.</param>
+    /// <param name="thresholdLimit">The fraction of the local mean used as the threshold.</param>
+    private static void AssertMatchesLocalMean(Image<L8> source, Image<L8> actual, Rectangle rectangle, int radius, float thresholdLimit)
+    {
+        for (int y = 0; y < source.Height; y++)
+        {
+            for (int x = 0; x < source.Width; x++)
+            {
+                L8 expected = source[x, y];
+                if (rectangle.Contains(x, y))
+                {
+                    int sum = 0;
+                    int count = 0;
+                    for (int sampleY = Math.Max(y - radius, rectangle.Top); sampleY <= Math.Min(y + radius, rectangle.Bottom - 1); sampleY++)
+                    {
+                        for (int sampleX = Math.Max(x - radius, rectangle.Left); sampleX <= Math.Min(x + radius, rectangle.Right - 1); sampleX++)
+                        {
+                            sum += source[sampleX, sampleY].PackedValue;
+                            count++;
+                        }
+                    }
+
+                    // Include both window endpoints. A one-pixel window therefore has count = 1.
+                    expected = new L8(source[x, y].PackedValue * count <= sum * thresholdLimit ? (byte)0 : (byte)255);
+                }
+
+                Assert.Equal(expected, actual[x, y]);
+            }
+        }
+    }
 }