射线与三角形求交的计算

将射线定义为一个基点和向量，那么射线上的任意点可以表示为P = P0 + tV,其中P0为基点，V为向量，t满足 t>=0.采用重心坐标定义三角形，三角形上的任意点可以表示为：P = w V0 + u V1 + v V2, 其中u+v+w = 1. 由此，将射线方程代入上式，则：

P0 + tV = (1 - (u+v)) V0 + uV1+vV2.

利用克莱姆法则计算得出t,u,v. 交点位于三角形内的条件为: 0 <= u <= 1, 0 <= v <= 1, 0 <= u+v <= 1. 否则位于三角形外部的平面上。

该方法是Moller, Trumbore与1997年提出的。

下面我们用代码实现（c++）： 

#ifdef _HIGH_PRECISIONtypedef double valuetype;#elsetypedef float valuetype;#endif#define EPSILON 10E-6struct CVector {valuetype x,y,z;};struct Vector3D{valuetype x,y,z;Vector3D();Vector3D(valuetype vx, valuetype vy, valuetype vz);Vector3D(const Vector3D& vec);Vector3D operator* (valuetype t);};Vector3D operator*(valuetype t, const Vector3D& v);struct Point3D{valuetype x,y,z;Point3D();Point3D(valuetype vx, valuetype vy, valuetype vz);Point3D(const Point3D& point);Point3D& operator=(const Point3D& p);Point3D operator+(const Vector3D& v) const;Vector3D operator-(const Point3D& other) const;};struct Ray3D{Point3D origin;Vector3D direction;Ray3D();Ray3D(const Point3D& p, const Vector3D& v);};struct Triangle3D{Point3D v0, v1, v2;Triangle3D();Triangle3D(const Point3D& p0, const Point3D& p1, const Point3D& p2);};//typedefstypedef std::vector<Triangle3D> VecTriangle3D;typedef std::vector<Ray3D> VecRay3D;typedef std::vector<Point3D> VecPoint3D;void Cross(const Vector3D& v1, const Vector3D& v2, Vector3D& v){//v = v1 x v2v.x = v1.y*v2.z - v1.z*v2.y;v.y = -(v1.x*v2.z - v1.z*v2.x);v.z= v1.x*v2.y - v1.y*v2.x;}valuetype Dot(const Vector3D& v1, const Vector3D& v2){return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;}//Get intersect point of line and trianglebool RayTriangleIntersect(const Triangle3D& tri, const Ray3D& line, Point3D& intersection){//Does ray parallel with triangle?Vector3D e1, e2, p, s, q;valuetype t, u, v, tmp;e1 = tri.v1 - tri.v0;e2 = tri.v2 - tri.v0;Cross(line.direction, e2, p);tmp = Dot(p, e1);if(tmp > -EPSILON && tmp < EPSILON){return false;}tmp = 1.0/tmp;s = line.origin - tri.v0;u = tmp * Dot(s,p);if(u < 0.0 || u > 1.0)return false;Cross(s, e1, q);v = tmp * Dot(line.direction, q);if(v < 0.0 || v > 1.0)return false;valuetype uv = u+v;if( uv > 1.0)return false;t = tmp * Dot(e2, q);if(t < 0.0)return false;intersection = line.origin + t * line.direction;return true;}//Compute intersect points of ray and trianglevoid Sequential_ComputeIntersections(const VecRay3D& vecRays, const VecTriangle3D& vecTriangles, VecPoint3D& intersectPoints){Point3D point;for(size_t i = 0; i < vecTriangles.size(); i++){const Triangle3D& tri = vecTriangles[i];Vector3D v0v1 = tri.v1 - tri.v0;Vector3D v0v2 = tri.v2 - tri.v0;Vector3D normalTri;Cross(v0v1, v0v2, normalTri);valuetype tmp1, tmp2;for (size_t j = 0; j < vecRays.size(); j++){const Ray3D& ray = vecRays[j];Vector3D v0p0 = ray.origin - tri.v0;tmp1 = Dot(normalTri, v0p0);tmp2 = Dot(normalTri, ray.direction);if(tmp1 * tmp2 > 0.0)continue;if(RayTriangleIntersect(tri, ray, point)){intersectPoints.push_back(point);}}}}

 

