// Shave and a Haircut
// (c) 2019 Epic Games
// US Patent 6720962

// ray poly adapted from Graphics Gems I
#define PVX1 37
#define PVX 36
static int      vox_been_done = 0;
static float    cx1 = 10000, cx2 = -10000, cy1 = 10000, cy2 = -10000, cz1 = 10000, cz2 = -10000;

//typedef struct {
//  int totalp;
//  int *plist;
//} POLYVOX;
static POLYVOX  pvox[PVX1 + 2][PVX1 + 2][PVX1 + 2];

//POLYVOX ***pvox=NULL;

static int      intersect_triangle( double orig[3], double dir[3], double vert0[3], double vert1[3], double vert2[3], double *t, double *u, double *v );
static float    ray_poly( VERT O, VERT ray, VERT * pp, VERT * nn, int i1, int i2, double *uu, double *vv );
static void     collision_setup( void );
static void     collision_destroy( void );
static float    ray_trace( VERT * O, VERT * ray, WFTYPE * csurf, double *uu, double *vv );
static int      lastDtotalfaces = 0;
static float   *Dface_equation = NULL;
static void
Dcalc_plane_eq( void )			// norm is a unit vector
{
	int             x;

	if( Gdrag_hair == 0 )
		if( COLLISION_METHOD == 1 )
#ifndef LW3D
//      if (skull_sphere>=0)
#endif
			if( DOCOLLIDE == 1 )
			{
				if( Dtotalfaces > 0 )
//if (lastDtotalfaces!=Dtotalfaces)
				{
					if( Dface_equation != NULL )
						Nfree( Dface_equation );
					Dface_equation = NULL;
					Dface_equation = ( float * ) malloc( Dtotalfaces * sizeof( float ) );
//printf ("malloc face equation\n");
					lastDtotalfaces = Dtotalfaces;
				}
#ifndef LW3D
//      if (skull_sphere>=0)
#endif
				if( COLLISION_METHOD == 1 )
					if( DOCOLLIDE == 1 )
						if( Dtotalfaces > 0 )
							for( x = 0; x < Dtotalfaces; x++ )
							{
								VERT            a, b;

								b = Dface_norm[x];
								a = Dv[Dfacelist[Dface_start[x]]];
								Dface_equation[x] = -( a.x * b.x + a.y * b.y + a.z * b.z );
							}
			}
}
static void
init_poly_vox( void )
{
	int             x, y, z;
	int             i, j;

//  if (COLLISION_METHOD==1)
//  if (DOCOLLIDE)
#ifdef crap
	if( pvox == NULL )
	{
		pvox = ( POLYVOX *** ) malloc( sizeof( POLYVOX ** ) * ( PVX1 + 2 ) );
		for( i = 0; i < PVX1 + 2; i++ )
		{
			pvox[i] = ( POLYVOX ** ) malloc( sizeof( POLYVOX * ) * ( PVX1 + 2 ) );

			for( j = 0; j < PVX1 + 2; j++ )
				pvox[i][j] = ( POLYVOX * ) malloc( sizeof( POLYVOX ) * ( PVX1 + 2 ) );
		}
	}
#endif
	if( Gdrag_hair == 0 )
		if( vox_been_done == 0 )
			for( x = 0; x < PVX1; x++ )
				for( y = 0; y < PVX1; y++ )
					for( z = 0; z < PVX1; z++ )
					{
						pvox[x][y][z].totalp = 0;
						pvox[x][y][z].plist = NULL;
					}
	vox_been_done = 1;
}
static void freevox( void )
{
	int             x, y, z;
	int             i, j;

//  if (COLLISION_METHOD==1)
//  if (DOCOLLIDE)
//      if (skull_sphere>0)

	if( Gdrag_hair == 0 )
		if( vox_been_done == 1 )
			for( x = 0; x < PVX1; x++ )
				for( y = 0; y < PVX1; y++ )
					for( z = 0; z < PVX1; z++ )
					{
						if( pvox[x][y][z].totalp > 0 )
							if( pvox[x][y][z].plist != NULL )
								free( pvox[x][y][z].plist );
						pvox[x][y][z].plist = NULL;
						pvox[x][y][z].totalp = 0;
//              if (pvox[x][y][z].totalp>0)
					}



}


static void initvox( void )
{
	int             x, y, z;
	int             i, j;


			for( x = 0; x < PVX1; x++ )
				for( y = 0; y < PVX1; y++ )
					for( z = 0; z < PVX1; z++ )
					{
						pvox[x][y][z].plist = NULL;
						pvox[x][y][z].totalp = 0;
					}



}


