#include "render.h"

#define EPSILON 1.0E-20

/* Clip.c -- intersection tests and closest approach tests for line
segments. */


/* compute intersection of a line segment with a plane.
   return values:
    1: both points on positive side of plane
    2: both points on negative side
    3: transition from positive to negative
    4: transition from negative to positive
 */

static int clip3d(PG_PLANE *a,PG_VECTOR *x0,PG_VECTOR *x1,PG_VECTOR *ret)
{
float t0,t1;
    t0= (x0->x) * (a->x) +
        (x0->y) * (a->y) +
        (x0->z) * (a->z) - (a->n);

    t1= (x1->x) * (a->x) +
        (x1->y) * (a->y) +
        (x1->z) * (a->z) - (a->n);

    if ((t0>=0.0F)&&(t1>=0.0F))
        return 1;
    if ((t0<=0.0F)&&(t1<=0.0F))
        return 2;

    /* intersection point is
        X=X0 + (X1-X0)*[ n - <(x,y,z),X0> ] / <(x,y,z),(X1-X0)>
       the denominator just happens to be t1-t0, and the
       [bracketed term] is -t0. */

/* this should never be true if t0*t1<0, but I'm paranoid about roundoff
    error. */

    if (t1-t0 == 0.0F) {
        return 1;
		}

    ret->x = (x0->x) - ((x1->x) - (x0->x))*t0/(t1-t0);
    ret->y = (x0->y) - ((x1->y) - (x0->y))*t0/(t1-t0);
    ret->z = (x0->z) - ((x1->z) - (x0->z))*t0/(t1-t0);

    if (t0>0.0F)
        return 3;
    return 4;
}

/* clippoly clips a polygon specified by the array pts to the intersection
    of half spaces bounded by plns.  The return is in out, and the
function returns the number of points.*/

#define MAXPOINTS 100
PG_VECTOR buf1[MAXPOINTS];
PG_VECTOR buf2[MAXPOINTS];

int ClipPoly(int npts,PG_VECTOR *pts,PG_VECTOR *out,int nplns,PG_PLANE *plns)
{
int i,j,s,os,fs;
int ns,nd;
PG_VECTOR *src;
PG_VECTOR *dest;
PG_VECTOR ret;

    src=buf2;
    dest=buf1;
    ns=npts;

    for (i=0;i<npts;i++) {
        src[i]=pts[i];
    }

    for (i=0;i<nplns;i++) {
        if (ns<2)
            return 0;
        nd=0;
        os=0;
        for (j=0;j<ns-1;j++) {
            s=clip3d(&(plns[i]),&(src[j]),&(src[j+1]),&ret);
            switch (s) {
              case 1:
                if ((os==2)||(os==3)) {
                    dest[nd]=src[j];
                    nd++;
                }
                dest[nd]=src[j+1];
                nd++;
                break;
              case 2:
                break;
              case 3:
                dest[nd]=ret;
                nd++;
                break;
              case 4:
                dest[nd]=ret;
                nd++;
                dest[nd]=src[j+1];
                nd++;
                break;
            }
            os=s;
            if (j==0)
                fs=s;
        }
        /* the closing line of the polygon is a special case. */
        s=clip3d(&(plns[i]),&(src[ns-1]),&(src[0]),&ret);
        switch (s) {
          case 1:
            if ((os==2)||(os==3)) {
                dest[nd]=src[ns-1];
                nd++;
            }
            dest[nd]=src[0];
            nd++;
            break;
          case 2:
            if (fs==1) {
                dest[nd]=src[0];
                nd++;
            }
            break;
          case 3:
            dest[nd]=ret;
            nd++;
            if (fs==1) {
                dest[nd]=src[0];
                nd++;
            }
            break;
          case 4:
            dest[nd]=ret;
            nd++;
            dest[nd]=src[0];
            nd++;
            break;
        }
        src=dest;
        ns=nd;
        if (src==buf1)
            dest=buf2;
        else
            dest=buf1;
    }
    if (ns<2)
        return 0;
    for(i=0;i<ns;i++) {
        out[i]=src[i];
    }
    return ns;
}

#if 0
/* this returns the square of the closest approach for two trajectories,
each being a straight line traveled at constant velocity.  The two
paths are
(from x0 at time t0 to x1 at time t1) and
(from y0 at time t0 to y1 at time t1).
note that t0 and t1 are not specified, because the answer doesn't depend
on them.  Also, note this is NOT the same as the distance between the
two line segments.
I'm returning the square of the distance, because it's usually faster
to square the number you're comparing to than to take a sqrt(). */

float sqrdist3d(PG_VECTOR *x0,PG_VECTOR *x1,PG_VECTOR *y0,PG_VECTOR *y1)
{
float t;
PG_VECTOR v;
PG_VECTOR d;

/* d is X0-Y0 */
    d.x=(x0->x)-(y0->x);
    d.y=(x0->y)-(y0->y);
    d.z=(x0->z)-(y0->z);

/* v is the relative velocity of the two points. */

    v.x=(x1->x) - (y1->x) - d.x;
    v.y=(x1->y) - (y1->y) - d.y;
    v.z=(x1->z) - (y1->z) - d.z;

    /* the time of closest approach is
        t= <(X0-Y0),V> / <V,V>
    */

    t= (v.x*v.x) + (v.y*v.y) + (v.z*v.z);

/* if |V|==0, it doesn't matter what time we choose, so pick 0. */
    if (t!=0.0F) {
        t=(d.x*v.x + d.y*v.y + d.z*v.z)/t;
    }

/* no I've got the time of closest approach.  Is it in the interval
    we're testing? */

    if (t<0.0F) {
        /* starting point is closest. */
        return (d.x*d.x) + (d.y*d.y) + (d.z*d.z);
    }

    if (t>1.0F) {
        /* endpoint is closest */
        return  (d.x+v.x)*(d.x+v.x) +
                (d.y+v.y)*(d.y+v.y) +
                (d.z+v.z)*(d.z+v.z);
    }
    return  (d.x+t*v.x)*(d.x+t*v.x) +
            (d.y+t*v.y)*(d.y+t*v.y) +
            (d.z+t*v.z)*(d.z+t*v.z);

}
#endif