射线与三角形求交的计算
2010-01-09 12:24
302 查看
将射线定义为一个基点和向量,那么射线上的任意点可以表示为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++):
当有大量射线和三角形求交的时候,计算是很耗时的。下面我们用TBB库来实现并行计算。我们使用TBB提供的parallel_for. 为了使用parallel_for,我们需要定义一个函数对象类来包装计算的代码:
在ParallelComputeIntersections2::operator()中,我们加入了一些代码来快速判断一条射线和一个三角形有没有可能相交。这些代码对于提高运算速度是很有用的。经过在我本地的数据集测试中(约2万个射线和2万个三角形),这个判断使原来18s的计算时间减少到15s.
下面是主函数
代码中使用了TBB提供的一个并行向量来存储运算结果。
使用上面提到的数据集测试,在我的双核笔记本上结果如下.
在没有加入快速判断代码的情况下:
顺序计算:38.25s
并行计算:18.57s
加入快速判断代码的情况下:
顺序计算:26.31s
并行计算:15.11s
默认TBB会根据当前CPU核的数量来设置并行任务数量。从结果来看,效果还是比较明显的。
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_PRECISION typedef double valuetype; #else typedef float valuetype; #endif #define EPSILON 10E-6 struct 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); }; //typedefs typedef 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 v2 v.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 triangle bool 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 triangle void 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 //Secondary class 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 //Master class 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_FOR parallel_for(blocked_range<size_t>(0, mVecRays.size()), ParallelComputeIntersections2(mVecTriangles[i], mVecRays, mIntersectPoints), auto_partitioner()); #else const 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 point void 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 file ifstream 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 file ifstream 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 file int 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 file int 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_TIME int main(int argc, char* argv[]) { #ifdef _SHOW_USED_TIME DWORD 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 computing VecPoint3D intersectPoints; Sequential_ComputeIntersections(vecRays, vecTriangles, intersectPoints); WriteIntersectPoints(outputFile, intersectPoints); #else //Parallel computing task_scheduler_init init; tbb::concurrent_vector<Point3D> vecIntersectPoints; Parallel_ComputeIntersections(vecRays, vecTriangles,vecIntersectPoints); WriteIntersectPoints(outputFile, vecIntersectPoints); #endif #ifdef _SHOW_USED_TIME DWORD dwEnd = GetTickCount(); printf("End time:%d/n", dwEnd); DWORD dwDuration = dwEnd - dwStart; printf("Time used(sec):%f", dwDuration * 1.0/1000.0); #endif return 0; }
代码中使用了TBB提供的一个并行向量来存储运算结果。
使用上面提到的数据集测试,在我的双核笔记本上结果如下.
在没有加入快速判断代码的情况下:
顺序计算:38.25s
并行计算:18.57s
加入快速判断代码的情况下:
顺序计算:26.31s
并行计算:15.11s
默认TBB会根据当前CPU核的数量来设置并行任务数量。从结果来看,效果还是比较明显的。
相关文章推荐
- 使用射线矢量对三角形图元求交 D3DXIntersect()函数说明
- 射线与三角形求交,并判断是否在三角形内的完整代码(带测试)
- 再次修订后的版本。。。。。。1.0(发布版,射线求交三角形)
- 使用射线矢量对三角形图元求交 D3DXIntersect()函数说明
- 探索射线与三角形求交的方法
- 射线与三角形求交
- 三角形跟射线求交
- 三角形与射线相交
- 分治 - 计算几何 - BZOJ2458,[BeiJing2011]最小三角形
- 计算三角形的面积,并判断三角形的类型?
- VB中函数的用法,计算三角形面积
- HDU 4667 Building Fence(2013多校7 1002题 计算几何,凸包,圆和三角形)
- 创建一个几何类型类,其中有计算面积getArea()和周长getPerimeter()抽象方法,然后通过它派生出三角形类、圆形类、矩形类,并通过测试类进行测试
- [计算几何]POJ2079 求点集中面积最大的三角形
- 计算边长为abc三角形的面积
- 判断空间射线是否穿过空间三角形的程序实现
- C++ 计算三角形面积
- 创建一个三角形类并且使用成员函数计算三角形的周长和面积《2》
- 3835. 计算三角形的周长
- 计算任意三边三角形的面积的算法