vtk学习笔记 --- 判断三角形相交

原文：/content/3091992.html

在使用三角网连接矿体的时候，需要判断当前连接的三角形和已经连接的三角形是否相交，所以，就需要进行三角形相交判断。

看了一些算法的文章，两个三角形相交的判断规则大体如下：

假设这两个三角形为A(a1,a2,a3),B(b1,b2,b3)，三角形A所在的平面为PA，法向量为NA，三角形B所在的平面为PB，法向量为NB。

1、将三角形A的所有顶点投影到平面PB上，投影得到的点为proja1,proja2,proja3

2、计算三角形A的所有顶点到平面PB的有符号距离：

|a1| = |a1-proja1|

|a2| = |a2-proja2|

|a3| = |a3-proja3|

其中，符号取决于 顶点与投影点构成的向量与平面PB法向量的方向是否一致，如果一致，则大于0，不一致则小于0

3、根据有符号距离来判断三角形相交情况：

a、如果三个有符号距离的符号一致（都大于0或者都小于0），说明两个三角形不相交

b、如果一个符号距离为0，另外两个符号距离乘积大于0，则说明两个三角形相交于一个顶点

c、 如果一个符号距离为0，另外两个符号距离乘积小于0，则说明符号距离为0的顶点位于平面PB内，而另外两个顶点位于平面PB两侧，需要计算交线来判断是否相交。

d、如果一个符号距离和另外两个符号距离的符号相反，说明一个顶点位于平面PB的一侧，两外两个顶点位于平面PB的另外一侧，需要计算交线来判断是否相交。

4、根据上面的符号距离，来计算三角形A与平面PB的交点，然后得到两个三角形与对应平面的交线，最后通过判断交线是否相交或者重合来判断三角形是否相交。

具体的C++代码如下：

Cpp代码


#include "vtkPolyData.h"

#include "vtkPoints.h"

#include "vtkTriangle.h"

#include "vtkCellArray.h"

#include "vtkActor.h"

#include "vtkPolyDataMapper.h"

#include "vtkRenderer.h"

#include "vtkRenderWindow.h"

#include "vtkRenderWindowInteractor.h"

#include "vtkMath.h"

#include "vtkPlane.h"

#include "vtkInteractorStyleTrackballCamera.h"

#include "vtkLine.h"

using namespace std;

//准备测试数据

vtkPolyData* prepareData(){

//交线位于其中一个三角形内部

//triangle 1

//vtkPoints* pts = vtkPoints::New();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.5 ,0.1,-1);

//pts->InsertNextPoint(0.5,0,1);

//pts->InsertNextPoint(0,0.5,1);

//部分相交，两个三角形相互咬住对方

//vtkPoints* pts = vtkPoints::New();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.25,0.25,0.5);

//pts->InsertNextPoint(0.25,-0.5,0.5);

//pts->InsertNextPoint(0,0,-0.5);

//边重合，不作为三角形相交

//vtkPoints* pts = vtkPoints::New();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.25,0.25,0.5);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

//一个点位于三角形内部，另外2个点位于两侧，有交线

//vtkPoints* pts = vtkPoints::New();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.25,0.25,0);

//pts->InsertNextPoint(0.45,0.45,1);

//pts->InsertNextPoint(0.45,0.45,-1);

//一个点位于三角形所在面内，另外2个点位于两侧，无交线

//vtkPoints* pts = vtkPoints::New();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.75,0.75,0);

//pts->InsertNextPoint(1.5,1.5,1);

//pts->InsertNextPoint(1.5,1.5,-1);

//一个点位于三角形所在面内，另外2个点位于同侧，有交点无交线

vtkPoints* pts = vtkPoints::New();

pts->InsertNextPoint(0,0,0);

pts->InsertNextPoint(1,0,0);

pts->InsertNextPoint(0,1,0);

//triangle 2

pts->InsertNextPoint(0.25,0.25,0);

pts->InsertNextPoint(1,0,1);

pts->InsertNextPoint(0,1,1);

//两个点位于另外一个三角形内部，另外一个点位于外面

//vtkPoints* pts = vtkPoints::New();

//pts->SetDataTypeToDouble();

//pts->InsertNextPoint(0,0,0);

//pts->InsertNextPoint(1,0,0);

//pts->InsertNextPoint(0,1,0);

////triangle 2

//pts->InsertNextPoint(0.25,0.25,1);

//pts->InsertNextPoint(0.1,0.5,0);

//pts->InsertNextPoint(0.5,0.1,0);

vtkTriangle* triangle1 = vtkTriangle::New();

triangle1->GetPointIds()->SetNumberOfIds(3);

triangle1->GetPointIds()->SetId(0, 0);

triangle1->GetPointIds()->SetId(1, 1);

triangle1->GetPointIds()->SetId(2, 2);

vtkTriangle* triangle2 = vtkTriangle::New();

triangle2->GetPointIds()->SetNumberOfIds(3);

triangle2->GetPointIds()->SetId(0, 3);

triangle2->GetPointIds()->SetId(1, 4);

triangle2->GetPointIds()->SetId(2, 5);

vtkCellArray* cells = vtkCellArray::New();

cells->InsertNextCell(triangle1);

cells->InsertNextCell(triangle2);

vtkPolyData* polyData = vtkPolyData::New();

polyData->SetPoints(pts);

polyData->SetStrips(cells);

polyData->BuildCells();

polyData->Update();

cells->Delete();

triangle1->Delete();

triangle2->Delete();

pts->Delete();

return polyData;

}

