Marching Square algorithm in Python

Question

The following Python source code is the implementation of Marching Square algorithm:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import cv2
from PIL import Image, ImageDraw

class Square():
    A = [0, 0]
    B = [0, 0]
    C = [0, 0]
    D = [0, 0]
    A_data = 0.0
    B_data = 0.0
    C_data = 0.0
    D_data = 0.0

    def GetCaseId(self, threshold):
        caseId = 0
        if (self.A_data >= threshold):
            caseId |= 1
        if (self.B_data >= threshold):
            caseId |= 2
        if (self.C_data >= threshold):
            caseId |= 4
        if (self.D_data >= threshold):
            caseId |= 8
            
        return caseId

    def GetLines(self, Threshold):
        lines = []
        caseId = self.GetCaseId(Threshold)

        if caseId in (0, 15):
            return []

        if caseId in (1, 14, 10):
            pX = (self.A[0] + self.B[0]) / 2
            pY = self.B[1]
            qX = self.D[0]
            qY = (self.A[1] + self.D[1]) / 2

            line = (pX, pY, qX, qY)

            lines.append(line)

        if caseId in (2, 13, 5):
            pX = (self.A[0] + self.B[0]) / 2
            pY = self.A[1]
            qX = self.C[0]
            qY = (self.A[1] + self.D[1]) / 2

            line = (pX, pY, qX, qY)

            lines.append(line)

        if caseId in (3, 12):
            pX = self.A[0]
            pY = (self.A[1] + self.D[1]) / 2
            qX = self.C[0]
            qY = (self.B[1] + self.C[1]) / 2

            line = (pX, pY, qX, qY)

            lines.append(line)

        if caseId in (4, 11, 10):
            pX = (self.C[0] + self.D[0]) / 2
            pY = self.D[1]
            qX = self.B[0]
            qY = (self.B[1] + self.C[1]) / 2

            line = (pX, pY, qX, qY)

            lines.append(line)

        elif caseId in (6, 9):
            pX = (self.A[0] + self.B[0]) / 2
            pY = self.A[1]
            qX = (self.C[0] + self.D[0]) / 2
            qY = self.C[1]

            line = (pX, pY, qX, qY)

            lines.append(line)

        elif caseId in (7, 8, 5):
            pX = (self.C[0] + self.D[0]) / 2
            pY = self.C[1]
            qX = self.A[0]
            qY = (self.A[1] + self.D[1]) / 2

            line = (pX, pY, qX, qY)

            lines.append(line)

        return lines

def marching_square(xVector, yVector, Data, threshold):
    linesList = []

    Height = len(Data)  # rows
    Width = len(Data[1])  # cols

    if ((Width == len(xVector)) and (Height == len(yVector))):
        squares = np.full((Height - 1, Width - 1), Square())

        sqHeight = squares.shape[0]  # rows count
        sqWidth = squares.shape[1]  # cols count

        for j in range(sqHeight):  # rows
            for i in range(sqWidth):  # cols
                a = Data[j + 1, i]
                b = Data[j + 1, i + 1]
                c = Data[j, i + 1]
                d = Data[j, i]
                A = [xVector[i], yVector[j + 1]]
                B = [xVector[i + 1], yVector[j + 1]]
                C = [xVector[i + 1], yVector[j]]
                D = [xVector[i], yVector[j]]

                squares[j, i].A_data = a
                squares[j, i].B_data = b
                squares[j, i].C_data = c
                squares[j, i].D_data = d

                squares[j, i].A = A
                squares[j, i].B = B
                squares[j, i].C = C
                squares[j, i].D = D

                list = squares[j, i].GetLines(threshold)

                linesList = linesList + list
    else:
        raise AssertionError

    return [linesList]

def main():
    x = [i for i in range(256)]
    y = [i for i in range(256)]
    
    example_l = [[0 for i in range(256)] for j in range(256)]
    
    for i in range(len(example_l)):
        for j in range(len(example_l[0])):
            example_l[i][j] = math.sin(i / 10.0)*math.cos(j / 10.0)
    example = np.array(example_l)

    im = Image.new('RGB', (256, 256), (128, 128, 128))

    collection = marching_square(x, y, example, 0.9)

    draw = ImageDraw.Draw(im)

    for ln in collection:
        for toup in ln:
            draw.line(toup, fill=(255, 255, 0), width=1)

    plt.axis("off")
    plt.imshow( im )
    plt.show()

if __name__ == '__main__':
    main()

Output:

My question is: Why does this source code generate circular patterns?

what is the mathematical explanation of this circular pattern?

Answer 1

You are applying the marching squares algorithm to a sample of the surface defined by:

F(x,y) = sin(x)*cos(y)

If you plot it, e.g. with google here

You get a surface that looks like eggs boxes, and your marching squares algorithm is finding the isolines ("lines following a single data level, or isovalue.") at 0.9, which are circles. You can imagine this as the intersection of that surface with a plane parallel to the XY plane and at Z = 0.9.