Interpolation in the Marching Square Algorithm

Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
Sign up Log in
1 reply · 2K views
user366312
Gold Member
Messages
88
Reaction score
3
TL;DR
I have implemented the Marching Square algorithm. But, it is not giving the correct output.
The cause is probably has something to do with the interpolation part.
In the "Linear Interpolation" section This article discusses how to interpolate the values when the lines are oblique.

For example, for Case#2 it has the following calculation:

tiEsh.png

HAQOb.png


I have written the following source code to implement the Marching Square algorithm:

[CODE lang="python" title="Marching Square"]import numpy as np

from PIL import Image, ImageDraw

import math

import pandas as pd

import copy

import cv2
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):

linesList = [];
caseId = self.GetCaseId(Threshold);
if (caseId == 0):

pass;

if (caseId == 15):

pass;
if ((caseId == 1) or (caseId == 14)):

pX = self.B[0] + (self.A[0] - self.B[0]) * ((1 - self.B_data) / (self.A_data - self.B_data));

pY = self.B[1];

qX = self.D[0];

qY = self.D[1] + (self.A[1] - self.D[1]) * ((1 - self.D_data) / (self.A_data - self.D_data));
line = (pX, pY, qX, qY);
linesList.append(line);
if ((caseId == 2) or (caseId == 13)):

pX = self.A[0] + (self.B[0] - self.A[0]) * ((1 - self.A_data) / (self.B_data - self.A_data));

pY = self.A[1];

qX = self.C[0];

qY = self.C[1] + (self.B[1] - self.C[1]) * ((1 - self.C_data) / (self.B_data - self.C_data));
line = (pX, pY, qX, qY);
linesList.append(line);
if ((caseId == 3) or (caseId == 12)):

pX = self.A[0];

pY = self.A[1] + (self.D[1] - self.A[1]) * ((1 - self.A_data) / (self.D_data - self.A_data));

qX = self.C[0];

qY = self.C[1] + (self.B[1] - self.C[1]) * ((1 - self.C_data) / (self.B_data - self.C_data));
line = (pX, pY, qX, qY);
linesList.append(line);
if ((caseId == 4) or (caseId == 11)):

pX = self.D[0] + (self.C[0] - self.D[0]) * ((1 - self.D_data) / (self.C_data - self.D_data));

pY = self.D[1];

qX = self.B[0];

qY = self.B[1] + (self.C[1] - self.B[1]) * ((1 - self.B_data) / (self.C_data - self.B_data));
line = (pX, pY, qX, qY);
linesList.append(line);
if ((caseId == 6) or (caseId == 9)):

pX = self.A[0] + (self.B[0] - self.A[0]) * ((1 - self.A_data) / (self.B_data - self.A_data));

pY = self.A[1];

qX = self.C[0] + (self.D[0] - self.C[0]) * ((1 - self.C_data) / (self.D_data - self.C_data));

qY = self.C[1];
line = (pX, pY, qX, qY);
linesList.append(line);
if ((caseId == 7) or (caseId == 8)):

pX = self.C[0] + (self.D[0] - self.C[0]) * ((1 - self.C_data) / (self.D_data - self.C_data));

pY = self.C[1];

qX = self.A[0];

qY = self.A[1] + (self.D[1] - self.A[1]) * ((1 - self.A_data) / (self.D_data - self.A_data));
line = (pX, pY, qX, qY);
linesList.append(line);
if (caseId == 5):

pX1 = self.A[0] + (self.B[0] - self.A[0]) * ((1 - self.A_data) / (self.B_data - self.A_data));

pY1 = self.A[1];

qX1 = self.C[0];

qY1 = self.C[1] + (self.B[1] - self.C[1]) * ((1 - self.C_data) / (self.B_data - self.C_data));
line1 = (pX1, pY1, qX1, qY1);
pX2 = self.C[0] + (self.D[0] - self.C[0]) * ((1 - self.C_data) / (self.D_data - self.C_data));

pY2 = self.C[1];

qX2 = self.A[0];

qY2 = self.A[1] + (self.D[1] - self.A[1]) * ((1 - self.A_data) / (self.D_data - self.A_data));
line2 = (pX2, pY2, qX2, qY2);
linesList.append(line1);

linesList.append(line2);
if (caseId == 10):

pX1 = self.B[0] + (self.A[0] - self.B[0]) * ((1 - self.B_data) / (self.A_data - self.B_data));

pY1 = self.B[1];

qX1 = self.D[0];

qY1 = self.D[1] + (self.A[1] - self.D[1]) * ((1 - self.D_data) / (self.A_data - self.D_data));
line1 = (pX1, pY1, qX1, qY1);
pX2 = self.D[0] + (self.C[0] - self.D[0]) * ((1 - self.D_data) / (self.C_data - self.D_data));

pY2 = self.D[1];

qX2 = self.B[0];

qY2 = self.B[1] + (self.C[1] - self.B[1]) * ((1 - self.B_data) / (self.C_data - self.B_data));
line2 = (pX2, pY2, qX2, qY2);
linesList.append(line1);

linesList.append(line2);
return linesList;def marching_square(xVector, yVector, Data, threshold):

linesList = [];
# Height = Data.shape[0];#rows

# Width = Data.shape[1];#cols
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];
squares[j, i].A_data = a;

squares[j, i].B_data = b;

squares[j, i].C_data = c;

squares[j, i].D_data = d;
A = [xVector, yVector[j + 1]];

B = [xVector[i + 1], yVector[j + 1]];

C = [xVector[i + 1], yVector[j]];

D = [xVector, yVector[j]];
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];[/CODE]

The following main() function gives the correct output:

[CODE lang="python" title="main()"]def main():

example = np.array([

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],

[0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],

[0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1],

[0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1],

[0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0]

]);
x = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100];

y = [0, 10, 20, 30, 40];
collection = marching_square(x, y, example, 1);
for ln in collection:

for toup in ln:

draw.line(toup, fill=(255, 255, 0), width=5)
im.save('output.jpg', quality=95)[/CODE]
HQ7oc.png


The following is not giving the correct output:

[CODE lang="python" title="main2()"]def main2():
x = [i for i in range(628)];

y = [i for i in range(628)];

example_l = [[0 for i in range(628)] for j in range(628)]

f = open("o.txt","w")

for i in range(len(example_l)):

for j in range(len(example_l[0])):

example_l[j] = math.sin(i / 100.0)*math.cos(j / 100.0)

print(i,j,example_l[j], file=f)

example = np.array(example_l)

f.close()
print("filled")
collection = marching_square(x, y, example, 0.9);
print("done")

I am = Image.new('RGB', (628, 628), (128, 128, 128))

draw = ImageDraw.Draw(im)
for ln in collection:

for toup in ln:

draw.line(toup, fill=(255, 255, 0), width=1)
im.save('output.jpg', quality=95)[/CODE]

0iqnA.png


The expected output should be probably as follows:

APQbx.png


I think the interpolation is not working correctly.

(https://github.com/mns-csharp/MarchingSquare is the C# GitRepository of the source code)
 
Physics news on Phys.org
This is 150 lines of code, without comments, with unnecessary empty lines that make scrolling through the code more annoying. I suggest you do more debugging yourself first. Disable more cases to identify the cases that give wrong lines. Then look at individual lines in the problematic cases. Most likely the calculated end points are wrong. This could be an incorrect formula, rounding errors, or python treating things as integer where you expect floats.
Your expectation has x and y swapped relative to the output, but that's probably a minor concern.
 
  • Like
Likes   Reactions: Filip Larsen