Zhang-Suen thinning algorithm C#

Question

I am trying to write a Zhang-Suen thinning algorithm in C# following this guideline, without processing the margins.

In the function 'zhangsuen', I am reading from the image 'imgUndo' and writting to the image 'img'. The pointers dataPtrOrigin_aux inside the for cycles are used to read the 9 pixels inside a 3x3 window such that dataPtrOrigin_aux5 is the central pixel of this window, and that window will move along the whole image, moving from left to right and from top to bottom. In each pixel, if the if the statements are verified to be true, the corresponding changes are made in the image to be written by the pointer dataPtrFinal.

Note that I stored the neighbours of the current pixel inside a 8 element array. As such, they will be stored following this order:

        internal static void zhangsuen(Image<Bgr, byte> img, Image<Bgr, byte> imgUndo)
    {
        unsafe
        {

            MIplImage m = img.MIplImage; //Image to be written.
            MIplImage mUndo = imgUndo.MIplImage; //Image to be read.
            byte* dataPtrFinal = (byte*)m.imageData.ToPointer();
            byte* dataPtrUndo = (byte*)mUndo.imageData.ToPointer();

            int width = img.Width; //Width of the image.
            int height = img.Height; //Height of the image.
            int nChan = m.nChannels; //3 channels (R, G, B).
            int wStep = m.widthStep; //Total width of the image (including padding).
            int padding = wStep - nChan * width; //Padding at the end of each line.

            int x, y, i;

            int[] neighbours = new int[8]; //Store the value of the surrounding neighbours in this array.

            int step; //Step 1 or 2.
            int[] sequence = { 1, 2, 4, 7, 6, 5, 3, 0, 1 };
            int blackn = 0; //Number of black neighbours.
            int numtransitions = 0; //Number of transitions from white to black in the sequence specified by the array sequence.
            int changed = 1; //Just so it enters the while.

            bool isblack = false;

            int counter = 0;


            while(changed > 0)
            {
                changed = 0;

                if (counter % 2 == 0) //We want to read all the pixels in the image before going to the next step
                    step = 1;
                else
                    step = 2;

                for (y = 0; y < height; y++)
                {
                    for (x = 0; x < width; x++)
                    {

                            byte* dataPtrOrigin_aux1 = (byte*)(dataPtrUndo + (y - 1) * m.widthStep + (x - 1) * m.nChannels);
                            byte* dataPtrOrigin_aux2 = (byte*)(dataPtrUndo + (y - 1) * m.widthStep + (x) * m.nChannels);
                            byte* dataPtrOrigin_aux3 = (byte*)(dataPtrUndo + (y - 1) * m.widthStep + (x + 1) * m.nChannels);
                            byte* dataPtrOrigin_aux4 = (byte*)(dataPtrUndo + (y) * m.widthStep + (x - 1) * m.nChannels);
                            byte* dataPtrOrigin_aux5 = (byte*)(dataPtrUndo + (y) * m.widthStep + (x) * m.nChannels);
                            byte* dataPtrOrigin_aux6 = (byte*)(dataPtrUndo + (y) * m.widthStep + (x + 1) * m.nChannels);
                            byte* dataPtrOrigin_aux7 = (byte*)(dataPtrUndo + (y + 1) * m.widthStep + (x - 1) * m.nChannels);
                            byte* dataPtrOrigin_aux8 = (byte*)(dataPtrUndo + (y + 1) * m.widthStep + (x) * m.nChannels);
                            byte* dataPtrOrigin_aux9 = (byte*)(dataPtrUndo + (y + 1) * m.widthStep + (x + 1) * m.nChannels);


                        if (x > 0 && y > 0 && x < width - 1 && y < height - 1)
                        {
                            if (dataPtrOrigin_aux5[0] == 0)
                                isblack = true;

                            if (isblack)
                            {

                                neighbours[0] = dataPtrOrigin_aux1[0];
                                neighbours[1] = dataPtrOrigin_aux2[0];
                                neighbours[2] = dataPtrOrigin_aux3[0];
                                neighbours[3] = dataPtrOrigin_aux4[0];
                                neighbours[4] = dataPtrOrigin_aux6[0];
                                neighbours[5] = dataPtrOrigin_aux7[0];
                                neighbours[6] = dataPtrOrigin_aux8[0];
                                neighbours[7] = dataPtrOrigin_aux9[0];

                                for(i = 0; i <= 7; i++)
                                {
                                    if (neighbours[i] == 0)
                                        blackn++;

                                    if (neighbours[sequence[i]] - neighbours[sequence[i + 1]] == 255) //número de transições de branco para preto, seguindo a ordem do vector sequence
                                        numtransitions++;
                                }


                                if ((blackn >= 2 && blackn <= 6) && numtransitions == 1)
                                {
                                        if (step == 1 && (neighbours[1] == 255 || neighbours[4] == 255 || neighbours[6] == 255) && (neighbours[4] == 255 || neighbours[6] == 255 || neighbours[3] == 255))
                                        {
                                            dataPtrFinal[0] = 255;
                                            dataPtrFinal[1] = 255;
                                            dataPtrFinal[2] = 255;

                                            changed++;
                                        }

                                        if (step == 2 && (neighbours[1] == 255 || neighbours[4] == 255 || neighbours[3] == 255) && (neighbours[1] == 255 || neighbours[6] == 255 || neighbours[3] == 255))
                                        {
                                            dataPtrFinal[0] = 255;
                                            dataPtrFinal[1] = 255;
                                            dataPtrFinal[2] = 255;

                                            changed++;
                                        }

                                }
                            }
                        }


                        dataPtrFinal += nChan;

                        isblack = false;
                        blackn = 0;
                        numtransitions = 0;

                    }
                    dataPtrFinal += padding;

                }


                dataPtrUndo = (byte*)m.imageData.ToPointer(); //Change the image to be read to the one that has just been written.

                counter++;

            }


        }
    }

