Interpolation in the Marching Square Algorithm

[user366312]

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.

• Marching Square - article

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:

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]

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]

The expected output should be probably as follows:

I think the interpolation is not working correctly.

(https://github.com/mns-csharp/MarchingSquare is the C# GitRepository of the source code)

 

[mfb]

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.