Improve astro stretch background removal

.github/copilot-instructions.md

@@ -0,0 +1,2 @@
+* Ask questions about the prompt to clarify any uncertainties.
+* Ensure that new code compiles without warnings or errors.

PhotoLocator.sln

@@ -1,7 +1,7 @@
 
 Microsoft Visual Studio Solution File, Format Version 12.00
-# Visual Studio Version 17
-VisualStudioVersion = 17.1.32210.238
+# Visual Studio Version 18
+VisualStudioVersion = 18.10.12224.181 stable
 MinimumVisualStudioVersion = 10.0.40219.1
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "PhotoLocator", "PhotoLocator\PhotoLocator.csproj", "{D9DD58A6-DE57-4DDB-92A8-2AF08B1D8309}"
 EndProject
@@ -11,6 +11,7 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution
 		.editorconfig = .editorconfig
 		.gitignore = .gitignore
 		.sonarcloud.properties = .sonarcloud.properties
+		.github\copilot-instructions.md = .github\copilot-instructions.md
 		Directory.Build.props = Directory.Build.props
 		.github\workflows\dotnet.yml = .github\workflows\dotnet.yml
 		LICENSE = LICENSE

PhotoLocator/BitmapOperations/AstroStretchOperation.cs

@@ -6,7 +6,7 @@ class AstroStretchOperation : OperationBase
     {
         public double Stretch { get; set; } = 10;
 
-        public double BackgroundSmooth { get; set; } = 8;
+        public double BackgroundSmooth { get; set; } = 2;
 
         public double BlackPoint { get; set; }
 
@@ -21,17 +21,24 @@ public override void Apply()
             }
             if (BackgroundSmooth > 0)
             {
-                var background = ConvertToGrayscaleOperation.ConvertToGrayscale(DstBitmap);
-                IIRSmoothOperation.Apply(background, (float)Math.Exp(BackgroundSmooth));
-                if (BlackPoint > 0)
-                {
-                    var bp = (float)BlackPoint;
-                    DstBitmap.ProcessElementWise(background, (p, b) => Math.Max(p - b - bp, 0));
-                }
-                else
-                    DstBitmap.ProcessElementWise(background, (p, b) => Math.Max(p - b, 0));
+                const int BackgroundWidth = 128;
+
+                var resizeOp = new LanczosResizeOperation { FilterFunc = LanczosResizeOperation.Lanczos2, FilterWindow = 1 };
+                var backgroundReduced = resizeOp.Apply(DstBitmap, BackgroundWidth, (int)Math.Ceiling((double)BackgroundWidth / DstBitmap.Width * DstBitmap.Height));
+                var grayBackground = ConvertToGrayscaleOperation.ConvertToGrayscale(backgroundReduced);
+                var mask = new FloatBitmap(grayBackground);
+                IIRSmoothOperation.Apply(grayBackground, BackgroundWidth * BackgroundSmooth);
+                mask.ProcessElementWise(grayBackground, (m, g) => m - g);
+                //DstBitmap.Assign(mask); DstBitmap.ProcessElementWise(p => p > 0 ? 1 : 0); return;
+
+                HoleClosingOperation.CloseHolesIteratively(backgroundReduced, mask, 0);
+                IIRSmoothOperation.Apply(backgroundReduced, BackgroundWidth * BackgroundSmooth);
+                var background = resizeOp.Apply(backgroundReduced, DstBitmap.Width, DstBitmap.Height);
+
+                var bp = (float)BlackPoint;
+                DstBitmap.ProcessElementWise(background, (p, b) => Math.Max(p - b - bp, 0));
             }
