C++ 中的 Delaunay 三角剖分

Question

以下步骤假设您从两个点开始 - A 和 B - 并试图确定要用于形成的点 C 三角形：

一个。创建一个成员函数，该函数将确定给定点 C 是在由两点 A 和 B 形成的线的左侧还是右侧。提示：为此，取向量的叉积 两点 A 和 B 之间以及 A 和 C 之间的向量。由于叉积 与两个向量之间夹角的正弦成正比，对于 0 到 180 度之间的夹角，它将是一个正值（即，如果点 C 位于从 A 到 B 的直线的左侧）。

b.创建一个成员函数，该函数将确定给定点是否在由其他三个点形成的圆内（在这种情况下，三个点将是三角形的点，生成的圆将是外接圆）。提示：一个非常优雅的实现 函数可以在 Delaunay 三角测量的 Wikipedia 条目中找到。如果你借 这个实现，确保彻底测试它，引用它，并解释它是如何实现的 在你的实验报告中起作用。

c。创建一个成员函数，给定两个点，找到下一个 Delaunay 的点 三角形。这个函数可以搜索整个列表（不是最有效的实现， 但对于本实验来说已经足够了）并且很可能会调用 4a 和 4b 中定义的函数 以上。

d。创建调用函数的递归成员函数 Delaunay(Point A, Point B) 从上面的 4c 找到下一个 C 点。当找到 C 点时，这个函数应该 在三角形向量中插入一个新三角形，同时更新布尔变量 点列表。然后，该函数应使用点 A 和 C 递归调用自身，并再次使用点 C 和 B。确保该函数也有两个基本情况；一种是当找到的点 C 已经被使用时，在这种情况下应该添加三角形，但不应该发生递归调用，并且；另一种没有找到点 C 的情况（因为 A 和 B 形成数据集的外边缘）

我已经完成了代码，但似乎在调试过程中不断出现错误，我相信我已经将其范围缩小到读取功能，所以这是我的代码：

PointList::PointList()
{
totPt = 0;
totTri = 0;
}


void PointList::read(const char* filename)
{
FILE* file;
file = fopen ( filename, "r" );

fscanf ( file, "%d  \n", &totPt);
datapoint.resize( totPt );

for ( unsigned int i=0; i < datapoint.size(); i++ )
{
    fscanf ( file, " %d %lf %lf \n", &datapoint.at(i).ptnum, &datapoint.at(i).xCoord, &datapoint.at(i).yCoord  );
}

}

void PointList::TriPrint()
{
cout << "# of triangles created: " << triPoint.size() << endl;
cout << "Triangle # : Vertice1   Vertice2  Vertice3" << endl;

for ( unsigned int i = 0; i < triPoint.size() ; i++)
{
    cout << triPoint.at(i).triNum << " : " << triPoint.at(i).vert1.ptnum << "   " << triPoint.at(i).vert2.ptnum << "   " << triPoint.at(i).vert3.ptnum << endl;
}
return;
 }

//TASK 4a：将确定给定点C是在由两点a和b形成的线的左侧还是右侧

bool PointList::isRight( double a, double b, double c )
{
 Point A = datapoint.at(a-1);
Point B = datapoint.at(b-1);
Point C = datapoint.at(c-1);

//Cross product of AB and AC
Point AB;
Point AC;

 AB.xCoord = B.xCoord - A.xCoord;
AB.yCoord = B.yCoord - A.yCoord;

 AC.xCoord = C.xCoord - A.xCoord;
AC.yCoord = C.yCoord - A.yCoord;

 double det = (AB.xCoord*AC.yCoord) - (AC.xCoord*AB.yCoord);

if ( det >= 0 )
{  
    return true;
}
else
{
    return false;
}
}

//任务4b：判断给定点是否在其他三个点形成的圆内

