FindResinTraps.cs

using Emgu.CV;
using Emgu.CV.CvEnum;
using Emgu.CV.Structure;
using Emgu.CV.Util;
using System;
using System.Collections;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Diagnostics;
using System.Drawing;
using System.IO;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Windows.Forms;

namespace ResinTrapPlayground
{
    public partial class FindResinTraps : Form
    {
        /* The following is a novel implementation of resin/suction trap detection for MSLA sliced images *
         * It largely consists of 2 passes through the layers. 
         * 
         * The 1st pass processes the layers from the bottom to the top, keeping track of what areas in 
         * the previous layer are air (or connected to air in some way). The initial air map is derived from
         * the starting layer, by inverting the layer and then filling all external contours black (not air)
         * which results in an inital air map of everything that is not solid, or inside a solid, being considered
         * air. While traversing the layers, the air map is updated to include any new air regions, or hollow areas
         * that intersect with the air map (since while they are not air, they are connected to air). Any regions 
         * that do not overlap the air map are declared resin traps.
         * 
         * The 2nd pass processes the layers from the top to the bottom, the initial air map differs from
         * the 1st pass. For the 2nd pass, anything that is not solid is considered to be air (a bitwise 
         * not of the top layer is used). The airmap is updated per layer just like it was in the 1st pass,
         * however we do not test all contours of the layer, only those that were marked as resin traps by 
         * the 1st pass. This is acceptable because in the event we did test all hollow areas, and one of those
         * was *not* detected on 1st pass, it must be air-accessible from the bottom. In the even that a contour
         * we are testing does overlap the air map during the 2nd pass, it is cleared from the "resin trap" list
         * and instead added to the "suction trap" list. This is to indicate that while there is access to air from
         * the top, there are regions that could be affected by vacuum/suction forces during the print
         */

        public FindResinTraps()
        {
            //string[] files = Directory.GetFiles(@"C:\Users\tslat\Downloads\RERF\R_E_R_F", "layer*.png");

            /** Just a bunch of code to load layer images from disk (extracted by UVTools) and store as PNG data **/
            /** This is done this way so the PNG data in memory is of a grayscale image, not full color in case that matters **/
            string[] files = Directory.GetFiles(@"C:\Users\tslat\Downloads\ResinTrapTests\hollow_yoda_no_print", "layer*.png");
            int layerCount = files.Length;

            List<byte[]> layerData = new List<byte[]>();
            var loadedMats = new Mat[files.Length];
            string windowName = "Debug"; //The name of the window
            CvInvoke.NamedWindow(windowName);
            for (var x = 0; x < loadedMats.Length; x++)
            {
                loadedMats[x] = new Mat(files[x], ImreadModes.Grayscale);
                CvInvoke.Threshold(loadedMats[x], loadedMats[x], 1, byte.MaxValue, ThresholdType.Binary);
                layerData.Add(CvInvoke.Imencode(".png", loadedMats[x]));
                loadedMats[x].Dispose();
            }

            Stopwatch sw = new Stopwatch();
            sw.Start();
            var (traps, suctions) = FindResinAndSuctionTraps(layerData, layerCount);
            sw.Stop();

            int layersWithTraps = traps.Where(t => t != null && t.Count > 0).Count();
            int layersWithSuction = suctions.Where(t => t != null && t.Count > 0).Count();

            InitializeComponent();
        }

