Python中的Marching Square算法

问题描述

以下Python源代码是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,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,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,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,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,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()

输出：

我的问题是：为什么此源代码会生成圆形图案？

这种圆形图案的数学解释是什么？

解决方法

您正在将行进平方算法应用于以下对象定义的表面样本：

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

如果将其绘制出来，例如与Google here

您得到的表面看起来像鸡蛋盒，行进的正方形算法正在寻找等值线（“遵循单个数据级别的线或等值线”），其值为0.9。您可以将其想象为该表面与平行于XY平面且Z = 0.9的平面的交点。