Interpolation in the Marching Square Algorithm

In summary: MarchingSquare(self, x, y, lines): x_start, y_start, x_end, y_end = x, y, x_start + lines[0], y_start + lines[1], x_end + lines[2], y_end + lines[3]; dX = (x_start - x_end) / lines[0];dY = (y_start - y_end) / lines[1];
  • #1
user366312
Gold Member
89
3
TL;DR Summary
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:

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[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 = 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];

The following main() function gives the correct output:

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)
HQ7oc.png


The following is not giving the correct output:

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[i][j] = math.sin(i / 100.0)*math.cos(j / 100.0)

            print(i,j,example_l[i][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)

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)
 
Technology news on Phys.org
  • #2
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 Filip Larsen

1. What is interpolation in the Marching Square Algorithm?

Interpolation in the Marching Square Algorithm is a method used to estimate values within a set of data points. In the context of the Marching Square Algorithm, it is used to generate smooth contours between grid cells by calculating the location of intersection points between contour lines and grid edges.

2. Why is interpolation important in the Marching Square Algorithm?

Interpolation is important in the Marching Square Algorithm because it allows for the creation of continuous and visually appealing contour lines between grid cells. Without interpolation, the resulting contours may appear jagged and less accurate.

3. What types of interpolation are commonly used in the Marching Square Algorithm?

The two most commonly used types of interpolation in the Marching Square Algorithm are linear and quadratic interpolation. Linear interpolation calculates the location of intersection points by taking a linear average of the values at the grid corners, while quadratic interpolation takes into account the slope and curvature of the contour lines.

4. How does interpolation impact the accuracy of the Marching Square Algorithm?

The accuracy of the Marching Square Algorithm is greatly improved with the use of interpolation. By estimating values between grid points, it creates a more precise representation of the data and produces smoother and more visually appealing contours.

5. Can interpolation be used in other algorithms besides the Marching Square Algorithm?

Yes, interpolation is a commonly used method in various algorithms and fields of study, such as computer graphics, image processing, and data analysis. It is a powerful tool for estimating values within a set of data points and is widely used in scientific and technical fields.

Similar threads

Replies
28
Views
2K
Replies
25
Views
2K
Back
Top