        public (List<VectorOfVectorOfPoint>[], List<VectorOfVectorOfPoint>[]) FindResinAndSuctionTraps(List<byte[]> layerData, int layerCount)
        {
            /* these hold an array (one element for each layer), that has a List of contours, or null */
            List<VectorOfVectorOfPoint>[] resinTraps = new List<VectorOfVectorOfPoint>[layerCount];
            List<VectorOfVectorOfPoint>[] suctionTraps = new List<VectorOfVectorOfPoint>[layerCount];

            /* this maintains a Mat that reflects the current air-accessible region for the layer that is being processed, black is no air, white is air (or air accessible) */
            Mat currentAirMap = null;

            /* the first pass does bottom to top, and tracks anything it thinks is a resin trap */
            for (var layerIndex = 0; layerIndex < layerData.Count; layerIndex++)
            {
                using Mat curLayer = new();
                CvInvoke.Imdecode(layerData[layerIndex], ImreadModes.Grayscale, curLayer);

                var layerAirMap = GenerateAirMap(curLayer);
                if (layerIndex == 0)
                {
                    currentAirMap = layerAirMap.Clone();
                }                

                /* remove solid areas of current layer from the air map */
                CvInvoke.Subtract(currentAirMap, curLayer, currentAirMap);

                /* add in areas of air in current layer to air map */
                CvInvoke.BitwiseOr(layerAirMap, currentAirMap, currentAirMap);

                /* all done with this mat, clean it up */
                layerAirMap.Dispose();

                /* find hollows of current layer */
                using VectorOfVectorOfPoint contours = new();
                using Mat hierarchy = new();
                CvInvoke.FindContours(curLayer, contours, hierarchy, RetrType.Tree, ChainApproxMethod.ChainApproxSimple);
                var arr = hierarchy.GetData();

                /* hollow contour grouping */
                var hollowGroups = new List<VectorOfVectorOfPoint>();
                var processedContours = new bool[contours.Size];
                for (int i = 0; i < contours.Size; i++)
                {
                    if (processedContours[i]) continue;
                    processedContours[i] = true;
                    if ((int)arr.GetValue(0, i, 3) == -1) continue;
                    hollowGroups.Add(new VectorOfVectorOfPoint(contours[i]));
                    for (int n = i + 1; n < contours.Size; n++)
                    {
                        if (processedContours[n] || (int)arr.GetValue(0, n, 3) != i) continue;
                        processedContours[n] = true;
                        hollowGroups[^1].Push(contours[n]);
                    }

                }

                for (var i = 0; i < hollowGroups.Count; i++)
                {
                    /* intersect current contour, with the current airmap. */
                    Mat currentContour = Mat.Zeros(curLayer.Height, curLayer.Width, curLayer.Depth, curLayer.NumberOfChannels);
                    CvInvoke.DrawContours(currentContour, hollowGroups[i], -1, new MCvScalar(255, 255, 255, 255), -1);
                    Mat airOverlap = currentAirMap.Clone();
                    CvInvoke.BitwiseAnd(currentAirMap, currentContour, airOverlap);
                    int overlapCount = CvInvoke.CountNonZero(airOverlap);
                    if (overlapCount == 0)
                    {
                        /* this countour does *not* overlap known air */
                        if (resinTraps[layerIndex] == null)
                        {
                            resinTraps[layerIndex] = new List<VectorOfVectorOfPoint>();
                        }

                        /* add a resin trap (for now... will be revisited in part 2) */
                        resinTraps[layerIndex].Add(hollowGroups[i]);
                    }
                    else
                    {
                        /* this contour does overlap air, add it to the current air map */
                        CvInvoke.BitwiseOr(currentContour, currentAirMap, currentAirMap);
                    }

                    /* cleanup */
                    currentContour.Dispose();
                    airOverlap.Dispose();
                };
            }

            /* starting over again but this time from the top to the bottom */
            currentAirMap = null;
            
            for (var layerIndex = resinTraps.Length - 1; layerIndex >= 0; layerIndex--)
            {
                using Mat curLayer = new();
                CvInvoke.Imdecode(layerData[layerIndex], ImreadModes.Grayscale, curLayer);

                if (layerIndex == resinTraps.Length - 1)
                {
                    /* this is subtly different that for the first pass, we don't use GenerateAirMap for the initial airmap */
                    /* instead we use a bitwise not, this way anything that is open/hollow on the top layer is treated as air */
                    currentAirMap = curLayer.Clone();
                    CvInvoke.BitwiseNot(curLayer, currentAirMap);
                }

                /* we still modify the airmap like normal, where we account for the air areas of the layer, and any contours that might overlap...*/
                var layerAirMap = GenerateAirMap(curLayer);

                /* remove solid areas of current layer from the air map */
                CvInvoke.Subtract(currentAirMap, curLayer, currentAirMap);

                /* add in areas of air in current layer to air map */
                CvInvoke.BitwiseOr(layerAirMap, currentAirMap, currentAirMap);
                layerAirMap.Dispose();

                if (resinTraps[layerIndex] != null)
                {
                    /* here we don't worry about finding contours on the layer, the bottom to top pass did that already */
                    /* all we care about is contours the first pass thought were resin traps, since there was no access to air from the bottom */
                    for (var x = 0; x < resinTraps[layerIndex].Count; x++)
                    {
                        /* check if each contour overlaps known air */
                        Mat currentContour = Mat.Zeros(curLayer.Height, curLayer.Width, curLayer.Depth, curLayer.NumberOfChannels);
                        CvInvoke.DrawContours(currentContour, resinTraps[layerIndex][x], -1, new MCvScalar(255, 255, 255, 255), -1);

                        Mat airOverlap = currentAirMap.Clone();
                        CvInvoke.BitwiseAnd(currentAirMap, currentContour, airOverlap);
                        int overlapCount = CvInvoke.CountNonZero(airOverlap);

                        /* if this trap still doesn't overlap with air from the top, its an actual resin trap - leave it as is*/

                        if (overlapCount > 0) { 
                            /* this contour does overlap air, add this it our air map */
                            CvInvoke.BitwiseOr(currentContour, currentAirMap, currentAirMap);

                            /* if we haven't defined a suctionTrap list for this layer, do so */
                            if (suctionTraps[layerIndex] == null)
                            {
                                suctionTraps[layerIndex] = new List<VectorOfVectorOfPoint>();
                            }

                            /* since we know it isn't a resin trap, it becomes a suction trap */
                            suctionTraps[layerIndex].Add(resinTraps[layerIndex][x]);

                            /* to keep things tidy while we iterate resin traps, it will be left in the list for now, and removed later */
                        }
                        currentContour.Dispose();
                        airOverlap.Dispose();
                    };

                    /* anything that converted to a suction trap needs to removed from resinTraps. Loop backwards so indexes don't shift */
                    if (suctionTraps[layerIndex] != null)
                    {
                        for (var x = suctionTraps[layerIndex].Count - 1; x >= 0; x--)
                        {
                            resinTraps[layerIndex].Remove(suctionTraps[layerIndex][x]);
                            if (resinTraps[layerIndex].Count == 0)
                            {
                                resinTraps[layerIndex] = null;
                            }
                        }
                    }
                }
            }

            /* some rendering code that draws resin traps blue, and suction traps green */
            for (var layerIndex = 0; layerIndex < resinTraps.Length; layerIndex++)
            {
                using Mat curLayer = new();
                CvInvoke.Imdecode(layerData[layerIndex], ImreadModes.Color, curLayer);

                if (resinTraps[layerIndex] != null && resinTraps[layerIndex].Count > 0)
                {
                    foreach (var r in resinTraps[layerIndex])
                    {
                        CvInvoke.DrawContours(curLayer, r, -1, new MCvScalar(255, 0, 0, 255), -1);
                    }
                }

                if (suctionTraps[layerIndex] != null)
                {
                    foreach (var s in suctionTraps[layerIndex])
                    {
                        CvInvoke.DrawContours(curLayer, s, -1, new MCvScalar(0, 255, 0, 255), -1);
                    }
                }

                ShowImage(curLayer);
            }

            CvInvoke.DestroyAllWindows();

            /* return the arrays we found */
            return (resinTraps,suctionTraps);
        }

        public void ShowImage(Mat img)
        {
            Mat resized = Mat.Zeros(img.Height, img.Width, img.Depth, img.NumberOfChannels);
            //CvInvoke.Resize(img, resized, new Size(1024, 576));
            CvInvoke.Resize(img, resized, new Size((int)(1024 * ((double)img.Width/img.Height)), 1024));
            CvInvoke.Imshow("Debug", resized);
            resized.Dispose();
            CvInvoke.WaitKey(40);
        }


        public Mat GenerateAirMap(Mat input)
        {
            var clone = input.Clone();
            CvInvoke.BitwiseNot(input, clone);
            using VectorOfVectorOfPoint contours = new();
            CvInvoke.FindContours(input, contours, null, RetrType.External, ChainApproxMethod.ChainApproxSimple);
            CvInvoke.DrawContours(clone, contours, -1, new MCvScalar(0,0,0, 255), -1);
            return clone;
        }
    }
}