/*已知直线（方向），和空间3个点构成的三角形，求直线与该三角形有效交点的算法。*/

 

//在Visual C++6.0及Turbo C3.0下编译通过
#include "stdio.h"
#include "conio.h"
#include "math.h"

#define MINERROR 0.0001
//定义空间点坐标结构
typedef struct SpacePoint
{
 float x;
 float y;
 float z;
}SpacePoint;

//定义空间向量结构
typedef struct SpaceVector
{
 float m;
 float n;
 float p;
}SpaceVector;

float Distance(SpacePoint p1, SpacePoint p2);
float Area(float a, float b, float c);
int ValidPoint(SpacePoint LinePoint, SpaceVector LineV, SpacePoint TrianglePoint1,
      SpacePoint TrianglePoint2, SpacePoint TrianglePoint3);

void main()
{
 //空间直线上的一个点
 SpacePoint LinePoint;
 //构成三角形的空间三个点
 SpacePoint TrianglePoint1, TrianglePoint2, TrianglePoint3;
 //空间直线的方向向量
 SpaceVector LineV;

 /*------------输入数据------------*/
 printf("Input the coordinate of the point in the line:\n");
 scanf("%f%f%f", &LinePoint.x, &LinePoint.y, &LinePoint.z);
 printf("Input the vector of the line:\n");
 scanf("%f%f%f", &LineV.m, &LineV.n, &LineV.p);
 printf("Input the the coordinate of three point to conpose a triangle:\n");
 printf("Point1:\n");
 scanf("%f%f%f", &TrianglePoint1.x, &TrianglePoint1.y, &TrianglePoint1.z);
 printf("Point2:\n");
 scanf("%f%f%f", &TrianglePoint2.x, &TrianglePoint2.y, &TrianglePoint2.z);
 printf("Point3:\n");
 scanf("%f%f%f", &TrianglePoint3.x, &TrianglePoint3.y, &TrianglePoint3.z);
 
 /*------调用ValidPoint求有效交点-----------*/
 ValidPoint(LinePoint, LineV, TrianglePoint1, TrianglePoint2, TrianglePoint3);
 //按任意建结束
 getch();
}

//计算p1到p2的距离的平方
float Distance(SpacePoint p1, SpacePoint p2)
{
 float dist;
 dist = ((p2.x-p1.x)*(p2.x-p1.x)
  + (p2.y-p1.y)*(p2.y-p1.y)
  + (p2.z-p1.z)*(p2.z-p1.z));
 return (float)sqrt(dist);
}

//利用海伦公式求变成为a,b,c的三角形的面积
float Area(float a, float b, float c)
{
 float s = (a+b+c)/2;
 return (float)sqrt(s*(s-a)*(s-b)*(s-c));
}

int ValidPoint(SpacePoint LinePoint, SpaceVector LineV,
      SpacePoint TrianglePoint1, SpacePoint TrianglePoint2,
      SpacePoint TrianglePoint3)
{
 //三角形所在平面的法向量
 SpaceVector TriangleV;
 //三角形的边方向向量
 SpaceVector VP12, VP13;
 //直线与平面的交点
 SpacePoint CrossPoint;
 //平面方程常数项
 float TriD;

 /*-------计算平面的法向量及常数项-------*/
 //point1->point2
 VP12.m = TrianglePoint2.x - TrianglePoint1.x;
 VP12.n = TrianglePoint2.y - TrianglePoint1.y;
 VP12.p = TrianglePoint2.z - TrianglePoint1.z;
 //point1->point3
 VP13.m = TrianglePoint3.x - TrianglePoint1.x;
 VP13.n = TrianglePoint3.y - TrianglePoint1.y;
 VP13.p = TrianglePoint3.z - TrianglePoint1.z;
 //VP12xVP13
 TriangleV.m = VP12.n*VP13.p - VP12.p*VP13.n;
 TriangleV.n = -(VP12.m*VP13.p - VP12.p*VP13.m);
 TriangleV.p= VP12.m*VP13.n - VP12.n*VP13.m;
 //计算常数项
 TriD = -(TriangleV.m*TrianglePoint1.x
   + TriangleV.n*TrianglePoint1.y
   + TriangleV.p*TrianglePoint1.z);

 /*-------求解直线与平面的交点坐标---------*/
 /* 思路：
  *     首先将直线方程转换为参数方程形式，然后代入平面方程，求得参数t，
  * 将t代入直线的参数方程即可求出交点坐标
 */
 float tempU, tempD;  //临时变量
 tempU = TriangleV.m*LinePoint.x + TriangleV.n*LinePoint.y
   + TriangleV.p*LinePoint.z + TriD;
 tempD = TriangleV.m*LineV.m + TriangleV.n*LineV.n + TriangleV.p*LineV.p;
 //直线与平面平行或在平面上
 if(tempD == 0.0)
 {
  printf("The line is parallel with the plane.\n");
  return 0;
 }

 //计算参数t
 float t = -tempU/tempD;
 //计算交点坐标
 CrossPoint.x = LineV.m*t + LinePoint.x;
 CrossPoint.y = LineV.n*t + LinePoint.y;
 CrossPoint.z = LineV.p*t + LinePoint.z;

 /*----------判断交点是否在三角形内部---------*/
 //计算三角形三条边的长度
 float d12 = Distance(TrianglePoint1, TrianglePoint2);
 float d13 = Distance(TrianglePoint1, TrianglePoint3);
 float d23 = Distance(TrianglePoint2, TrianglePoint3);
 //计算交点到三个顶点的长度
 float c1 = Distance(CrossPoint, TrianglePoint1);
 float c2 = Distance(CrossPoint, TrianglePoint2);
 float c3 = Distance(CrossPoint, TrianglePoint3);
 //求三角形及子三角形的面积
 float areaD = Area(d12, d13, d23);  //三角形面积
 float area1 = Area(c1, c2, d12);    //子三角形1
 float area2 = Area(c1, c3, d13);    //子三角形2
 float area3 = Area(c2, c3, d23);    //子三角形3
 //根据面积判断点是否在三角形内部
 if(fabs(area1+area2+area3-areaD) > MINERROR)
 {
  printf("There is no valid point of intersection\n");
  return 0;
 }
 printf("\nThe valid point of intersection:\n");
 printf("(%f, %f, %f)\n", CrossPoint.x, CrossPoint.y, CrossPoint.z); 
 return 1;
}

暂无评论