Marching Square 算法中的插值

Interpolation in Marching Square Algorithm

Here是源代码的GitRepository

本文在“线性插值”部分讨论了当直线倾斜时如何对值进行插值。

例如,案例#2 的计算如下:

我实现了如下插值:

public class Square
{
    public Point A { get; set; }//bottom left point
    public Point B { get; set; }//bottom right point
    public Point C { get; set; }//top right point
    public Point D { get; set; }//top left point

    public double A_data { get; set; }//bottom left data
    public double B_data { get; set; }//bottom right data
    public double C_data { get; set; }//top roght data
    public double D_data { get; set; }//top left data

    public Square()
    {
        A = new Point();
        B = new Point();
        C = new Point();
        D = new Point();
    }

    private double GetCaseId(double threshold)
    {
        int caseId = 0;

        if (A_data >= threshold)
        {
            caseId |= 1;
        }
        if (B_data >= threshold)
        {
            caseId |= 2;
        }
        if (C_data >= threshold)
        {
            caseId |= 4;
        }
        if (D_data >= threshold)
        {
            caseId |= 8;
        }

        return caseId;
    }

    public List<Line> GetLines(double Threshold)
    {
        List<Line> linesList = new List<Line>();

        double caseId = GetCaseId(Threshold);

        if (caseId == 0) {/*do nothing*/ }
        if (caseId == 15) {/*do nothing*/ }

        if ((caseId == 1) || (caseId == 14))
        {
            double pX = B.X + (A.X - B.X) * ((1 - B_data) / (A_data - B_data));
            double pY = B.Y;
            Point p = new Point(pX, pY);

            double qX = D.X;
            double qY = D.Y + (A.Y - D.Y) * ((1 - D_data) / (A_data - D_data));
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }
        /*2==13*/
        if ((caseId == 2) || (caseId == 13))//B
        {
            double pX = A.X + (B.X - A.X) * ((1 - A_data) / (B_data - A_data));
            double pY = A.Y;
            Point p = new Point(pX, pY);

            double qX = C.X;
            double qY = C.Y + (B.Y - C.Y) * ((1 - C_data) / (B_data - C_data));
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }
        /*3==12*/
        if ((caseId == 3) || (caseId == 12))
        {
            double pX = A.X;
            double pY = A.Y + (D.Y - A.Y) * ((1 - A_data) / (D_data - A_data));
            Point p = new Point(pX, pY);

            double qX = C.X;
            double qY = C.Y + (B.Y - C.Y) * ((1 - C_data) / (B_data - C_data));
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }
        /*4==11*/
        if ((caseId == 4) || (caseId == 11))
        {
            double pX = D.X + (C.X - D.X) * ((1 - D_data) / (C_data - D_data));
            double pY = D.Y;
            Point p = new Point(pX, pY);

            double qX = B.X;
            double qY = B.Y + (C.Y - B.Y) * ((1 - B_data) / (C_data - B_data));
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }
        /*6==9*/
        if ((caseId == 6) || (caseId == 9))
        {
            double pX = A.X + (B.X - A.X) * ((1 - A_data) / (B_data - A_data));
            double pY = A.Y;
            Point p = new Point(pX, pY);

            double qX = C.X + (D.X - C.X) * ((1 - C_data) / (D_data - C_data));
            double qY = C.Y;
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }

        /*7==8*/
        if ((caseId == 7) || (caseId == 8))
        {
            double pX = C.X + (D.X - C.X) * ((1 - C_data) / (D_data - C_data));
            double pY = C.Y;
            Point p = new Point(pX, pY);

            double qX = A.X;
            double qY = A.Y + (D.Y - A.Y) * ((1 - A_data) / (D_data - A_data));
            Point q = new Point(qX, qY);

            Line line = new Line(p, q);

            linesList.Add(line);
        }

        /*ambiguous case*/
        if (caseId == 5)
        {
            double pX1 = A.X + (B.X - A.X) * ((1 - A_data) / (B_data - A_data));
            double pY1 = A.Y;
            Point p1 = new Point(pX1, pY1);
            double qX1 = C.X;
            double qY1 = C.Y + (B.Y - C.Y) * ((1 - C_data) / (B_data - C_data));
            Point q1 = new Point(qX1, qY1);
            Line line1 = new Line(p1, q1);

            double pX2 = C.X + (D.X - C.X) * ((1 - C_data) / (D_data - C_data));
            double pY2 = C.Y;
            Point p2 = new Point(pX2, pY2);
            double qX2 = A.X;
            double qY2 = A.Y + (D.Y - A.Y) * ((1 - A_data) / (D_data - A_data));
            Point q2 = new Point(qX2, qY2);
            Line line2 = new Line(p2, q2);

            linesList.Add(line1);
            linesList.Add(line2);
        }
        if (caseId == 10)
        {
            double pX1 = B.X + (A.X - B.X) * ((1 - B_data) / (A_data - B_data));
            double pY1 = B.Y;
            Point p1 = new Point(pX1, pY1);
            double qX1 = D.X;
            double qY1 = D.Y + (A.Y - D.Y) * ((1 - D_data) / (A_data - D_data));
            Point q1 = new Point(qX1, qY1);
            Line line1 = new Line(p1, q1);

            double pX2 = D.X + (C.X - D.X) * ((1 - D_data) / (C_data - D_data));
            double pY2 = D.Y;
            Point p2 = new Point(pX2, pY2);
            double qX2 = B.X;
            double qY2 = B.Y + (C.Y - B.Y) * ((1 - B_data) / (C_data - B_data));
            Point q2 = new Point(qX2, qY2);
            Line line2 = new Line(p2, q2);

            linesList.Add(line1);
            linesList.Add(line2);
        }

        return linesList;
    }

但是,这不能正常工作。

任何人都可以检查插值部分并告诉我哪里出了问题吗?

我检查了你的 repo,错误不在你发布的代码中:|

实际错误在于您在 MarchingSquare.cs

中初始化数据的方式

替换以下内容:

double a = Data[j + 1, i];
double b = Data[j + 1, i + 1];
double c = Data[j, i + 1];
double d = Data[j, i];

Point A = new Point(xVector[i], yVector[j + 1]);//A
Point B = new Point(xVector[i + 1], yVector[j + 1]);//B
Point C = new Point(xVector[i + 1], yVector[j]);//C
Point D = new Point(xVector[i], yVector[j]);//D

与:

double a = Data[j, i];
double b = Data[j, i + 1];
double c = Data[j+1, i + 1];
double d = Data[j+1, i];

Point A = new Point(j, i);//A
Point B = new Point(j, i+1);//B
Point C = new Point(j+1, i+1);//C
Point D = new Point(j+1,i);//D

你会得到想要的结果(可能是转置)

编辑:我想我把它调换了:Plot[Sin[x/100.0] * Cos[y/100.0]] = 0.9, x = 0 to 250, y = 0 to 500

看了一点插值代码:

  • 对案例使用枚举而不是神秘常量
  • 不要将 caseId 转换为 double
  • 提取辅助方法的插值!
using System;
using System.Collections.Generic;

namespace G__Marching_Sqaure
{
    public class Square
    {
        public Point A { get; set; } //bottom left point
        public Point B { get; set; } //bottom right point
        public Point C { get; set; } //top right point
        public Point D { get; set; } //top left point

        public double A_data { get; set; } //bottom left data
        public double B_data { get; set; } //bottom right data
        public double C_data { get; set; } //top roght data
        public double D_data { get; set; } //top left data

        public Square()
        {
            A = new Point();
            B = new Point();
            C = new Point();
            D = new Point();
        }

        private LineShapes GetCaseId(double threshold)
        {
            int caseId = 0;

            if (A_data >= threshold)
            {
                caseId |= 1;
            }

            if (B_data >= threshold)
            {
                caseId |= 2;
            }

            if (C_data >= threshold)
            {
                caseId |= 4;
            }

            if (D_data >= threshold)
            {
                caseId |= 8;
            }

            return (LineShapes)caseId;
        }

        public List<Line> GetLines(double Threshold)
        {
            List<Line> linesList = new List<Line>();

            LineShapes caseId = GetCaseId(Threshold);

            if (caseId == LineShapes.Empty)
            {
                /*do nothing*/
            }

            if (caseId == LineShapes.All)
            {
                /*do nothing*/
            }

            if ((caseId == LineShapes.BottomLeft) || (caseId == LineShapes.AllButButtomLeft))
            {
                var p = InterpolateHorizonal(B, A, B_data, A_data);
                var q = InterpolateVertical(D, A, D_data, A_data);
                Line line = new Line(p, q);
                linesList.Add(line);
            }

            /*2==13*/
            if ((caseId == LineShapes.BottomRight) || (caseId == LineShapes.AllButButtomRight)) //B
            {
                var p = InterpolateHorizonal(A, B, A_data, B_data);
                var q = InterpolateVertical(C, B, C_data, B_data);
                Line line = new Line(p, q);

                linesList.Add(line);
            }

            /*3==12*/
            if ((caseId == LineShapes.Bottom) || (caseId == LineShapes.Top))
            {
                // interpolate vertical
                var p = InterpolateVertical(A, D, A_data, D_data);
                var q = InterpolateVertical(C, B, C_data, B_data);
                Line line = new Line(p, q);
                linesList.Add(line);
            }

            /*4==11*/
            if ((caseId == LineShapes.TopRight) || (caseId == LineShapes.AllButTopRight))
            {
                var p = InterpolateHorizonal(D, C, D_data, C_data);
                var q = InterpolateVertical(B, C, B_data, C_data);
                Line line = new Line(p, q);
                linesList.Add(line);
            }

            /*6==9*/
            if ((caseId == LineShapes.Right) || (caseId == LineShapes.Left))
            {
                var p = InterpolateHorizonal(A, B, A_data, B_data);
                var q = InterpolateHorizonal(C, D, C_data, D_data);
                Line line = new Line(p, q);

                linesList.Add(line);
            }

            /*7==8*/
            if ((caseId == LineShapes.AllButTopLeft) || (caseId == LineShapes.TopLeft))
            {
                var p = InterpolateHorizonal(C, D, C_data, D_data);
                var q = InterpolateVertical(A, D, A_data, D_data);
                Line line = new Line(p, q);
                linesList.Add(line);
            }

            /*ambiguous case*/
            if (caseId == LineShapes.TopRightBottomLeft)
            {
                var p1 = InterpolateHorizonal(A, B, A_data, B_data);
                var q1 = InterpolateVertical(C, B, C_data, B_data);
                Line line1 = new Line(p1, q1);

                var p2 = InterpolateHorizonal(C, D, C_data, D_data);
                var q2 = InterpolateVertical(A, D, A_data, D_data);
                Line line2 = new Line(p2, q2);

                linesList.Add(line1);
                linesList.Add(line2);
            }

            if (caseId == LineShapes.TopLeftBottomRight)
            {
                var p1 = InterpolateHorizonal(B, A, B_data, A_data);
                var q1 = InterpolateVertical(D, A, D_data, A_data);
                Line line1 = new Line(p1, q1);

                var p2 = InterpolateHorizonal(D, C, D_data, C_data);
                var q2 = InterpolateVertical(B, C, B_data, C_data);
                Line line2 = new Line(p2, q2);

                linesList.Add(line1);
                linesList.Add(line2);
            }

            return linesList;
        }

        private static Point InterpolateVertical(Point point, Point point1, double data, double data1)
        {
            double qX = point.X;
            double qY = point.Y + (point1.Y - point.Y) * ((1 - data) / (data1 - data));
            Point q = new Point(qX, qY);
            return q;
        }


        private static Point InterpolateHorizonal(Point start, Point end, double startForce, double endForce)
        {
            double pX = start.X + (end.X - start.X) * ((1 - startForce) / (endForce - startForce));
            double pY = start.Y;
            Point p = new Point(pX, pY);
            return p;
        }
    }

    internal enum LineShapes
    {
        Empty = 0,
        //  ○----○
        //  |    |
        //  |    |
        //  ○----○

        BottomLeft = 1,
        //  ○----○
        //  |    |
        //  |    |
        //  ●----○

        BottomRight = 2,
        //  ○----○
        //  |    |
        //  |    |
        //  ○----●

        Bottom = 3,
        //  ○----○
        //  |    |
        //  |    |
        //  ●----●

        TopRight = 4,
        //  ○----●
        //  |    |
        //  |    |
        //  ○----○

        TopRightBottomLeft = 5,
        //  ○----●
        //  |    |
        //  |    |
        //  ●----○

        Right = 6,
        //  ○----●
        //  |    |
        //  |    |
        //  ○----●

        AllButTopLeft = 7,
        //  ○----●
        //  |    |
        //  |    |
        //  ●----●

        TopLeft = 8,
        //  ●----○
        //  |    |
        //  |    |
        //  ○----○

        Left = 9,
        //  ●----○
        //  |    |
        //  |    |
        //  ●----○

        TopLeftBottomRight = 10,
        //  ●----○
        //  |    |
        //  |    |
        //  ○----●

        AllButTopRight = 11,
        //  ●----○
        //  |    |
        //  |    |
        //  ●----●

        Top = 12,
        //  ●----●
        //  |    |
        //  |    |
        //  ○----○

        AllButButtomRight = 13,
        //  ●----●
        //  |    |
        //  |    |
        //  ●----○

        AllButButtomLeft = 14,
        //  ●----●
        //  |    |
        //  |    |
        //  ○----●

        All = 15,
        //  ●----●
        //  |    |
        //  |    |
        //  ●----●
    }
}

另外,清理package.json

更新

再看文章,发现插值公式中的值1就是阈值

如果您使用任意阈值检查案例但不修改插值公式(因此在将阈值设置为 1 的情况下进行插值),您可能会得到位于给定正方形之外的点!

更新的插值例程:

private static Point InterpolateVertical(Point start, Point end, double startForce, double endForce, double threshold)
{
   
    double a = ((threshold - startForce) / (endForce - startForce));
    double qX = start.X;
    double qY = start.Y + (end.Y - start.Y) * a;
    Point q = new Point(qX, qY);
    return q;
}


private static Point InterpolateHorizonal(Point start, Point end, double startForce, double endForce, double threshold)
{

    double a = ((threshold - startForce) / (endForce - startForce));
    double pX = start.X + (end.X - start.X) * a;
    double pY = start.Y;
    Point p = new Point(pX, pY);
    return p;
}

您正确地注意到,在初始示例中,您的网格过于密集,因此 运行 每个像素的测试同样有效。

我检查了你的 sin-cos 示例,你仍在对每个点进行采样(正方形的宽度和高度等于 1),因此行进点算法没有真正的好处。

检查以下初始化代码,网格中有较大的方块(分辨率 x 分辨率)

int height = 1000;
int width = height;
int resolution = 20;
double threshold = 0.9;

int xLen = width / resolution;
int yLen = height / resolution;
int[] x = new int[xLen];
int[] y = new int[yLen];
for (int i = 0; i < xLen; i++)
{
    x[i] = i * resolution;
}

for (int j = 0; j < yLen; j++)
{
    y[j] = j * resolution;
}

/////////////////////SIN-COS/////////////////////////////
           
double[,] example = new double[xLen, yLen];
for (int j = 0; j < yLen; j++)
{
    for (int i = 0; i < xLen; i++)
    {
        double example_l = Math.Sin(i * resolution /100.0 ) * Math.Cos(j * resolution /100.0);
        example[j, i] = example_l;
    }
}

调整分辨率(例如尝试分辨率 = 10)以获得更平滑的形状。