当有大量射线和三角形求交的时候，计算是很耗时的。下面我们用TBB库来实现并行计算。我们使用TBB提供的parallel_for. 为了使用parallel_for,我们需要定义一个函数对象类来包装计算的代码：

 

//Parallel computation class//Secondaryclass ParallelComputeIntersections2{public:ParallelComputeIntersections2(const Triangle3D& tri, const VecRay3D& vecRays, concurrent_vector<Point3D>& vecIntersectPoints):mTri(tri), mVecRays(vecRays),mIntersectPoints(vecIntersectPoints){}void operator() ( const blocked_range<size_t>& r ) const{Point3D point;Vector3D v0v1 = mTri.v1 - mTri.v0;Vector3D v0v2 = mTri.v2 - mTri.v0;Vector3D normalTri;Cross(v0v1, v0v2, normalTri);valuetype tmp1, tmp2;for (size_t j = r.begin(); j != r.end(); ++j){const Ray3D& ray = mVecRays[j];Vector3D v0p0 = ray.origin - mTri.v0;tmp1 = Dot(normalTri, v0p0);tmp2 = Dot(normalTri, ray.direction);if(tmp1 * tmp2 > 0.0)continue;if(RayTriangleIntersect(mTri, ray, point)){mIntersectPoints.push_back(point);}}}private:const VecRay3D& mVecRays;const Triangle3D& mTri;concurrent_vector<Point3D>& mIntersectPoints;};#define _DOUBLE_PARALLEL_FOR//Masterclass ParallelComputeIntersections{public:ParallelComputeIntersections(const VecRay3D& vecRays, const VecTriangle3D& vecTriangles, concurrent_vector<Point3D>& vecIntersectPoints): mVecRays(vecRays), mVecTriangles(vecTriangles),mIntersectPoints(vecIntersectPoints){}void operator() ( const blocked_range<size_t>& r ) const{Point3D point;for(size_t i = r.begin(); i != r.end(); ++i){#ifdef _DOUBLE_PARALLEL_FORparallel_for(blocked_range<size_t>(0, mVecRays.size()),ParallelComputeIntersections2(mVecTriangles[i], mVecRays, mIntersectPoints), auto_partitioner());#elseconst Triangle3D& tri = mVecTriangles[i];Vector3D v0v1 = tri.v1 - tri.v0;Vector3D v0v2 = tri.v2 - tri.v0;Vector3D normalTri;Cross(v0v1, v0v2, normalTri);valuetype tmp1, tmp2;for (size_t j = 0; j < mVecTriangles.size(); j++){const Ray3D& ray = mVecRays[j];Vector3D v0p0 = ray.origin - tri.v0;tmp1 = Dot(normalTri, v0p0);tmp2 = Dot(normalTri, ray.direction);if(tmp1 * tmp2 > 0.0)continue;if(RayTriangleIntersect(tri, ray, point)){mIntersectPoints.push_back(point);}}#endif}}private:const VecRay3D& mVecRays;const VecTriangle3D& mVecTriangles;concurrent_vector<Point3D>& mIntersectPoints;};//Parallel compute entry pointvoid Parallel_ComputeIntersections(const VecRay3D& vecRays, const VecTriangle3D& vecTriangles, concurrent_vector<Point3D>& vecIntersectPoints){parallel_for(blocked_range<size_t>(0, vecTriangles.size()), ParallelComputeIntersections(vecRays, vecTriangles, vecIntersectPoints), auto_partitioner());} 

 

在ParallelComputeIntersections2::operator()中，我们加入了一些代码来快速判断一条射线和一个三角形有没有可能相交。这些代码对于提高运算速度是很有用的。经过在我本地的数据集测试中（约2万个射线和2万个三角形），这个判断使原来18s的计算时间减少到15s.

下面是主函数

