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

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

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

假设这两个三角形为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++代码如下：

#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，则说明相交，隐约的觉得这种办法不够好，但一时也找不到更好的判断方法，后面继续完善！