-            else if (BlackPoint > 0)
+            else if (BlackPoint != 0)
             {
                 var bp = (float)BlackPoint;
                 DstBitmap.ProcessElementWise(p => Math.Max(p - bp, 0));

PhotoLocator/BitmapOperations/BilinearResizeOperation.cs

@@ -162,6 +162,8 @@ public static void ApplyToPlane(BitmapPlaneInt16 srcPlane, BitmapPlaneInt16 dstP
 
         public static void ApplyToPlaneParallel(FloatBitmap srcPlane, FloatBitmap dstPlane, ROI dstROI)
         {
+            Debug.Assert(srcPlane.PlaneCount == 1);
+            Debug.Assert(dstPlane.PlaneCount == 1);
             Debug.Assert(dstROI.Right <= dstPlane.Width && dstROI.Bottom <= dstPlane.Height, "Invalid ROI");
             float scaleX = 0, scaleY = 0;
             if (dstPlane.Width > 1)

PhotoLocator/BitmapOperations/CombineFramesOperationBase.cs

@@ -157,77 +157,23 @@ protected byte[] PrepareFrame(BitmapSource image)
             return pixels;
         }
 
-        readonly struct HotPixel
-        {
-            public int Target { get; init; }
-            public int Source { get; init; }
-        }
-
-        List<HotPixel>? _hotPixels;
+        List<HoleClosingOperation.HolePatch>? _hotPixels;
 
         private void SubtractDarkFrame(byte[] pixels)
         {
             if (_darkFramePixels is null)
                 return;
             //Parallel.For(0, pixels.Length, i => pixels[i] = (byte)Math.Max(0, pixels[i] - _darkFramePixels[i]));
 
-            _hotPixels ??= FindHotPixels();
+            if (_hotPixels is null)
+            {
+                _hotPixels = HoleClosingOperation.FindHolePatches(_darkFramePixels, Width, Height, PixelSize, 32, _ct);
+                Log.Write($"Found {_hotPixels.Count} hot pixels in dark frame");
+                if (_hotPixels.Count > Width * Height * PixelSize / 2)
+                    throw new UserMessageException("Bad dark frame, too many hot pixels");
+            }
             foreach (var hotpixel in _hotPixels)
                 pixels[hotpixel.Target] = pixels[hotpixel.Source];
         }
-
-        private List<HotPixel> FindHotPixels()
-        {
-            const int HotPixelThreshold = 32;
-
-            if (_darkFramePixels is null)
-                throw new InvalidOperationException("Dark frame not set");
-            var hotPixels = new List<HotPixel>();
-            int stride = Width * PixelSize;
-            Parallel.For(0, Height, new ParallelOptions { CancellationToken = _ct }, y =>
-            {
-                var rowStart = y * stride;
-                for (int x = 0; x < stride; x++)
-                    if (_darkFramePixels[rowStart + x] > HotPixelThreshold)
-                    {
-                        bool found = false;
-                        for (int radius = 1; !found && radius < 100; radius++)
-                            for (int d = 0; !found && d < radius; d++)
-                                // ^ <1 5>
-                                // 7     ^
-                                // 3  *  4
-                                // ·     8
-                                // <6 2> ·
-                                found =
-                                    CheckCandidate(x - d * PixelSize, y - radius) ||       // 1
-                                    CheckCandidate(x + d * PixelSize, y + radius) ||       // 2
-                                    CheckCandidate(x - radius * PixelSize, y + d) ||       // 3
-                                    CheckCandidate(x + radius * PixelSize, y - d) ||       // 4
-                                    CheckCandidate(x + (d + 1) * PixelSize, y - radius) || // 5
-                                    CheckCandidate(x - (d + 1) * PixelSize, y + radius) || // 6
-                                    CheckCandidate(x - radius * PixelSize, y - (d + 1)) || // 7
-                                    CheckCandidate(x + radius * PixelSize, y + d + 1);     // 8
-                        if (!found)
-                            throw new UserMessageException("Bad dark frame, unable to patch hot pixel");
-
-                        bool CheckCandidate(int sx, int sy)
-                        {
-                            if (sx < 0 || sy < 0 || sx >= stride || sy >= Height || _darkFramePixels[sy * stride + sx] > HotPixelThreshold)
-                                return false;
-                            lock (hotPixels)
-                                hotPixels.Add(new HotPixel
-                                {
-                                    Target = rowStart + x,
-                                    Source = sy * stride + sx,
-                                });
-                            return true;
-                        }
-                    }
-            });
-            Log.Write($"Found {hotPixels.Count} hot pixels in dark frame");
-            if (hotPixels.Count > Width * Height * PixelSize / 2)
-                throw new UserMessageException("Bad dark frame, too many hot pixels");
-            return hotPixels;
-        }
     }
 }

PhotoLocator/BitmapOperations/FloatBitmap.cs

@@ -287,7 +287,7 @@ public BitmapSource ToBitmapSource16(double dpiX, double dpiY, double gamma)
                 });
             }
 
