pl2plane.c

图形处理算法合集3：包括图形处理、帧缓存技术、渲染、矩阵运算、建模方法

/******************************************************************************
                       Plane-to-Plane Intersection
                       Author: Priamos Georgiades
******************************************************************************/
#include <stdio.h>
#include <math.h>

#define TRUE 1
#define FALSE 0
#define MRG_ZERO 1.0e-8

typedef struct vect3str {
    float x, y, z;
    } vect3;
typedef struct vect4str {
    float x, y, z, w;
    } vect4;

#define Vect3Init(a, b, c, v) { \
    (v).x = a; (v).y = b; (v).z = c; \
    }
#define Vect3Muls(s, v1, v2) { \
    (v2).x = s * (v1).x; \
    (v2).y = s * (v1).y; \
    (v2).z = s * (v1).z; \
    }
#define Vect3Cross(v1, v2, v3) { \
    (v3).x = (v1).y * (v2).z - (v1).z * (v2).y; \
    (v3).y = (v1).z * (v2).x - (v1).x * (v2).z; \
    (v3).z = (v1).x * (v2).y - (v1).y * (v2).x; \
    }

/*
Calculate the line of intersection between two planes. The two planes are
specified by their equations in the form
              P->x * X + P->y * Y + P->z * Z + P->w = 0.
Initialize the unit direction vector of the line of intersection in xdir.
Pick the point on the line of intersection on the coordinate plane most normal
to xdir. Return TRUE if successful, FALSE otherwise (indicating that the planes
don't intersect). The order in which the planes are given determines the choice
of direction of xdir.
*/
int GetXLine(vect4 *pl1, vect4 *pl2, vect3 *xdir, vect3 *xpt)
{
float invdet;  /* inverse of 2x2 matrix determinant */
vect3 dir2;    /* holds the squares of the coordinates of xdir */

Vect3Cross(*pl1, *pl2, *xdir)
dir2.x = xdir->x * xdir->x;
dir2.y = xdir->y * xdir->y;
dir2.z = xdir->z * xdir->z;


if (dir2.z > dir2.y && dir2.z > dir2.x && dir2.z > MRG_ZERO)
    {
    /* then get a point on the XY plane */
    invdet = 1.0 / xdir->z;
    /* solve < pl1.x * xpt.x + pl1.y * xpt.y = - pl1.w >
             < pl2.x * xpt.x + pl2.y * xpt.y = - pl2.w > */
    Vect3Init(pl1->y * pl2->w - pl2->y * pl1->w,
              pl2->x * pl1->w - pl1->x * pl2->w, 0.0, *xpt)
    }
else if (dir2.y > dir2.x && dir2.y > MRG_ZERO)
    {
    /* then get a point on the XZ plane */
    invdet = -1.0 / xdir->y;	/*** correction ***/
    /* solve < pl1.x * xpt.x + pl1.z * xpt.z = -pl1.w >
             < pl2.x * xpt.x + pl2.z * xpt.z = -pl2.w > */
    Vect3Init(pl1->z * pl2->w - pl2->z * pl1->w, 0.0,
              pl2->x * pl1->w - pl1->x * pl2->w, *xpt)
    }
else if (dir2.x > MRG_ZERO)
    {
    /* then get a point on the YZ plane */
    invdet = 1.0 / xdir->x;
    /* solve < pl1.y * xpt.y + pl1.z * xpt.z = - pl1.w >
             < pl2.y * xpt.y + pl2.z * xpt.z = - pl2.w > */
    Vect3Init(0.0, pl1->z * pl2->w - pl2->z * pl1->w,
              pl2->y * pl1->w - pl1->y * pl2->w, *xpt)
    }
else /* xdir is zero, then no point of intersection exists */
    return FALSE;
Vect3Muls(invdet, *xpt, *xpt)
invdet = 1.0 / (float)sqrt(dir2.x + dir2.y + dir2.z);
Vect3Muls(invdet, *xdir, *xdir)
return TRUE;
}