As I end reading the first image and writing the changes to the image 'img' (As soon as the cycle for (y = 0; y < height; y++) ends I want the image I have just written to be the one I will read in the next cycle so that further thinning is made. I tried to accomplish this with the line

dataPtrUndo = (byte*)m.imageData.ToPointer();

Although at some value of counter that is greater than 0 (depends on the image that is read) I get an error that says that protected memory has been tried to be written which indicates I have tried to write outside of the image limits, but I don't understand why. Is it the last attribution to dataPtrUndo that I am doing erroneously?

Answer 1

bwang22's answer works. Sort of. But with two issues: It doesn't do the iterations until no more changes happen. And it does a shallow copy of the Array.. The two issues cooperate so to speak, cancelling each other out, resulting in an thinning, but not the best looking one.

Here is the corrected code, which gives a nicer looking result:

First two methods to convert from Image to bool[][] and back; the functions are not optimzed for speed; if you need that go for lockbits/unsafe..:

public static bool[][] Image2Bool(Image img)
{
    Bitmap bmp = new Bitmap(img);
    bool[][] s = new bool[bmp.Height][];
    for (int y = 0; y < bmp.Height; y++ )
    {
        s[y] = new bool[bmp.Width];
        for (int x = 0; x < bmp.Width; x++)
            s[y][x] = bmp.GetPixel(x, y).GetBrightness() < 0.3;
    }
    return s;

}

public static Image Bool2Image(bool[][] s)
{
    Bitmap bmp = new Bitmap(s[0].Length, s.Length);
    using (Graphics g = Graphics.FromImage(bmp)) g.Clear(Color.White);
    for (int y = 0; y < bmp.Height; y++)
        for (int x = 0; x < bmp.Width; x++)
            if (s[y][x]) bmp.SetPixel(x, y, Color.Black);

    return (Bitmap)bmp;
}

Now the corrected thinning code, much of it more or less unchanged from bwang22's answer:

public static bool[][] ZhangSuenThinning(bool[][] s)
{
    bool[][] temp = ArrayClone(s);  // make a deep copy to start.. 
    int count = 0;
    do  // the missing iteration
    {
        count  = step(1, temp, s);
        temp = ArrayClone(s);      // ..and on each..
        count += step(2, temp, s);
        temp = ArrayClone(s);      // ..call!
    }
    while (count > 0);

    return s;
}

static int step(int stepNo, bool[][] temp, bool[][] s)
{
    int count = 0;

    for (int a = 1; a < temp.Length - 1; a++)
    {
        for (int b = 1; b < temp[0].Length - 1; b++)
        {
            if (SuenThinningAlg(a, b, temp, stepNo == 2))
            {
                // still changes happening?
                if (s[a][b]) count++; 
                s[a][b] = false; 
            }
        }
    }
    return count;
}

static bool SuenThinningAlg(int x, int y, bool[][] s, bool even)
{
    bool p2 = s[x][y - 1];
    bool p3 = s[x + 1][y - 1];
    bool p4 = s[x + 1][y];
    bool p5 = s[x + 1][y + 1];
    bool p6 = s[x][y + 1];
    bool p7 = s[x - 1][y + 1];
    bool p8 = s[x - 1][y];
    bool p9 = s[x - 1][y - 1];


    int bp1 = NumberOfNonZeroNeighbors(x, y, s);
    if (bp1 >= 2 && bp1 <= 6) //2nd condition
    {
        if (NumberOfZeroToOneTransitionFromP9(x, y, s) == 1)
        {
            if (even)
            {
                if (!((p2 && p4) && p8))
                {
                    if (!((p2 && p6) && p8))
                    {
                        return true;
                    }
                }
            }
            else
            {
                if (!((p2 && p4) && p6))
                {
                    if (!((p4 && p6) && p8))
                    {
                        return true;
                    }
                }
            }
        }
    }
    return false;
}

static int NumberOfZeroToOneTransitionFromP9(int x, int y, bool[][] s)
{
    bool p2 = s[x][y - 1];
    bool p3 = s[x + 1][y - 1];
    bool p4 = s[x + 1][y];
    bool p5 = s[x + 1][y + 1];
    bool p6 = s[x][y + 1];
    bool p7 = s[x - 1][y + 1];
    bool p8 = s[x - 1][y];
    bool p9 = s[x - 1][y - 1];

    int A = Convert.ToInt32((!p2  && p3 )) + Convert.ToInt32((!p3  && p4 )) +
            Convert.ToInt32((!p4  && p5 )) + Convert.ToInt32((!p5  && p6 )) +
            Convert.ToInt32((!p6  && p7 )) + Convert.ToInt32((!p7  && p8 )) +
            Convert.ToInt32((!p8  && p9 )) + Convert.ToInt32((!p9  && p2 ));
    return A;
}
static int NumberOfNonZeroNeighbors(int x, int y, bool[][] s)
{
    int count = 0;
    if (s[x - 1][y])     count++;
    if (s[x - 1][y + 1]) count++;
    if (s[x - 1][y - 1]) count++;
    if (s[x][y + 1])     count++;
    if (s[x][y - 1])     count++;
    if (s[x + 1][y])     count++;
    if (s[x + 1][y + 1]) count++;
    if (s[x + 1][y - 1]) count++;
    return count;
}

I have kept the original even flag, but call it by comparing a step number. And I have saved a few characters by using the bools directly..

Finally a function to get a deep copy of the nested 2d array:

public static T[][] ArrayClone<T>(T [][] A) 
     { return A.Select(a => a.ToArray()).ToArray(); }

This is how to call it, using two PictureBoxes:

pictureBox1.Image = Image.FromFile("D:\\RCdemo.png");
bool[][] t = Image2Bool(pictureBox1.Image);
t = ZhangSuenThinning(t);
pictureBox2.Image = Bool2Image(t);

I append a test image.

Answer 2

bwang22's answer is very slow. Try this instead:

public readonly struct ConnectivityData
{
    public readonly int[] N;
    public readonly int NumNeighbors;
    public readonly int NumChanges;

    public ConnectivityData(in int[] n, in int numNeighbors, in int numChanges)
    {
        N = n;
        NumNeighbors = numNeighbors;
        NumChanges = numChanges;
    }
}


public static void ZhangSuen(in HashSet<Pixel> pixels)
{        
    while (true)
    {
        // Pass #1:
        List<Pixel> mark1 = new List<Pixel>();
        foreach (Pixel p in pixels)
        {
            ConnectivityData conn = ComputeConnectivity(p, pixels);

            if (conn.NumNeighbors > 1 && 
                conn.NumNeighbors < 7 && 
                conn.NumChanges == 1 &&
                conn.N[0] * conn.N[2] * conn.N[4] == 0 &&
                conn.N[2] * conn.N[4] * conn.N[6] == 0)
            {
                mark1.Add(p);
            }
        }

        //delete all marked:
        foreach (Pixel p in mark1)
        {
            pixels.Remove(p);
        }

        // PASS #2:
        List<Pixel> mark2 = new List<Pixel>();
        foreach (Pixel p in pixels)
        {
            ConnectivityData conn = ComputeConnectivity(p, pixels);

            if (conn.NumNeighbors > 1 &&
                conn.NumNeighbors < 7 &&
                conn.NumChanges == 1 &&
                conn.N[0] * conn.N[2] * conn.N[6] == 0 &&
                conn.N[0] * conn.N[4] * conn.N[6] == 0)
            {
                mark2.Add(p);
            }
        }

        //delete all marked:
        foreach (Pixel p in mark2)
        {
            pixels.Remove(p);
        }

        if (mark1.Count == 0 && mark2.Count == 0)
        {
            break;
        }
    }

}

private static ConnectivityData ComputeConnectivity(
    in Pixel p, 
    in HashSet<Pixel> pixels)
{
    // calculate #neighbors and number of changes:
    int[] n = new int[8];
    if (pixels.Contains(new Pixel(p.X, p.Y - 1)))
    {
        n[0] = 1;
    }
    if (pixels.Contains(new Pixel(p.X + 1, p.Y - 1)))
    {
        n[1] = 1;
    }
    if (pixels.Contains(new Pixel(p.X + 1, p.Y)))
    {
        n[2] = 1;
    }
    if (pixels.Contains(new Pixel(p.X + 1, p.Y + 1)))
    {
        n[3] = 1;
    }
    if (pixels.Contains(new Pixel(p.X, p.Y + 1)))
    {
        n[4] = 1;
    }
    if (pixels.Contains(new Pixel(p.X - 1, p.Y + 1)))
    {
        n[5] = 1;
    }
    if (pixels.Contains(new Pixel(p.X - 1, p.Y)))
    {
        n[6] = 1;
    }
    if (pixels.Contains(new Pixel(p.X - 1, p.Y - 1)))
    {
        n[7] = 1;
    }

    return new ConnectivityData(
            n,
            n[0] + n[1] + n[2] + n[3] + n[4] + n[5] + n[6] + n[7],
            ComputeNumberOfChanges(n));
}

private static int ComputeNumberOfChanges(in int[] n)
{
    int numberOfChanges = 0;

    // Iterate over each location and see if it is has changed from 0 to 1:
    int current = n[0];
    for (int i = 1; i < 8; i++)
    {
        if (n[i] == 1 && current == 0)
        {
             numberOfChanges++;
        }
        current = n[i];
    }

    // Also consider the change over the discontinuity between n[7] and n[0]:
    if (n[0] == 1 && n[7] == 0)
    {
        numberOfChanges++;
    }

    return numberOfChanges;
}

To use:

From your Bitmap etc, create a hash set of type Pixel, (which contains all the black pixels you want to thin) eg:

public class Pixel
{
    public int X;
    public int Y;

    public Pixel(in int x, in int y)
    {
        X = x;
        Y = y;
    }

    public override bool Equals(object pixel)
    {
        Pixel b = pixel as Pixel;
        return X == b.X && Y == b.Y;
    }

    public override int GetHashCode()
    {
        //return (a.X << 2) ^ a.Y; // this is also commonly used as a pixel hash code
        return X * 100000 + Y; // a bit hacky [will fail if bitmap width is > 100000]
    }
}

...then call ZhangSuen(pixels). This will delete the appropriate pixels from the set.

Note that this method does not work perfectly on all images. It makes parts of some images disappear. Specifically, I am having problems with downward-right pointing diagonal lines of thickness around 11 pixels wide.

I am currently working on a way to improve this, but it performs better than the similar Staniford algorithm on most files I have tested it with (CAD files).

Answer 3

Here is my C# implementation of Zhang-Suen thinning algorithm

public static bool[][] ZhangSuenThinning(bool[][] s)
    {
        bool[][] temp = s;
        bool even = true;

        for (int a = 1; a < s.Length-1; a++)
        {
            for (int b = 1; b < s[0].Length-1; b++)
            {
                if (SuenThinningAlg(a, b, temp, even))
                {
                    temp[a][b] = false;
                }
                even = !even;
            }
        }

        return temp;
    }
static bool SuenThinningAlg(int x, int y, bool[][] s, bool even)
    {
        bool p2 = s[x][y - 1];
        bool p3 = s[x + 1][y - 1];
        bool p4 = s[x + 1][y];
        bool p5 = s[x + 1][y + 1];
        bool p6 = s[x][y + 1];
        bool p7 = s[x - 1][y + 1];
        bool p8 = s[x - 1][y];
        bool p9 = s[x - 1][y - 1];


            int bp1 = NumberOfNonZeroNeighbors(x, y, s);
            if (bp1 >= 2 && bp1 <= 6)//2nd condition
            {
                if (NumberOfZeroToOneTransitionFromP9(x, y, s) == 1)
                {
                    if (even)
                    {
                        if (!((p2 && p4) && p8))
                        {
                            if (!((p2 && p6) && p8))
                            {
                                return true;
                            }
                        }
                    }
                    else
                    {
                        if (!((p2 && p4) && p6))
                        {
                            if (!((p4 && p6) && p8))
                            {
                                return true;
                            }
                        }
                    }
                }
            }


        return false;
    }
    static int NumberOfZeroToOneTransitionFromP9(int x, int y, bool[][]s)
    {
        bool p2 = s[x][y - 1];
        bool p3 = s[x + 1][y - 1];
        bool p4 = s[x + 1][y];
        bool p5 = s[x + 1][y + 1];
        bool p6 = s[x][y + 1];
        bool p7 = s[x - 1][y + 1];
        bool p8 = s[x - 1][y];
        bool p9 = s[x - 1][y - 1];

        int A = Convert.ToInt32((p2 == false && p3 == true)) + Convert.ToInt32((p3 == false && p4 == true)) +
                 Convert.ToInt32((p4 == false && p5 == true)) + Convert.ToInt32((p5 == false && p6 == true)) +
                 Convert.ToInt32((p6 == false && p7 == true)) + Convert.ToInt32((p7 == false && p8 == true)) +
                 Convert.ToInt32((p8 == false && p9 == true)) + Convert.ToInt32((p9 == false && p2 == true));
        return A;
    }
    static int NumberOfNonZeroNeighbors(int x, int y, bool[][]s)
    {
        int count = 0;
        if (s[x-1][y])
            count++;
        if (s[x-1][y+1])
            count++;
        if (s[x-1][y-1])
            count++;
        if (s[x][y+1])
            count++;
        if (s[x][y-1])
            count++;
        if (s[x+1][y])
            count++;
        if (s[x+1][y+1])
            count++;
        if (s[x+1][y-1])
            count++;
        return count;
    }