void ReadTriangles(const char* szTriangleFile, VecTriangle3D& vecTriangles){// read the geometry input fileifstream streamGeom;streamGeom.open(szTriangleFile, ios_base::in);if(!streamGeom.is_open()){printf("Failed to open triangle file/n");return;}int numOfTriangles = 0;streamGeom >> numOfTriangles;vecTriangles.resize(numOfTriangles);//Triangle3D tri;for (int i = 0; i < numOfTriangles; i++){Triangle3D& tri = vecTriangles[i];streamGeom >> tri.v0.x >> tri.v0.y >> tri.v0.z;streamGeom >> tri.v1.x >> tri.v1.y >> tri.v1.z;streamGeom >> tri.v2.x >> tri.v2.y >> tri.v2.z;}streamGeom.close();}void ReadRays(const char* szRayFile, VecRay3D& vecRays){// read the ray input fileifstream streamRay;streamRay.open(szRayFile, ios_base::in);if(!streamRay.is_open()){printf("Failed to open ray file/n");return;}int numOfRays = 0;streamRay >> numOfRays;vecRays.resize(numOfRays);for (int i = 0; i < numOfRays; i++){Ray3D& ray = vecRays[i];streamRay >> ray.origin.x >> ray.origin.y >> ray.origin.z;streamRay >> ray.direction.x >> ray.direction.y >> ray.direction.z; }streamRay.close();}void WriteIntersectPoints(const char* szOutputFile, const VecPoint3D& intersectPoints){ofstream streamOut;streamOut.open(szOutputFile, ios_base::out);if(!streamOut.is_open()){printf("Failed to open output file/n");return;}// we now output the results into the specified fileint numOfIntersect = intersectPoints.size();streamOut << numOfIntersect << endl;for (int i = 0; i < numOfIntersect; i++){const Point3D& currentIntersect = intersectPoints[i];streamOut << currentIntersect.x << " " << currentIntersect.y << " " << currentIntersect.z << endl;}streamOut.close();}void WriteIntersectPoints(const char* szOutputFile, const concurrent_vector<Point3D>& intersectPoints){ofstream streamOut;streamOut.open(szOutputFile, ios_base::out);if(!streamOut.is_open()){printf("Failed to open output file/n");return;}// we now output the results into the specified fileint numOfIntersect = intersectPoints.size();streamOut << numOfIntersect << endl;for (int i = 0; i < numOfIntersect; i++){const Point3D& currentIntersect = intersectPoints[i];streamOut << currentIntersect.x << " " << currentIntersect.y << " " << currentIntersect.z << endl;}streamOut.close();}#define _SHOW_USED_TIMEint main(int argc, char* argv[]){#ifdef _SHOW_USED_TIMEDWORD dwStart = GetTickCount();printf("Start time: %d/n", dwStart);#endif if (argc != 4) { printf("Usage: InterSect.exe geometry_input.txt ray_input.txt output.txt/n"); return 0; } const char* geomFile = argv[1];const char* rayFile = argv[2];const char* outputFile = argv[3];VecTriangle3D vecTriangles;VecRay3D vecRays;ReadTriangles(geomFile, vecTriangles);ReadRays(rayFile, vecRays);#ifndef _PARALLEL_COMPUTE //Sequential computingVecPoint3D intersectPoints; Sequential_ComputeIntersections(vecRays, vecTriangles, intersectPoints);WriteIntersectPoints(outputFile, intersectPoints);#else//Parallel computingtask_scheduler_init init;tbb::concurrent_vector<Point3D> vecIntersectPoints;Parallel_ComputeIntersections(vecRays, vecTriangles,vecIntersectPoints); WriteIntersectPoints(outputFile, vecIntersectPoints);#endif#ifdef _SHOW_USED_TIMEDWORD dwEnd = GetTickCount();printf("End time:%d/n", dwEnd);DWORD dwDuration = dwEnd - dwStart;printf("Time used(sec):%f", dwDuration * 1.0/1000.0);#endifreturn 0;} 

 

代码中使用了TBB提供的一个并行向量来存储运算结果。

使用上面提到的数据集测试,在我的双核笔记本上结果如下.

在没有加入快速判断代码的情况下：

 

顺序计算：38.25s

并行计算：18.57s

 

加入快速判断代码的情况下：

顺序计算：26.31s

并行计算：15.11s

 

默认TBB会根据当前CPU核的数量来设置并行任务数量。从结果来看，效果还是比较明显的。