bool PointList::inCircle ( double a, double b, double c, double d )
{
bool inCirc = false;
Point A = datapoint.at(a-1);
Point B = datapoint.at(b-1);
Point C = datapoint.at(c-1);
Point D = datapoint.at(d-1);

vector <vector<double>> Matrix;
Matrix.resize(3);
for ( int i = 0; i < 3 ; i++)
{
    Matrix.at(i).resize(3);
}

Matrix.at(0).at(0) = A.xCoord - D.xCoord;
Matrix.at(1).at(0) = B.xCoord - D.xCoord;
Matrix.at(2).at(0) = C.xCoord - D.xCoord;

Matrix.at(0).at(1) = A.yCoord - D.yCoord;
Matrix.at(1).at(1) = B.yCoord - D.yCoord;
Matrix.at(2).at(1) = C.yCoord - D.yCoord;

Matrix.at(0).at(2) = ((A.xCoord*A.xCoord) - (D.xCoord*D.xCoord)) + ((A.yCoord*A.yCoord) - (D.yCoord*D.yCoord));
Matrix.at(1).at(2) = ((B.xCoord*B.xCoord) - (D.xCoord*D.xCoord)) + ((B.yCoord*B.yCoord) - (D.yCoord*D.yCoord));
Matrix.at(2).at(2) = ((C.xCoord*C.xCoord) - (D.xCoord*D.xCoord)) + ((C.yCoord*C.yCoord) - (D.yCoord*D.yCoord));


double det = 0;

det = (Matrix.at(0).at(0) * Matrix.at(1).at(1) * Matrix.at(2).at(2)) + (Matrix.at(0).at(1) * Matrix.at(1).at(2) * Matrix.at(2).at(0)) + (Matrix.at(0).at(2) * Matrix.at(1).at(0) * Matrix.at(2).at(1));
det = det - (Matrix.at(2).at(0) * Matrix.at(1).at(1) * Matrix.at(0).at(2)) - (Matrix.at(2).at(1) * Matrix.at(1).at(2) * Matrix.at(0).at(0)) - (Matrix.at(2).at(2) * Matrix.at(1).at(0) * Matrix.at(0).at(1));

if ( det >= 0 )     // determinant is positive if and only if D lies inside the circumcircle
{
    inCirc = false;
}
else
{
    inCirc = true;
}

return inCirc;
}

//任务4c：给定两个点，找到下一个德劳内三角形的点

double PointList::nextDelaunay ( double a, double b )
{
bool incircle = false;
bool isright = false;
 bool noRight = false;
 int f = -1;


//check all points in file
for ( int i = 0; i < datapoint.size() ; i++)
{
    if ( i == a || i == b)
    {
        continue;
    }
     else
    {
         isright = isRight( a, b, i);   //checks to see if its to the right
        if(isright)
        {
            noRight = false;
            for ( int j = 1; j < datapoint.size(); j++ )
            {
                if ( j == a || j == b || j == i )
                                   {
                    continue;
                }
                else
                {
                    incircle = inCircle ( a, b, i, j );    
                                           //checks to see if point in circle
                     if(incircle)
                    {
                        break;
                    }
                }
            }
            if ( !incircle)
            {
                return i;
            }

        }
        else
        {
            continue;
        }
    }
}
if (noRight)
{
    return f;
 }

}

//任务4d：递归成员函数，Delaunay，上面4c的调用，寻找下一个点C

void PointList::Delaunay( int a, int b )
{
Triangle x;
Point A = datapoint.at(a - 1);
Point B = datapoint.at(b - 1);

int c = nextDelaunay(a,b);

if ( c == -1)
{
    return;
}
else
{
    Point C = datapoint.at(c-1);


    if ( C.usedPt == false)
    {
        x.vert1.ptnum = a;
        x.vert1.xCoord = A.xCoord;
        x.vert1.yCoord = A.yCoord;

        x.vert2.ptnum = b;
        x.vert2.xCoord = B.xCoord;
        x.vert2.yCoord = B.yCoord;

        x.vert3.ptnum = c;
        x.vert3.xCoord = C.xCoord;
        x.vert3.yCoord = C.yCoord;

        x.triNum = triPoint.size()+1;
        triPoint.push_back(x);
        datapoint.at(c-1).usedPt = true;

        Delaunay( a, c );
        Delaunay( c, b );

    return;
    }

    if ( C.usedPt == true )
    {
        x.vert1.ptnum = a;
        x.vert1.xCoord = A.xCoord;
        x.vert1.yCoord = A.yCoord;

        x.vert2.ptnum = b;
        x.vert2.xCoord = B.xCoord;
        x.vert2.yCoord = B.yCoord;

        x.vert3.ptnum = c;
        x.vert3.xCoord = C.xCoord;
        x.vert3.yCoord = C.yCoord;

        x.triNum = triPoint.size()+1;
        triPoint.push_back(x);

    return;
    }
}
}

Answer 1

当我为 ENGO 333 做同样的任务时，我可以告诉你，你的错误在于你的 nextDelaunay 函数。通过实现每一个函数来分别测试每一个函数。

Answer 2

有一个优秀的图书馆可以做你想做的事：http://www.vtk.org/。