python - Python 中的 Marching Square 算法

Question

以下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, 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()

输出:

我的问题是:为什么这个源代码会生成圆形图案？

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

Answer 1

您正在将行进方 block 算法应用于由以下定义的表面样本:

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

如果你绘制它，例如与谷歌 here

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