[Неочевидные алгоритмы очевидных вещей] Алгоритм 2. Принадлежность точки треугольнику в пространстве

Серия постов [Неочевидные алгоритмы очевидных вещей] будет содержать алгоритмы действий, которые кажутся очевидными и простыми, но если задать себе вопрос «как это делается?», то ответ является далеко не очевидным. Разумеется, все эти алгоритмы можно найти в литературе. Под катом располагается алгоритм определения принадлежности точки P треугольнику ABC в пространстве.

Принадлежность точки P треугольнику ABC в пространстве

Идея

Иногда встречается необходимость определения принадлежности точки какой-нибудь геометрической фигуре. Самой простейшей является треугольник ABC.

Дано:

Точка P(P_x,P_y,P_z), треугольник с вершинами: A(A_x,A_y,A_z), B(B_x,B_y,B_z), C(C_x,C_y,C_z).

методика

Формируются три треугольника: ABP, ACP, BCP. После вычисляются их площади S_ABP,S_ACP,S_BCP. После этого сверяется сумма этих площадей с площадью треугольника S_ABC. Если точка лежит на треугольнике ABC, то треугольники ABP, ACP, BCP будут просто частями треугольника ABC, и сумма их площадей будет равна его площади S_ABC. Если же точка не принадлежит треугольнику, сумма площадей S_ABP,S_ACP,S_BCP превысит площадь треугольника ABC.
Легкое лирическое отступление для тех, кто не помнит, чему равна площадь треугольника: проще всего использовать формулу Герона, которая позволяет найти площадь треугольника, зная только его стороны, оч. удобно

Код функции:

#define EPS 1e-10

float triangle_square(float a,float b, float c){
   float p=(a+b+c)/2;
   return sqrt(p*(p-a)*(p-b)*(p-c));
}
float inside_triangle(float P_x,float P_y, float P_z, float A_x, float A_y, float A_z, float B_x, float B_y, float B_z, float C_x, float C_y, float C_z){
  int inside=0;
  float AB=sqrt( (A_x-B_x)*(A_x-B_x) + (A_y-B_y)*(A_y-B_y) + (A_z-B_z)*(A_z-B_z) );
  float BC=sqrt( (B_x-C_x)*(B_x-C_x) + (B_y-C_y)*(B_y-C_y) + (B_z-C_z)*(B_z-C_z) );
  float CA=sqrt( (A_x-C_x)*(A_x-C_x) + (A_y-C_y)*(A_y-C_y) + (A_z-C_z)*(A_z-C_z) );

  float AP=sqrt( (P_x-A_x)*(P_x-A_x) + (P_y-A_y)*(P_y-A_y) + (P_z-A_z)*(P_z-A_z) );
  float BP=sqrt( (P_x-B_x)*(P_x-B_x) + (P_y-B_y)*(P_y-B_y) + (P_z-B_z)*(P_z-B_z) );
  float CP=sqrt( (P_x-C_x)*(P_x-C_x) + (P_y-C_y)*(P_y-C_y) + (P_z-C_z)*(P_z-C_z) );
  float diff=(triangle_square(AP,BP,AB)+triangle_square(AP,CP,CA)+triangle_square(BP,CP,BC))-triangle_square(AB,BC,CA);
        if (fabs(diff)<EPS) inside=1;
  return inside;
}

Если вдруг кому лень прикреплять math.h, то можете воспользоваться функцией корня из этого поста.