static void
setup_poly_vox( void )
{
	int             x, y, z, q, yy;
	int             doit = 1;

#ifndef NOLIB
	if( Gdynamics_mode == 2 )
		doit = 0;
//if (Gdynamics_mode==0) doit=0;
#endif


	if( Gdrag_hair == 0 )
		if( COLLISION_METHOD == 1 )
			doit = 1;

	if( doit )
		if( DOCOLLIDE )
		{
			if( vox_been_done == 0 )
				init_poly_vox(  );
			if( Dtotalverts > 0 )
			{
				for( x = 0; x < PVX1; x++ )
					for( y = 0; y < PVX1; y++ )
						for( z = 0; z < PVX1; z++ )
						{
							pvox[x][y][z].totalp = 0;
						}
				cx1 = 1000000;
				cx2 = -1000000;
				cy1 = 1000000;
				cy2 = -1000000;
				cz1 = 1000000;
				cz2 = -1000000;
				for( q = 0; q < Dtotalverts; q++ )
				{
					VERT           *a = NULL;

					a = &Dv[q];
					if( a->x < cx1 )
						cx1 = a->x;
					if( a->x > cx2 )
						cx2 = a->x;
					if( a->y < cy1 )
						cy1 = a->y;
					if( a->y > cy2 )
						cy2 = a->y;
					if( a->z < cz1 )
						cz1 = a->z;
					if( a->z > cz2 )
						cz2 = a->z;
				}
// tag the voxels
				for( yy = 0; yy < Dtotalfaces; yy++ )
				{
					float           bx1 = 1000000, bx2 = -1000000, by1 = 1000000, by2 = -1000000, bz1 = 1000000, bz2 = -1000000;

#ifndef LW3D
					if( Dmtlgroup[yy] == skull_sphere )
#endif
						for( q = Dface_start[yy]; q < Dface_end[yy]; q++ )
						{
							int             g;
							VERT           *a = NULL;

							g = Dfacelist[q];
							a = &Dv[g];
							if( a->x < bx1 )
								bx1 = a->x;
							if( a->x > bx2 )
								bx2 = a->x;
							if( a->y < by1 )
								by1 = a->y;
							if( a->y > by2 )
								by2 = a->y;
							if( a->z < bz1 )
								bz1 = a->z;
							if( a->z > bz2 )
								bz2 = a->z;
						}		// computed bounding box for poly
					bx1 -= cx1;
					bx1 /= ( cx2 - cx1 );
					bx1 *= PVX;
					bx2 -= cx1;
					bx2 /= ( cx2 - cx1 );
					bx2 *= PVX;
					by1 -= cy1;
					by1 /= ( cy2 - cy1 );
					by1 *= PVX;
					by2 -= cy1;
					by2 /= ( cy2 - cy1 );
					by2 *= PVX;
					bz1 -= cz1;
					bz1 /= ( cz2 - cz1 );
					bz1 *= PVX;
					bz2 -= cz1;
					bz2 /= ( cz2 - cz1 );
					bz2 *= PVX;
					bx1 = floor( bx1 );	// - 1;
					bx2 = ceil( bx2 );	// + 1;
					by1 = floor( by1 );	// - 1;
					by2 = ceil( by2 );	// + 1;
					bz1 = floor( bz1 );	// - 1;
					bz2 = ceil( bz2 );	// + 1;
					if( bx1 < bx2 )
						if( by1 < by2 )
							if( bz1 < bz2 )
								for( x = ( int ) bx1; x <= ( int ) bx2; x++ )
									for( y = ( int ) by1; y <= ( int ) by2; y++ )
										for( z = ( int ) bz1; z <= ( int ) bz2; z++ )
											if( ( x >= 0 ) && ( y >= 0 ) && ( z >= 0 ) && ( x < PVX ) && ( y < PVX ) && ( z < PVX ) )
												pvox[x][y][z].totalp++;
				}
// allocate the voxels
				for( x = 0; x < PVX1; x++ )
					for( y = 0; y < PVX1; y++ )
						for( z = 0; z < PVX1; z++ )
						{
							if( pvox[x][y][z].plist != NULL )
								free( pvox[x][y][z].plist );
							pvox[x][y][z].plist = NULL;
							if( pvox[x][y][z].totalp > 0 )
								pvox[x][y][z].plist = ( int * ) malloc( pvox[x][y][z].totalp * sizeof( int ) );
							pvox[x][y][z].totalp = 0;
						}
// tag em again and fill in the lists
				for( yy = 0; yy < Dtotalfaces; yy++ )
				{
					float           bx1 = 1000000, bx2 = -1000000, by1 = 1000000, by2 = -1000000, bz1 = 1000000, bz2 = -1000000;

#ifndef LW3D
					if( Dmtlgroup[yy] == skull_sphere )
#endif
						for( q = Dface_start[yy]; q < Dface_end[yy]; q++ )
						{
							int             g;
							VERT           *a = NULL;

							g = Dfacelist[q];
							a = &Dv[g];
							if( a->x < bx1 )
								bx1 = a->x;
							if( a->x > bx2 )
								bx2 = a->x;
							if( a->y < by1 )
								by1 = a->y;
							if( a->y > by2 )
								by2 = a->y;
							if( a->z < bz1 )
								bz1 = a->z;
							if( a->z > bz2 )
								bz2 = a->z;
						}		// computed bounding box for poly
					bx1 -= cx1;
					bx1 /= ( cx2 - cx1 );
					bx1 *= PVX;
					bx2 -= cx1;
					bx2 /= ( cx2 - cx1 );
					bx2 *= PVX;
					by1 -= cy1;
					by1 /= ( cy2 - cy1 );
					by1 *= PVX;
					by2 -= cy1;
					by2 /= ( cy2 - cy1 );
					by2 *= PVX;
					bz1 -= cz1;
					bz1 /= ( cz2 - cz1 );
					bz1 *= PVX;
					bz2 -= cz1;
					bz2 /= ( cz2 - cz1 );
					bz2 *= PVX;
					bx1 = floor( bx1 );	// - 1;
					bx2 = ceil( bx2 );	// + 1;
					by1 = floor( by1 );	// - 1;
					by2 = ceil( by2 );	// + 1;
					bz1 = floor( bz1 );	// - 1;
					bz2 = ceil( bz2 );	// + 1;
					if( bx1 < bx2 )
						if( by1 < by2 )
							if( bz1 < bz2 )
								for( x = ( int ) bx1; x <= ( int ) bx2; x++ )
									for( y = ( int ) by1; y <= ( int ) by2; y++ )
										for( z = ( int ) bz1; z <= ( int ) bz2; z++ )
											if( ( x >= 0 ) && ( y >= 0 ) && ( z >= 0 ) && ( x < PVX ) && ( y < PVX ) && ( z < PVX ) )
											{
												pvox[x][y][z].plist[pvox[x][y][z].totalp] = yy;
												pvox[x][y][z].totalp++;
											}
				}
			}
		}
}
static void     surface_collide2( int hh, float radius );
static void
surface_collide( int hh, float radius )
{
	surface_collide2( hh, radius );
//surface_collide2(  hh, radius );

}
static void
surface_collide2( int hh, float radius )
{
	VERT            pp[3];
	VERT            nn[3];
	VERT            avgnorm;
	VERT            vec;
	int             doit = 1;
	int             ret = 0;
	int             collide = 0;
	int             q;
	int             x, y, xxx, yyy, zzz, yy;
	int             slg = 0;

//printf("testing surface_collide here\n");

				doit = 1;
#ifndef NOLIB
	if( Gdynamics_mode == 2 )
		doit = 0;
	if( Gdynamics_mode == 0 )
		doit = 0;
	if( Gsculpt_cursors )
		doit = 1;
#endif

//              printf ("doit %d  method %d docllide %d\n",doit,COLLISION_METHOD,DOCOLLIDE);

	if( doit )
		if( COLLISION_METHOD == 1 )
#ifndef LW3D
//      if (skull_sphere>=0)
#endif
			if( DOCOLLIDE )
			{
				vec.x = 0;
				vec.y = 0;
				vec.z = 0;
//  printf ("looking for collision\n");
				if( Dtotalverts > 0 )
								if( hair[hh].restlength > 0 )
//for (q=0;q<1;q++)
								{
									int             gotone = 0;

									slg = 0;
									if( head >= 0 )
										if( hair[hh].mtl == head )
										{
											slg = 0;
										}
									if( beard >= 0 )
										if( hair[hh].mtl == beard )
										{
											slg = 1;
										}
									if( eyebrow >= 0 )
										if( hair[hh].mtl == eyebrow )
										{
											slg = 2;
										}
									if( eyelash >= 0 )
										if( hair[hh].mtl == eyelash )
										{
											slg = 3;
										}
									if( splines >= 0 )
										if( hair[hh].mtl == splines )
										{
											slg = 4;
										}
									if( Gsurface_collide[slg] )
										for( x = 1; x < 15; x++ )
											if( hair[hh].pfID[x] == -1 )
											{
//  x=vrt;
												{

													VERT           *a = NULL;
													VERT            aa;
													float           smallest = 99999.0f;
													float           smallestd = 99999.0f;

													vec.x = 0;
													vec.y = 0;
													vec.z = 0;
													a = &hair[hh].hv[x];
													aa = hair[hh].hv[x];
													{
														float           tx, ty, tz;

														tx = hair[hh].hv[x].x;
														ty = hair[hh].hv[x].y;
														tz = hair[hh].hv[x].z;
														tx -= cx1;
														tx /= ( cx2 - cx1 );
														tx *= PVX;
														ty -= cy1;
														ty /= ( cy2 - cy1 );
														ty *= PVX;
														tz -= cz1;
														tz /= ( cz2 - cz1 );
														tz *= PVX;
														xxx = ( int ) floor( tx );
														yyy = ( int ) floor( ty );
														zzz = ( int ) floor( tz );
													}

													if( ( xxx >= 0 ) && ( yyy >= 0 ) && ( zzz >= 0 ) && ( xxx < PVX ) && ( yyy < PVX ) && ( zzz < PVX ) )
													{
														if( pvox[xxx][yyy][zzz].totalp > 0 )
															for( yy = 0; yy < pvox[xxx][yyy][zzz].totalp; yy++ )
															{
																int             q;
																float           d;
																float           fd;

																y = pvox[xxx][yyy][zzz].plist[yy];
//          y=0;
																// plane test;
																d = VDot( Dface_norm[y], aa ) + Dface_equation[y];
																fd = d;
																fd = ( float ) fabs( d );
//          if (0==1)
//          if ((fd> 0.0001)&&(fd< hair[hh].restlength*(float)x)&&(fd>hair[hh].restlength/20.0)) // it's behind the plane but closer than the rad
																if( ( fabs( fd ) > 0.000001 ) && ( fd > hair[hh].restlength / ( float ) 5.0 ) );	//&&(fd>hair[hh].restlength/120.0)) // it's behind the plane but closer than the rad
																if( Dface_end[y] - Dface_start[y] > 2 )
																{
																	int             g;
																	float           dr;
																	int             i0, i1, i2;
																	double          t, uu, vv;
																	double          vert0[3], vert1[3], vert2[3];
																	double          o[3], vdir[3];

																	g = Dfacelist[Dface_start[y]];
																	vert2[0] = ( double ) Dv[g].x;
																	vert2[1] = ( double ) Dv[g].y;
																	vert2[2] = ( double ) Dv[g].z;
																	g = Dfacelist[Dface_start[y] + 1];
																	vert1[0] = ( double ) Dv[g].x;
																	vert1[1] = ( double ) Dv[g].y;
																	vert1[2] = ( double ) Dv[g].z;
																	g = Dfacelist[Dface_start[y] + 2];
																	vert0[0] = ( double ) Dv[g].x;
																	vert0[1] = ( double ) Dv[g].y;
																	vert0[2] = ( double ) Dv[g].z;
																	o[0] = ( double ) aa.x;
																	o[1] = ( double ) aa.y;
																	o[2] = ( double ) aa.z;
																	vdir[0] = ( double ) Dface_norm[y].x;
																	vdir[1] = ( double ) Dface_norm[y].y;
																	vdir[2] = ( double ) Dface_norm[y].z;
																	if( intersect_triangle( o, vdir, vert0, vert1, vert2, &t, &uu, &vv ) == 1 )
																	{
//                  t= -t;

																		if( ( float ) fabs( d ) < smallest )
																		{
																			//          t-=.00005;
//                      t*=.999999;
																			if( test_stat == 0 )
																				t += ( double ) ( ( hair[hh].restlength * ( float ) 1 ) / 6.0 );
																			if( test_stat == 1 )
																				t += ( double ) ( ( hair[hh].restlength * ( float ) 1 ) / 6.0 );
																			smallest = ( float ) fabs( d );
																			smallestd = ( float ) ( d );
																			vec.x = ( float ) ( ( t ) * vdir[0] );
																			vec.y = ( float ) ( ( t ) * vdir[1] );
																			vec.z = ( float ) ( ( t ) * vdir[2] );
																			collide = 1;
                  hair[hh].velocity[x].x=0;
                  hair[hh].velocity[x].y=0;
                  hair[hh].velocity[x].z=0;
//                  printf("colliding here\n");

																		}
//                  ret=1;                  
																	}
																}
															}
													}

//if (smallest> hair[hh].restlength*x*.99999) // hair is straight
//{
//
//}
//recalc_hair(hh);
//if (0==1)
//      if (a!=NULL)


													{
														float           dd;
														VERT            df;

														df.x = hair[hh].hv[x - 1].x - hair[hh].hv[x].x;
														df.y = hair[hh].hv[x - 1].y - hair[hh].hv[x].y;
														df.z = hair[hh].hv[x - 1].z - hair[hh].hv[x].z;
														df = Vnorm( df );
														dd = sqrt( vec.x * vec.x + vec.y * vec.y + vec.z * vec.z );
														if( dd < hair[hh].restlength * 15.0f )	// don't correct for distances over a couple links in length
															if( VDot( df, vec ) )	// is the correction vector in the same hemisphere as the vector to the last vert?

																if( Gsurface_collide[slg] == 1 )
																	if( smallestd < ( float ) ( ( hair[hh].restlength * ( float ) 1 ) / 12.0 ) )
																		if( hair[hh].pfID[x] == -1 )
																		{
																			hair[hh].hv[x].x += vec.x * .83;
																			hair[hh].hv[x].y += vec.y * .83;
																			hair[hh].hv[x].z += vec.z * .83;
																			if( test_stat != 1 )
																			{
																				hair[hh].hv[x].x += vec.x * .3;
																				hair[hh].hv[x].y += vec.y * .3;
																				hair[hh].hv[x].z += vec.z * .3;
																			}
																			// correct the length
																			df.x = hair[hh].hv[x].x - hair[hh].hv[x - 1].x;
																			df.y = hair[hh].hv[x].y - hair[hh].hv[x - 1].y;
																			df.z = hair[hh].hv[x].z - hair[hh].hv[x - 1].z;
																			df = Vnorm( df );
																			hair[hh].hv[x].x = hair[hh].hv[x - 1].x + df.x * hair[hh].restlength;
																			hair[hh].hv[x].y = hair[hh].hv[x - 1].y + df.y * hair[hh].restlength;
																			hair[hh].hv[x].z = hair[hh].hv[x - 1].z + df.z * hair[hh].restlength;

																			Gcollision_hit = 1;
																			gotone = 1;
hair[hh].velocity[x].x*=.4; // drag
hair[hh].velocity[x].y*=.4;
hair[hh].velocity[x].z*=.4;
//      if (resthair[hh].pfID[x]== -1)
//          spare[hh].hv[x]=hair[hh].hv[x];
																			hair[hh].velocity[x].x = vec.x * .97;
																			hair[hh].velocity[x].y = vec.y * .97;
																			hair[hh].velocity[x].z = vec.z * .97;
																		}
													}

//if (0==1)     
													//  for(x=0;x<15;x++)
//if (collide==1)
//  spare[hh].hv[x]=hair[hh].hv[x];
												}
											}	// surf collide

								}	//q
			}
}