-            PixelFormat pixelFormat16 = PlaneCount switch
+            var pixelFormat16 = PlaneCount switch
             {
                 1 => PixelFormats.Gray16,
                 3 => PixelFormats.Rgb48,

PhotoLocator/BitmapOperations/HoleClosingOperation.cs

@@ -0,0 +1,139 @@
+using PhotoLocator.Helpers;
+using System;
+using System.Collections.Generic;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace PhotoLocator.BitmapOperations;
+
+static class HoleClosingOperation
+{
+    public readonly struct HolePatch
+    {
+        public int Target { get; init; }
+        public int Source { get; init; }
+    }
+
+    public static List<HolePatch> FindHolePatches(byte[] mask, int width, int height, int planeCount, int holeThreshold, CancellationToken ct)
+    {
+        var holePatches = new List<HolePatch>();
+        int stride = width * planeCount;
+        Parallel.For(0, height, new ParallelOptions { CancellationToken = ct }, y =>
+        {
+            var rowStart = y * stride;
+            for (int x = 0; x < stride; x++)
+                if (mask[rowStart + x] > holeThreshold)
+                {
+                    bool found = false;
+                    for (int radius = 1; !found && radius < 100; radius++)
+                        for (int d = 0; !found && d < radius; d++)
+                            // ^ <1 5>
+                            // 7     ^
+                            // 3  *  4
+                            // ·     8
+                            // <6 2> ·
+                            found =
+                                CheckCandidate(x - d * planeCount, y - radius) ||       // 1
+                                CheckCandidate(x + d * planeCount, y + radius) ||       // 2
+                                CheckCandidate(x - radius * planeCount, y + d) ||       // 3
+                                CheckCandidate(x + radius * planeCount, y - d) ||       // 4
+                                CheckCandidate(x + (d + 1) * planeCount, y - radius) || // 5
+                                CheckCandidate(x - (d + 1) * planeCount, y + radius) || // 6
+                                CheckCandidate(x - radius * planeCount, y - (d + 1)) || // 7
+                                CheckCandidate(x + radius * planeCount, y + d + 1);     // 8
+                    if (!found)
+                        throw new UserMessageException("Bad dark frame, unable to patch hot pixel");
+
+                    bool CheckCandidate(int sx, int sy)
+                    {
+                        if (sx < 0 || sy < 0 || sx >= stride || sy >= height || mask[sy * stride + sx] > holeThreshold)
+                            return false;
+                        lock (holePatches)
+                            holePatches.Add(new HolePatch
+                            {
+                                Target = rowStart + x,
+                                Source = sy * stride + sx,
+                            });
+                        return true;
+                    }
+                }
+        });
+        return holePatches;
+    }
+
+    /// <summary>
+    /// Fill holes in image where the corresponding pixel in the grayscale mask is above threshold
+    /// </summary>
+    /// <param name="image">Color or grayscale image</param>
+    /// <param name="mask">Grayscale mask</param>
+    /// <param name="threshold">Threshold in mask</param>
+    public static void CloseHolesIteratively(FloatBitmap image, FloatBitmap mask, float threshold)
+    {
+        if (image.Width != mask.Width || image.Height != mask.Height)
+            throw new ArgumentException("Image and mask dimensions must match.");
+
+        var width = image.Width;
+        var height = image.Height;
+        var planeCount = image.PlaneCount;
+        var holes = new bool[width * height];
+        var remaining = 0;
+
+        for (var y = 0; y < height; y++)
+            for (var x = 0; x < width; x++)
+                if (mask[x, y] > threshold)
+                {
+                    holes[y * width + x] = true;
+                    remaining++;
+                }
+        Log.Write($"Detected {remaining * 100 / holes.Length}% holes in mask");
+
+        var replacements = new List<(int X, int Y, float[] Values)>();
+        while (remaining > 0)
+        {
+            for (var y = 0; y < height; y++)
+                for (var x = 0; x < width; x++)
+                {
+                    if (!holes[y * width + x])
+                        continue;
+
+                    var sums = new float[planeCount];
+                    var count = 0;
+                    for (var dy = -1; dy <= 1; dy++)
+                        for (var dx = -1; dx <= 1; dx++)
+                        {
+                            if (dx == 0 && dy == 0)
+                                continue;
+
+                            var nx = x + dx;
+                            var ny = y + dy;
+                            if (nx < 0 || nx >= width || ny < 0 || ny >= height || holes[ny * width + nx])
+                                continue;
+
+                            for (var p = 0; p < planeCount; p++)
+                                sums[p] += image.Elements[ny, nx * planeCount + p];
+                            count++;
+                        }
+
+                    if (count > 0)
+                    {
+                        for (var p = 0; p < planeCount; p++)
+                            sums[p] /= count;
+                        replacements.Add((x, y, sums));
+                    }
+                }
+
+            if (replacements.Count == 0)
+                break;
+
+            foreach (var replacement in replacements)
+            {
+                var offset = replacement.X * planeCount;
+                for (var p = 0; p < planeCount; p++)
+                    image.Elements[replacement.Y, offset + p] = replacement.Values[p];
+                holes[replacement.Y * width + replacement.X] = false;
+            }
+            remaining -= replacements.Count;
+            replacements.Clear();
+        }
+    }
+}