vtkActor* makeActor(vtkPolyData* polydata)

{

vtkPolyDataMapper* mapper = vtkPolyDataMapper::New();

mapper->SetInput(polydata);

vtkActor* actor = vtkActor::New();

actor->SetMapper(mapper);

mapper->Delete();

return actor;

}

double calcSignDistance(double p[],double proj[],double normal[]){

double dis = (p[0]-proj[0])*(p[0]-proj[0])+(p[1]-proj[1])*(p[1]-proj[1])+(p[2]-proj[2])*(p[2]-proj[2]);

if( abs(dis) < 0.0001 )return 0;

//向量proj->p

double vec[] = {

p[0] - proj[0],

p[1] - proj[1],

p[2] - proj[2]

};

double dotRes = vtkMath::Dot(vec, normal);

if( dotRes >= 0 )

return sqrt(dis);

else

return -sqrt(dis);

}

void calcTriangleNormal(double a[],double b[],double c[],double normal[])

{

double v1[3];

vtkMath::Subtract(b,a,v1);

double v2[3];

vtkMath::Subtract(c,b,v2);

vtkMath::Cross(v1,v2,normal);

vtkMath::Normalize(normal);

}

//计算一个三角形和一个平面的交点

int calcIntersectPoints(double p2a[],double p2b[],double p2c[],vtkPlane* plane,double x1[],double x2[]){

//将三角形的顶点投影到plane上

double p2aProj[3], p2bProj[3], p2cProj[3];

plane->ProjectPoint(p2a, p2aProj);

plane->ProjectPoint(p2b, p2bProj);

plane->ProjectPoint(p2c, p2cProj);

double normal[3];

plane->GetNormal(normal);

//计算有符号距离

double d1a = calcSignDistance(p2a, p2aProj, normal);

double d1b = calcSignDistance(p2b, p2bProj, normal);

double d1c = calcSignDistance(p2c, p2cProj, normal);

cout<<"符号距离："<<d1a<<","<<d1b<<","<<d1c<<endl;

if( (d1a > 0 && d1b > 0 && d1c > 0) ||

(d1a < 0&& d1b < 0 && d1c < 0) ) {

cout<<"说明没有相交 "<<endl;

return 0;

} else if(d1a == 0 && d1b*d1c > 0){//交于顶点p2a

x1[0] = p2a[0];

x1[1] = p2a[1];

x1[2] = p2a[2];

return 1;

} else if(d1b == 0 && d1a*d1c > 0){//交于顶点p2b

x1[0] = p2b[0];

x1[1] = p2b[1];

x1[2] = p2b[2];

return 1;

} else if(d1c == 0 && d1a*d1b > 0){//交于顶点p2c

x1[0] = p2c[0];

x1[1] = p2c[1];

x1[2] = p2c[2];

return 1;

} else {

//说明三角形所在的面相交，还需要判断三角形是否相交

//找出符号相异的点，来计算triangle和另外一个三角形定义的面plane的交点

//a : d1a 和 d1b,d1c 位于不同侧

double a[3],b[3],c[3];

a[0] = b[0] = c[0] = 0;

a[1] = b[1] = c[1] = 0;

a[2] = b[2] = c[2] = 0;

if( (d1a < 0 && d1b > 0 && d1c > 0) ||

(d1a > 0 && d1b < 0 && d1c < 0)){

a[0] = p2a[0];a[1] = p2a[1];a[2] = p2a[2];

b[0] = p2b[0];b[1] = p2b[1];b[2] = p2b[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

if( (d1b < 0 && d1a > 0 && d1c > 0) ||

(d1b > 0 && d1a < 0 && d1c < 0)){

a[0] = p2b[0];a[1] = p2b[1];a[2] = p2b[2];

b[0] = p2a[0];b[1] = p2a[1];b[2] = p2a[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

if( (d1c < 0 && d1b > 0 && d1a > 0) ||

(d1c > 0 && d1b < 0 && d1a < 0)){

a[0] = p2c[0];a[1] = p2c[1];a[2] = p2c[2];

b[0] = p2b[0];b[1] = p2b[1];b[2] = p2b[2];

c[0] = p2a[0];c[1] = p2a[1];c[2] = p2a[2];

}

//一个点到平面距离为0，其它2个到平面距离符号相异

if(d1a == 0 && d1b*d1c < 0){

a[0] = p2b[0];a[1] = p2b[1];a[2] = p2b[2];

b[0] = p2c[0];b[1] = p2c[1];b[2] = p2c[2];

c[0] = p2a[0];c[1] = p2a[1];c[2] = p2a[2];

}

if(d1b == 0 && d1a*d1c < 0){

a[0] = p2a[0];a[1] = p2a[1];a[2] = p2a[2];

b[0] = p2b[0];b[1] = p2b[1];b[2] = p2b[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

if(d1c == 0 && d1a*d1b < 0){

a[0] = p2a[0];a[1] = p2a[1];a[2] = p2a[2];

b[0] = p2b[0];b[1] = p2b[1];b[2] = p2b[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

//说明有2个点位于triangle1所在平面内部

if(d1a != 0 && d1b == 0 && d1c == 0){

a[0] = p2a[0];a[1] = p2a[1];a[2] = p2a[2];

b[0] = p2b[0];b[1] = p2b[1];b[2] = p2b[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

if(d1b != 0 && d1a == 0 && d1c == 0){

a[0] = p2b[0];a[1] = p2b[1];a[2] = p2b[2];

b[0] = p2a[0];b[1] = p2a[1];b[2] = p2a[2];

c[0] = p2c[0];c[1] = p2c[1];c[2] = p2c[2];

}

if(d1c != 0 && d1a == 0 && d1b == 0){

a[0] = p2c[0];a[1] = p2c[1];a[2] = p2c[2];

b[0] = p2a[0];b[1] = p2a[1];b[2] = p2a[2];

c[0] = p2b[0];c[1] = p2b[1];c[2] = p2b[2];

}

//求直线 a - b 和 a - c 与 plane的交点 e,f

double t;

plane->IntersectWithLine(a,b,t,x1);

plane->IntersectWithLine(a,c,t,x2);

return 2;

}

//never

return -1;

}

void build(){

vtkRenderer* renderer = vtkRenderer::New();

vtkRenderWindow* renderWin = vtkRenderWindow::New();

renderWin->AddRenderer(renderer);

vtkRenderWindowInteractor* inter = vtkRenderWindowInteractor::New();

vtkInteractorStyleTrackballCamera* style = vtkInteractorStyleTrackballCamera::New();

inter->SetInteractorStyle(style);

inter->SetRenderWindow(renderWin);

style->Delete();

renderWin->Delete();

renderer->Delete();

vtkPolyData* trianglePolyData = prepareData();

//显示三角形

vtkActor* triActor1 = makeActor(trianglePolyData);

renderer->AddActor(triActor1);

//读取三角形

vtkIdList* triList1 = vtkIdList::New();

trianglePolyData->GetCellPoints(0, triList1);

vtkIdList* triList2 = vtkIdList::New();

trianglePolyData->GetCellPoints(1, triList2);

//取得三角形triangle1的所有顶点

double p1a[3];

trianglePolyData->GetPoint(triList1->GetId(0),p1a);

double p1b[3];

trianglePolyData->GetPoint(triList1->GetId(1),p1b);

double p1c[3];

trianglePolyData->GetPoint(triList1->GetId(2),p1c);

//取得三角形triangle2的所有顶点

double p2a[3];

trianglePolyData->GetPoint(triList2->GetId(0),p2a);

double p2b[3];

trianglePolyData->GetPoint(triList2->GetId(1),p2b);

double p2c[3];

trianglePolyData->GetPoint(triList2->GetId(2),p2c);

double n1[3], n2[3];

calcTriangleNormal(p1a, p1b, p1c,n1);

calcTriangleNormal(p2a, p2b, p2c,n2);

cout<<n1[0]<<","<<n1[1]<<","<<n1[2]<<endl;

cout<<n2[0]<<","<<n2[1]<<","<<n2[2]<<endl;

//以顶点p1a为原点 n1为法向量来构造vtkPlane

vtkPlane* plane1 = vtkPlane::New();

plane1->SetOrigin(p1a);

plane1->SetNormal(n1);

vtkPlane* plane2 = vtkPlane::New();

plane2->SetOrigin(p2a);

plane2->SetNormal(n2);

double x1[3], x2[3], x3[3], x4[4];

x1[0] = x1[1] = x1[2] = 0;

x2[0] = x2[1] = x2[2] = 0;

x3[0] = x3[1] = x3[2] = 0;

x4[0] = x4[1] = x4[2] = 0;

int numOfInter1 = calcIntersectPoints(p2a,p2b,p2c,plane1,x1,x2);

int numOfInter2 = calcIntersectPoints(p1a,p1b,p1c,plane2,x3,x4);

cout <<"两个三角形的交点个数："<<numOfInter1 <<","<<numOfInter2<<endl;

if(numOfInter1 == 1 || numOfInter2 == 1){

cout <<"两个三角形相交于一个顶点"<<endl;

} else if(numOfInter1 == 2 && numOfInter2 == 2){//交于两个点

//计算交点组成的线段有没有重合，如果有重合，则说明相交

double cp1[3],cp2[3],t1,t2;

double dis = vtkLine::DistanceBetweenLineSegments(x1,x2,x3,x4,cp1,cp2,t1,t2);

cout<<"两个三角形交线的距离："<<dis<<endl;

if(dis < 0.0001)

cout<<"两个三角形相交"<<endl;

else

cout<<"两个三角形未相交"<<endl;

}else{

cout<<"没有相交"<<endl;

}

inter->Initialize();

inter->Start();

trianglePolyData->Delete();

triActor1->Delete();

triList1->Delete();

triList2->Delete();

plane1->Delete();

plane2->Delete();

inter->Delete();

}

int main(){

build();

return 0;

}

下图是其中一个测试情况，即两个三角形相交于一个顶点：



备注：最后在判断两个三角形的交线是否相交时，采用的是判断两个交线之间的距离，如果距离为0，则说明相交，隐约的觉得这种办法不够好，但一时也找不到更好的判断方法，后面继续完善！