static float
ray_trace( VERT * O, VERT * ray, WFTYPE * csurf, double *uu, double *vv )
{
	int             x;
	VERT            oo, ro;
	float           mint = 1000000.0f;
	float           t = 1000000.0f;

//float t=1000000.0f;
	VERT            pp[3], nn[3];
	int             i1 = 1, i2 = 2;
	float           d;

	oo.x = O->x;
	oo.y = O->y;
	oo.z = O->z;
	ro.x = ray->x;
	ro.y = ray->y;
	ro.z = ray->z;
	d = sqrt( ro.x * ro.x + ro.y * ro.y + ro.z * ro.z );
	if( d != 0 )
	{
		ro.x /= d;
		ro.y /= d;
		ro.z /= d;
	}
	for( x = 0; x < csurf->totalfaces; x++ )
	{
		int             g;
		float           u, v;
		VERT            avgnorm;
		int             i0 = 0;

		g = csurf->facelist[csurf->face_start[x]];
		pp[0] = csurf->v[g];
//  nn[0]=norms->v[g];
		nn[0] = csurf->vn[g];
		g = csurf->facelist[csurf->face_start[x] + 1];
		pp[1] = csurf->v[g];
		nn[1] = csurf->vn[g];
//  csurf->
		g = csurf->facelist[csurf->face_start[x] + 2];
		pp[2] = csurf->v[g];
//  nn[2]=norms->v[g];
		nn[2] = csurf->vn[g];
		avgnorm.x = ( nn[0].x + nn[1].x + nn[2].x ) / 3.0f;
		avgnorm.y = ( nn[0].y + nn[1].y + nn[2].y ) / 3.0f;
		avgnorm.z = ( nn[0].z + nn[1].z + nn[2].z ) / 3.0f;
		avgnorm = Vnorm( avgnorm );
		avgnorm.x = fabs( avgnorm.x );
		avgnorm.y = fabs( avgnorm.y );
		avgnorm.z = fabs( avgnorm.z );
		if( ( avgnorm.x >= avgnorm.y ) && ( avgnorm.x >= avgnorm.z ) )
			i0 = 0;
		if( ( avgnorm.y >= avgnorm.x ) && ( avgnorm.y >= avgnorm.z ) )
			i0 = 1;
		if( ( avgnorm.z >= avgnorm.x ) && ( avgnorm.z >= avgnorm.y ) )
			i0 = 2;
		if( i0 == 0 )
		{
			i1 = 1;
			i2 = 2;
		}
		if( i0 == 1 )
		{
			i1 = 0;
			i2 = 2;
		}
		if( i0 == 2 )
		{
			i1 = 0;
			i2 = 1;
		}
		t = ray_poly( oo, ro, pp, nn, i1, i2, uu, vv );
		if( t < d )
			if( t > 0 )
				if( t < mint )
				{
//      char tmp[255];
					mint = t;
				}
	}
	if( mint >= 1000.0 )
		mint = -1;
	if( mint > 0 )
	{
		oo.x = oo.x + ro.x * mint;
		oo.y = oo.y + ro.y * mint;
		oo.z = oo.z + ro.z * mint;
	}
	return ( mint );
}



