First of all, we must reverse the original image, because all the operations after the algorithm take 0 as the foreground and 1 as the background.

The 8-near-neighbor definition of
center pixel X_1 (x, y) is shown below:

Consider the following two steps

• Step 1: Perform the grating scan and mark all the pixels that meet the following 5 conditions:
  • This is a black pixel;
  • clockwise viewing X2, X3, …, X9, X2, the number of changes from 0 to 1 is only 1;
  • x2, x3, …, x9, 1 of 1 or less than 2 or less;

    At least 1 of

  • x2, x4, and x6 is 1;
  • x4, x6, x8, at least one is 1;
    All pixels that meet the conditions are 1
• Step 2: Perform the grating scan and mark all the pixels that meet the following 5 conditions:
  • This is a black pixel;
  • Clockwise to view X2, X3, …, X9, X2, the number of changes from 0 to 1 is only 1;
  • x2, x3, …, x9, 1 of 1 or less than 2 or less;

    At least 1 of

  • x2, X4, X8 is 1;

    At least 1 of

  • x2, X6, X8 is 1;
    All pixels that meet the conditions are 1
• Repeat step 1 and step 2 until there is no change.
  python implementation:

import cv2
import numpy as np
import matplotlib.pyplot as plt

# Zhang Suen thining algorythm
def Zhang_Suen_thining(img):
    # get shape
    H, W, C = img.shape

    # prepare out image
    out = np.zeros((H, W), dtype=np.int)
    out[img[..., 0] > 0] = 1

    # inverse
    out = 1 - out

    while True:
        s1 = []
        s2 = []

        # step 1 ( rasta scan )
        for y in range(1, H-1):
            for x in range(1, W-1):
                
                # condition 1
                if out[y, x] > 0:
                    continue

                # condition 2
                f1 = 0
                if (out[y-1, x+1] - out[y-1, x]) == 1:
                    f1 += 1
                if (out[y, x+1] - out[y-1, x+1]) == 1:
                    f1 += 1
                if (out[y+1, x+1] - out[y, x+1]) == 1:
                    f1 += 1
                if (out[y+1, x] - out[y+1,x+1]) == 1:
                    f1 += 1
                if (out[y+1, x-1] - out[y+1, x]) == 1:
                    f1 += 1
                if (out[y, x-1] - out[y+1, x-1]) == 1:
                    f1 += 1
                if (out[y-1, x-1] - out[y, x-1]) == 1:
                    f1 += 1
                if (out[y-1, x] - out[y-1, x-1]) == 1:
                    f1 += 1

                if f1 != 1:
                    continue
                    
                # condition 3
                f2 = np.sum(out[y-1:y+2, x-1:x+2])
                if f2 < 2 or f2 > 6:
                    continue
                
                # condition 4
                # x2 x4 x6
                if (out[y-1, x] + out[y, x+1] + out[y+1, x]) < 1 :
                    continue

                # condition 5
                # x4 x6 x8
                if (out[y, x+1] + out[y+1, x] + out[y, x-1]) < 1 :
                    continue
                    
                s1.append([y, x])

        for v in s1:
            out[v[0], v[1]] = 1

        # step 2 ( rasta scan )
        for y in range(1, H-1):
            for x in range(1, W-1):
                
                # condition 1
                if out[y, x] > 0:
                    continue

                # condition 2
                f1 = 0
                if (out[y-1, x+1] - out[y-1, x]) == 1:
                    f1 += 1
                if (out[y, x+1] - out[y-1, x+1]) == 1:
                    f1 += 1
                if (out[y+1, x+1] - out[y, x+1]) == 1:
                    f1 += 1
                if (out[y+1, x] - out[y+1,x+1]) == 1:
                    f1 += 1
                if (out[y+1, x-1] - out[y+1, x]) == 1:
                    f1 += 1
                if (out[y, x-1] - out[y+1, x-1]) == 1:
                    f1 += 1
                if (out[y-1, x-1] - out[y, x-1]) == 1:
                    f1 += 1
                if (out[y-1, x] - out[y-1, x-1]) == 1:
                    f1 += 1

                if f1 != 1:
                    continue
                    
                # condition 3
                f2 = np.sum(out[y-1:y+2, x-1:x+2])
                if f2 < 2 or f2 > 6:
                    continue
                
                # condition 4
                # x2 x4 x8
                if (out[y-1, x] + out[y, x+1] + out[y, x-1]) < 1 :
                    continue

                # condition 5
                # x2 x6 x8
                if (out[y-1, x] + out[y+1, x] + out[y, x-1]) < 1 :
                    continue
                    
                s2.append([y, x])

        for v in s2:
            out[v[0], v[1]] = 1

        # if not any pixel is changed
        if len(s1) < 1 and len(s2) < 1:
            break

    out = 1 - out
    out = out.astype(np.uint8) * 255

    return out


# Read image
img = cv2.imread("../thin.png").astype(np.float32)

# Zhang Suen thining
out = Zhang_Suen_thining(img)


# Save result
cv2.imwrite("out.png", out)
cv2.imshow("result", out)
cv2.waitKey(0)
cv2.destroyAllWindows()

Experimental results: