行进平方算法中的插值

问题描述

Here是源代码的GitRepository。

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

例如,对于案例2,它具有以下计算:

enter image description here

enter image description here

我实现了如下插值:

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,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,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,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,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,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,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,q2);

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

        return linesList;
    }

但是,这不能正常工作。

谁能检查插值部分并告诉我哪里出了问题?

解决方法

我检查了您的仓库,错误不在您发布的代码中:|

实际错误在于您如何在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转换为两倍
  • 将插值提取为辅助方法!
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,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,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,D_data);
                var q = InterpolateVertical(C,q);
                linesList.Add(line);
            }

            /*4==11*/
            if ((caseId == LineShapes.TopRight) || (caseId == LineShapes.AllButTopRight))
            {
                var p = InterpolateHorizonal(D,C,C_data);
                var q = InterpolateVertical(B,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_data);
                var q = InterpolateHorizonal(C,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_data);
                var q = InterpolateVertical(A,q);
                linesList.Add(line);
            }

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

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

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

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

                var p2 = InterpolateHorizonal(D,C_data);
                var q2 = InterpolateVertical(B,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,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,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;
    }
}

以分辨率播放(例如,尝试使用resolution = 10)以获得更平滑的形状。