static float
ray_poly( VERT O, VERT ray, VERT * pp, VERT * nn, int i1, int i2, double *uu, double *vv )
{
	float           t;
	int             doit;
	double          orig[3];
	double          dir[3];
	double          vert0[3];
	double          vert1[3];
	double          vert2[3];
	double          tt;

	orig[0] = ( double ) O.x;
	orig[1] = ( double ) O.y;
	orig[2] = ( double ) O.z;
	dir[0] = ( double ) ray.x;
	dir[1] = ( double ) ray.y;
	dir[2] = ( double ) ray.z;
	vert0[0] = ( double ) pp[0].x;
	vert0[1] = ( double ) pp[0].y;
	vert0[2] = ( double ) pp[0].z;
	vert1[0] = ( double ) pp[1].x;
	vert1[1] = ( double ) pp[1].y;
	vert1[2] = ( double ) pp[1].z;
	vert2[0] = ( double ) pp[2].x;
	vert2[1] = ( double ) pp[2].y;
	vert2[2] = ( double ) pp[2].z;
	doit = intersect_triangle( orig, dir, vert0, vert1, vert2, &tt, uu, vv );
	if( !doit )
		t = -999;
	else
		t = ( float ) tt;

	return ( t );
}


//#define TEST_CULL 1
#define EPSILON 0.000001
#define CROSS(dest,v1,v2) \
          dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
          dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
          dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
