한 점의 삼각형포함문제

#include <stdio.h>

#include <math.h>

typedef struct _point {

    double x;

    double y;

    double z;

} point;

int inside_tri(point s, point a, point b, point c) {

    double as_x = s.x-a.x;

    double as_y = s.y-a.y;

    int s_ab = (b.x-a.x)*as_y-(b.y-a.y)*as_x > 0.0;

    if((c.x-a.x)*as_y-(c.y-a.y)*as_x > (double)(1- s_ab)) return 1;

    if((c.x-b.x)*(s.y-b.y)-(c.y-b.y)*(s.x-b.x) > (double)s_ab) return 1;

    return 0;

}

int main() {

    // point 1, 2, 3, 4

    point p1, p2, p3, p4;

    // vector 1, 2

    point v1, v2;

    // equation of the plane a*x+b*y+c*z=d 를 만족하는 a, b, c, d를 구한다.

    // a b c d and distance

    double a, b, c, d, distance;

    printf("point1 : "); // 점1을 읽는다.

    scanf("%lf %lf %lf", &p1.x, &p1.y, &p1.z);

    printf("point2 : "); // 점2를 읽는다.

    scanf("%lf %lf %lf", &p2.x, &p2.y, &p2.z);

    printf("point3 : "); // 점3을 읽는다

    scanf("%lf %lf %lf", &p3.x, &p3.y, &p3.z);

    printf("point4 : "); // 점4를 읽는다.

    scanf("%lf %lf %lf", &p4.x, &p4.y, &p4.z);

    // 삼각형의 3개의 점으로 부터 2개의 벡터를 구한다. v1과 v2

    v1.x=p3.x-p1.x;

    v1.y=p3.y-p1.y;

    v1.z=p3.z-p1.z;

    v2.x=p2.x-p1.x;

    v2.y=p2.y-p1.y;

    v2.z=p2.z-p1.z;

    // v1과 v2의 벡터곱으로 부터 a, b, c를 구한다.

    a=v1.y*v2.z-v1.z*v2.y;

    b=v1.z*v2.x-v1.x*v2.z;

    c=v1.x*v2.y-v1.y*v2.x;

    // a*x+b*y+c*z=d 에 p1을 대입하여 d를 구한다.

    d = a*p1.x + b*p1.y + c*p1.z;

    printf("a=%g b=%g c=%g d=%g\n", a, b, c, d);

    // p4가 삼각형을 이루는 평면상에 있는지 조사하기 위해서

    // a*x+b*y+c*z-d 식에 대입해서 0이 되면 평면상에 있는 것으로 판정한다.

    distance = a*p4.x + b*p4.y + c*p4.z - d;

    distance = distance / sqrt(a*a+b*b+c*c);

    // distance가 충분히 0에 가까우면 평면상에 있는 것으로 판정한다.

    if(fabs(distance) < 0.0000001) {

        printf("pont4 is on the plane : distance = %g\n", distance);

        // p4가 삼각형내에 있는지 판별한다.

        if(inside_tri(p4, p1, p2, p3)) {

            printf("point4 is inside the triangle\n");

        } else {

            printf("point4 is outside the triangle\n");

        }

    } else {

        printf("point4 is not on the plane : distance = %g\n", distance);

    }

    return 0;

}