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