#define SUB(dest,v1,v2) \
          dest[0]=v1[0]-v2[0]; \
          dest[1]=v1[1]-v2[1]; \
          dest[2]=v1[2]-v2[2];
int static
intersect_triangle( double orig[3], double dir[3], double vert0[3], double vert1[3], double vert2[3], double *t, double *uu, double *vv )
{
	double          edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
	double          det, inv_det;

	/* find vectors for two edges sharing vert0 */
	SUB( edge1, vert1, vert0 );
	SUB( edge2, vert2, vert0 );
	/* begin calculating determinant - also used to calculate U parameter */
	CROSS( pvec, dir, edge2 );
	/* if determinant is near zero, ray lies in plane of triangle */
	det = DOT( edge1, pvec );
#ifdef TEST_CULL				/* define TEST_CULL if culling is desired */
	if( det < EPSILON )
		return 0;
	/* calculate distance from vert0 to ray origin */
	SUB( tvec, orig, vert0 );
	/* calculate U parameter and test bounds */
	*uu = DOT( tvec, pvec );
	if( *uu < 0.0 || *uu > det )
		return 0;
	/* prepare to test V parameter */
	CROSS( qvec, tvec, edge1 );
	/* calculate V parameter and test bounds */
	*vv = DOT( dir, qvec );
	if( *vv < 0.0 || *uu + *vv > det )
		return 0;
	/* calculate t, scale parameters, ray intersects triangle */
	*t = DOT( edge2, qvec );
	inv_det = 1.0 / det;
	*t *= inv_det;
	*uu *= inv_det;
	*vv *= inv_det;
#else /* the non-culling branch */
	if( det > -EPSILON && det < EPSILON )
		return 0;
	inv_det = 1.0 / det;
	/* calculate distance from vert0 to ray origin */
	SUB( tvec, orig, vert0 );
	/* calculate U parameter and test bounds */
	*uu = DOT( tvec, pvec ) * inv_det;
	if( *uu < 0.0 || *uu > 1.0 )
		return 0;
	/* prepare to test V parameter */
	CROSS( qvec, tvec, edge1 );
	/* calculate V parameter and test bounds */
	*vv = DOT( dir, qvec ) * inv_det;
	if( *vv < 0.0 || *uu + *vv > 1.0 )
		return 0;
	/* calculate t, ray intersects triangle */
	*t = DOT( edge2, qvec ) * inv_det;
#endif
	return 1;
}


static void
test_view( void )
{
	int             x;
	VERT            rbase, rdir;
	VERT            sspace;
	char            tmpp[255];
	VERT            s;
	float           WX, WY;
	float           asp;
	float           sz;
	double          ray_base[3];
	double          ray_dir[3];
	VERT            lookat;
	int             tmp_curpoly;
	float           tmp_dist = 1000000.0f;

	tmp_curpoly = Gcurrent_poly;
//  WX=580-800+SXRES;
	WX = work2.x;
	WY = 570 - 580 + SYRES;
	asp = ( WX / WY );
	s.x = G_MOUSEX - 55;
	s.y = G_MOUSEY - 45;
//    s.x/=(580-55);
	s.x /= ( WX - 55 );
	s.y /= ( WY - 45 );
//    s.y/=(570-45);
	s.x *= 2.0 * asp;
	s.y *= 2.0f;
	s.x -= 1.0 * asp;
	s.y -= 1.0f;				// mouse position in view/screnespace
	s.z = -100.0f;
//      v1.x=(v1.x+cam_panx)*cam_scale;
//      v1.y=(v1.y+cam_pany)*cam_scale;
//      v1.z=(v1.z+cam_panz)*cam_scale;
//      v1=vxm(InvCamMatrix,v1);

//  s=vxm(InvCamMatrix,s);
	lookat.x = 0;
	lookat.y = 0;
	lookat.z = 1;
//      tvert.x=(tvert.x+cam_panx)*cam_scale;
	//      tvert.y=(tvert.y+cam_pany)*cam_scale;
	//    tvert.z=(tvert.z+cam_panz)*cam_scale;
//  s=vxm(InvCamMatrix,s);
	lookat = Vnorm( lookat );
	ray_base[0] = s.x;
	ray_base[1] = s.y;
	ray_base[2] = s.z;
	ray_dir[0] = lookat.x;
	ray_dir[1] = lookat.y;
	ray_dir[2] = lookat.z;

	if( Dtotalfaces > 0 )
		for( x = 0; x < Dtotalfaces; x++ )
			if( x != Gcurrent_poly )
			{
				int             g0, g1, g2;
				double          vert0[3], vert1[3], vert2[3];
				VERT            v1, v2, v3;
				double          t, uu, vv;
				int             fs;

				fs = Dface_start[x];
				g0 = Dfacelist[fs];
				g1 = Dfacelist[fs + 1];
				g2 = Dfacelist[fs + 2];
				v1 = Dv[g0];
				v2 = Dv[g1];
				v3 = Dv[g2];
				v1.x = ( v1.x + cam_panx ) * cam_scale;
				v1.y = ( v1.y + cam_pany ) * cam_scale;
				v1.z = ( v1.z + cam_panz ) * cam_scale;
				v1 = vxm( CamMatrix, v1 );
				v2.x = ( v2.x + cam_panx ) * cam_scale;
				v2.y = ( v2.y + cam_pany ) * cam_scale;
				v2.z = ( v2.z + cam_panz ) * cam_scale;
				v2 = vxm( CamMatrix, v2 );

				v3.x = ( v3.x + cam_panx ) * cam_scale;
				v3.y = ( v3.y + cam_pany ) * cam_scale;
				v3.z = ( v3.z + cam_panz ) * cam_scale;
				v3 = vxm( CamMatrix, v3 );

				vert0[0] = v1.x;
				vert0[1] = v1.y;
				vert0[2] = v1.z;
				vert1[0] = v2.x;
				vert1[1] = v2.y;
				vert1[2] = v2.z;
				vert2[0] = v3.x;
				vert2[1] = v3.y;
				vert2[2] = v3.z;
				if( intersect_triangle( ray_base, ray_dir, vert0, vert1, vert2, &t, &uu, &vv ) == 1 )
				{
					if( t < tmp_dist )
					{
						tmp_dist = t;
						Gcurrent_poly = x;
					}
				}
			}
	if( Gcurrent_poly != tmp_curpoly )
	{
		sprintf( tmpp, "base = %f %f %f  dir = %f %f %f  hit = %d", ray_base[0], ray_base[1], ray_base[2], ray_dir[0], ray_dir[1], ray_dir[2], Gcurrent_poly );
		DRAW_STATUS( tmpp );
		Glast_poly = tmp_curpoly;
	}
}




static void
Ssurface_collide( int hh, float radius )
{
	VERT            pp[3];
	VERT            nn[3];
	VERT            avgnorm;
	VERT            vec;
	int             ret = 0;
	int             collide = 0;
	int             q;
	int             x, y, xxx, yyy, zzz, yy;
	int             slg = 0;

				if( hh < total_splines )
					if( COLLISION_METHOD == 1 )
#ifndef LW3D
//      if (skull_sphere>=0)
#endif
						if( DOCOLLIDE )
						{
							vec.x = 0;
							vec.y = 0;
							vec.z = 0;
							if( Dtotalverts > 0 )
								if( group_stat[Shair[hh].mtl] == 0 )
									if( Shair[hh].hide == 0 )
										if( Shair[hh].select[0] )
											if( Shair[hh].restlength > 0 )
//for (q=0;q<1;q++)
											{
												slg = 0;
												if( head >= 0 )
													if( Sresthair[hh].mtl == head )
													{
														slg = 0;
													}
												if( beard >= 0 )
													if( Sresthair[hh].mtl == beard )
													{
														slg = 1;
													}
												if( eyebrow >= 0 )
													if( Sresthair[hh].mtl == eyebrow )
													{
														slg = 2;
													}
												if( eyelash >= 0 )
													if( Sresthair[hh].mtl == eyelash )
													{
														slg = 3;
													}
												if( splines >= 0 )
													if( Sresthair[hh].mtl == splines )
													{
														slg = 4;
													}
												if( Gsurface_collide[slg] )
													for( x = 1; x < 15; x++ )
														if( Shair[hh].pfID[x] == -1 )
														{
//  x=vrt;
															{

																VERT * a = NULL;
																VERT aa;
																float           smallest = 99999.0f;
																float           smallestd = 99999.0f;

																vec.x = 0;
																vec.y = 0;
																vec.z = 0;
																a = &Shair[hh].hv[x];
																aa = Shair[hh].hv[x];
																{
																	float           tx, ty, tz;

																	tx = Shair[hh].hv[x].x;
																	ty = Shair[hh].hv[x].y;
																	tz = Shair[hh].hv[x].z;
																	tx -= cx1;
																	tx /= ( cx2 - cx1 );
																	tx *= PVX;
																	ty -= cy1;
																	ty /= ( cy2 - cy1 );
																	ty *= PVX;
																	tz -= cz1;
																	tz /= ( cz2 - cz1 );
																	tz *= PVX;
																	xxx = ( int ) floor( tx );
																	yyy = ( int ) floor( ty );
																	zzz = ( int ) floor( tz );
																}

																if( ( xxx >= 0 ) && ( yyy >= 0 ) && ( zzz >= 0 ) && ( xxx < PVX ) && ( yyy < PVX ) && ( zzz < PVX ) )
																{
																	if( pvox[xxx][yyy][zzz].totalp > 0 )
																		for( yy = 0; yy < pvox[xxx][yyy][zzz].totalp; yy++ )
																		{
																			int             q;
																			float           d;
																			float           fd;

																			y = pvox[xxx][yyy][zzz].plist[yy];
//          y=0;
																			// plane test;
																			d = VDot( Dface_norm[y], aa ) + Dface_equation[y];
																			fd = d;
																			fd = ( float ) fabs( d );
//          if (0==1)
//          if ((fd> 0.0001)&&(fd< Shair[hh].restlength*(float)x)&&(fd>Shair[hh].restlength/20.0)) // it's behind the plane but closer than the rad
																			if( ( fabs( fd ) > 0.000001 ) && ( fd > Shair[hh].restlength / ( float ) 5.0 ) );	//&&(fd>Shair[hh].restlength/120.0)) // it's behind the plane but closer than the rad
																			if( Dface_end[y] - Dface_start[y] > 2 )
																			{
																				int             g;
																				float           dr;
																				int             i0, i1, i2;
																				double          t, uu, vv;
																				double          vert0[3], vert1[3], vert2[3];
																				double          o[3], vdir[3];

																				g = Dfacelist[Dface_start[y]];
																				vert2[0] = ( double ) Dv[g].x;
																				vert2[1] = ( double ) Dv[g].y;
																				vert2[2] = ( double ) Dv[g].z;
																				g = Dfacelist[Dface_start[y] + 1];
																				vert1[0] = ( double ) Dv[g].x;
																				vert1[1] = ( double ) Dv[g].y;
																				vert1[2] = ( double ) Dv[g].z;
																				g = Dfacelist[Dface_start[y] + 2];
																				vert0[0] = ( double ) Dv[g].x;
																				vert0[1] = ( double ) Dv[g].y;
																				vert0[2] = ( double ) Dv[g].z;
																				o[0] = ( double ) aa.x;
																				o[1] = ( double ) aa.y;
																				o[2] = ( double ) aa.z;
																				vdir[0] = ( double ) Dface_norm[y].x;
																				vdir[1] = ( double ) Dface_norm[y].y;
																				vdir[2] = ( double ) Dface_norm[y].z;
																				if( intersect_triangle( o, vdir, vert0, vert1, vert2, &t, &uu, &vv ) == 1 )
																				{
//                  t= -t;

																					if( ( float ) fabs( d ) < smallest )
																					{
																						//          t-=.00005;
//                      t*=.999999;
																						if( test_stat == 0 )
																							t += ( double ) ( ( Shair[hh].restlength * ( float ) x ) / 16.0 );
																						if( test_stat == 1 )
																							t += ( double ) ( ( Shair[hh].restlength * ( float ) x ) / 16.0 );
																						smallest = ( float ) fabs( d );
																						smallestd = ( float ) ( d );
																						vec.x = ( float ) ( ( t ) * vdir[0] );
																						vec.y = ( float ) ( ( t ) * vdir[1] );
																						vec.z = ( float ) ( ( t ) * vdir[2] );
																						collide = 1;
																					}
//                  ret=1;                  
																				}
																			}
																		}
																}

//if (smallest> Shair[hh].restlength*x*.99999) // Shair is straight
//{
//
//}
//recalc_Shair(hh);
//if (0==1)
//      if (a!=NULL)
																if( Gsurface_collide[slg] == 1 )
																	if( smallestd < ( float ) ( ( Shair[hh].restlength * ( float ) x ) / 16.0 ) )
																		if( Sresthair[hh].pfID[x] == -1 )
																		{
																			Shair[hh].hv[x].x += vec.x * .3;
																			Shair[hh].hv[x].y += vec.y * .3;
																			Shair[hh].hv[x].z += vec.z * .3;
																			if( test_stat != 1 )
																			{
																				Shair[hh].hv[x].x += vec.x * .7;
																				Shair[hh].hv[x].y += vec.y * .7;
																				Shair[hh].hv[x].z += vec.z * .7;
																			}
//hairvelocity[hh].hv[x].x*=.1;
//hairvelocity[hh].hv[x].y*=.1;
//hairvelocity[hh].hv[x].z*=.1;
																			Shairvelocity[hh].hv[x].x += vec.x * .7;
																			Shairvelocity[hh].hv[x].y += vec.y * .7;
																			Shairvelocity[hh].hv[x].z += vec.z * .7;
																		}
//      if (Sresthair[hh].pfID[x]== -1)
//          spare[hh].hv[x]=hair[hh].hv[x];
//if (0==1)     
																//  for(x=0;x<15;x++)
//if (collide==1)
//  spare[hh].hv[x]=hair[hh].hv[x];
															}
														}	// surf collide
//  for (x=0;x<15;x++)
//          if (Sresthair[hh].pfID[x]== -1)
//          spare[hh].hv[x]=Shair[hh].hv[x];
											}	//q
						}
}




#ifdef crap
static void
surface_collide3( int hh, float radius )
{
	VERT            pp[3];
	VERT            nn[3];
	VERT            avgnorm;
	VERT            vec;
	int             doit = 1;
	int             ret = 0;
	int             collide = 0;
	int             q;
	int             x, y, xxx, yyy, zzz, yy;
	int             slg = 0;

//printf("testing surface_collide here\n");

				doit = 1;
#ifndef NOLIB
	if( Gdynamics_mode == 2 )
		doit = 0;
	if( Gdynamics_mode == 0 )
		doit = 0;
	if( Gsculpt_cursors )
		doit = 1;
#endif

//              printf ("doit %d  method %d docllide %d\n",doit,COLLISION_METHOD,DOCOLLIDE);

	if( doit )
		if( COLLISION_METHOD == 1 )
#ifndef LW3D
//      if (skull_sphere>=0)
#endif
			if( DOCOLLIDE )
			{
				vec.x = 0;
				vec.y = 0;
				vec.z = 0;
//  printf ("looking for collision\n");
				if( Dtotalverts > 0 )
								if( hair[hh].restlength > 0 )
//for (q=0;q<1;q++)
								{
									int             gotone = 0;

									slg = 0;
									if( head >= 0 )
										if( hair[hh].mtl == head )
										{
											slg = 0;
										}
									if( beard >= 0 )
										if( hair[hh].mtl == beard )
										{
											slg = 1;
										}
									if( eyebrow >= 0 )
										if( hair[hh].mtl == eyebrow )
										{
											slg = 2;
										}
									if( eyelash >= 0 )
										if( hair[hh].mtl == eyelash )
										{
											slg = 3;
										}
									if( splines >= 0 )
										if( hair[hh].mtl == splines )
										{
											slg = 4;
										}
									if( Gsurface_collide[slg] )
										for( x = 1; x < 15; x++ )
											if( hair[hh].pfID[x] == -1 )
											{
//  x=vrt;
												{

													VERT           *a = NULL;
													VERT            aa;
													float           smallest = 99999.0f;
													float           smallestd = 99999.0f;
													int tp;
													int *list=NULL; // pointer to list

													vec.x = 0;
													vec.y = 0;
													vec.z = 0;
													a = &hair[hh].hv[x];
													aa = hair[hh].hv[x];

													tp=vlad_request_polys_to_test(&list,aa,hair[hh].restlength / ( float ) 5.0);

													{
															for( yy = 0; yy < tp; yy++ )
															{
																int             q;
																float           d;
																float           fd;

																y = list[yy];
//          y=0;
																// plane test;
																d = VDot( Dface_norm[y], aa ) + Dface_equation[y];
//          if (0==1)
//          if ((fd> 0.0001)&&(fd< hair[hh].restlength*(float)x)&&(fd>hair[hh].restlength/20.0)) // it's behind the plane but closer than the rad
																if( ( fabs( fd ) > 0.000001 )  );	//&&(fd>hair[hh].restlength/120.0)) // it's behind the plane but closer than the rad
																if( Dface_end[y] - Dface_start[y] > 2 )
																{
																	int             g;
																	float           dr;
																	int             i0, i1, i2;
																	double          t, uu, vv;
																	double          vert0[3], vert1[3], vert2[3];
																	double          o[3], vdir[3];

																	g = Dfacelist[Dface_start[y]];
																	vert2[0] = ( double ) Dv[g].x;
																	vert2[1] = ( double ) Dv[g].y;
																	vert2[2] = ( double ) Dv[g].z;
																	g = Dfacelist[Dface_start[y] + 1];
																	vert1[0] = ( double ) Dv[g].x;
																	vert1[1] = ( double ) Dv[g].y;
																	vert1[2] = ( double ) Dv[g].z;
																	g = Dfacelist[Dface_start[y] + 2];
																	vert0[0] = ( double ) Dv[g].x;
																	vert0[1] = ( double ) Dv[g].y;
																	vert0[2] = ( double ) Dv[g].z;
																	o[0] = ( double ) aa.x;
																	o[1] = ( double ) aa.y;
																	o[2] = ( double ) aa.z;
																	vdir[0] = ( double ) Dface_norm[y].x;
																	vdir[1] = ( double ) Dface_norm[y].y;
																	vdir[2] = ( double ) Dface_norm[y].z;
																	if( intersect_triangle( o, vdir, vert0, vert1, vert2, &t, &uu, &vv ) == 1 )
																	{
//                  t= -t;

																		if( ( float ) fabs( d ) < smallest )
																		{
																			//          t-=.00005;
//                      t*=.999999;
																			if( test_stat == 0 )
																				t += ( double ) ( ( hair[hh].restlength * ( float ) 1 ) / 6.0 );
																			if( test_stat == 1 )
																				t += ( double ) ( ( hair[hh].restlength * ( float ) 1 ) / 6.0 );
																			smallest = ( float ) fabs( d );
																			smallestd = ( float ) ( d );
																			vec.x = ( float ) ( ( t ) * vdir[0] );
																			vec.y = ( float ) ( ( t ) * vdir[1] );
																			vec.z = ( float ) ( ( t ) * vdir[2] );
																			collide = 1;
//                  hair[hh].velocity[x].x=0;
//                  hair[hh].velocity[x].y=0;
//                  hair[hh].velocity[x].z=0;
//                  printf("colliding here\n");

																		}
//                  ret=1;                  
																	}
																}
															}
													}

//if (smallest> hair[hh].restlength*x*.99999) // hair is straight
//{
//
//}
//recalc_hair(hh);
//if (0==1)
//      if (a!=NULL)


													{
														float           dd;
														VERT            df;

														df.x = hair[hh].hv[x - 1].x - hair[hh].hv[x].x;
														df.y = hair[hh].hv[x - 1].y - hair[hh].hv[x].y;
														df.z = hair[hh].hv[x - 1].z - hair[hh].hv[x].z;
														df = Vnorm( df );
														dd = sqrt( vec.x * vec.x + vec.y * vec.y + vec.z * vec.z );
														if( dd < hair[hh].restlength * 15.0f )	// don't correct for distances over a couple links in length
															if( VDot( df, vec ) )	// is the correction vector in the same hemisphere as the vector to the last vert?

																if( Gsurface_collide[slg] == 1 )
																	if( smallestd < ( float ) ( ( hair[hh].restlength * ( float ) 1 ) / 12.0 ) )
																		if( hair[hh].pfID[x] == -1 )
																		{
																			hair[hh].hv[x].x += vec.x * .83;
																			hair[hh].hv[x].y += vec.y * .83;
																			hair[hh].hv[x].z += vec.z * .83;
																			if( test_stat != 1 )
																			{
																				hair[hh].hv[x].x += vec.x * .3;
																				hair[hh].hv[x].y += vec.y * .3;
																				hair[hh].hv[x].z += vec.z * .3;
																			}
																			// correct the length
																			df.x = hair[hh].hv[x].x - hair[hh].hv[x - 1].x;
																			df.y = hair[hh].hv[x].y - hair[hh].hv[x - 1].y;
																			df.z = hair[hh].hv[x].z - hair[hh].hv[x - 1].z;
																			df = Vnorm( df );
																			hair[hh].hv[x].x = hair[hh].hv[x - 1].x + df.x * hair[hh].restlength;
																			hair[hh].hv[x].y = hair[hh].hv[x - 1].y + df.y * hair[hh].restlength;
																			hair[hh].hv[x].z = hair[hh].hv[x - 1].z + df.z * hair[hh].restlength;

																			Gcollision_hit = 1;
																			gotone = 1;
//hairvelocity[hh].hv[x].x*=.1;
//hairvelocity[hh].hv[x].y*=.1;
//hairvelocity[hh].hv[x].z*=.1;
//      if (resthair[hh].pfID[x]== -1)
//          spare[hh].hv[x]=hair[hh].hv[x];
																			hair[hh].velocity[x].x += vec.x * .97;
																			hair[hh].velocity[x].y += vec.y * .97;
																			hair[hh].velocity[x].z += vec.z * .97;
																		}
													}

//if (0==1)     
													//  for(x=0;x<15;x++)
//if (collide==1)
//  spare[hh].hv[x]=hair[hh].hv[x];
													free(list);
												}
											}	// surf collide

								}	//q
			}
}
#endif