+ v3f ab, ac, ap;
+ f32 d1, d2;
+
+ /* Region outside A */
+ v3_sub( tri[1], tri[0], ab );
+ v3_sub( tri[2], tri[0], ac );
+ v3_sub( p, tri[0], ap );
+
+ d1 = v3_dot(ab,ap);
+ d2 = v3_dot(ac,ap);
+ if( d1 <= 0.0f && d2 <= 0.0f )
+ {
+ v3_copy( tri[0], dest );
+ return k_contact_type_default;
+ }
+
+ /* Region outside B */
+ v3f bp;
+ f32 d3, d4;
+
+ v3_sub( p, tri[1], bp );
+ d3 = v3_dot( ab, bp );
+ d4 = v3_dot( ac, bp );
+
+ if( d3 >= 0.0f && d4 <= d3 )
+ {
+ v3_copy( tri[1], dest );
+ return k_contact_type_edge;
+ }
+
+ /* Edge region of AB */
+ f32 vc = d1*d4 - d3*d2;
+ if( vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f )
+ {
+ f32 v = d1 / (d1-d3);
+ v3_muladds( tri[0], ab, v, dest );
+ return k_contact_type_edge;
+ }
+
+ /* Region outside C */
+ v3f cp;
+ f32 d5, d6;
+ v3_sub( p, tri[2], cp );
+ d5 = v3_dot(ab, cp);
+ d6 = v3_dot(ac, cp);
+
+ if( d6 >= 0.0f && d5 <= d6 )
+ {
+ v3_copy( tri[2], dest );
+ return k_contact_type_edge;
+ }
+
+ /* Region of AC */
+ f32 vb = d5*d2 - d1*d6;
+ if( vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f )
+ {
+ f32 w = d2 / (d2-d6);
+ v3_muladds( tri[0], ac, w, dest );
+ return k_contact_type_edge;
+ }
+
+ /* Region of BC */
+ f32 va = d3*d6 - d5*d4;
+ if( va <= 0.0f && (d4-d3) >= 0.0f && (d5-d6) >= 0.0f )
+ {
+ f32 w = (d4-d3) / ((d4-d3) + (d5-d6));
+ v3f bc;
+ v3_sub( tri[2], tri[1], bc );
+ v3_muladds( tri[1], bc, w, dest );
+ return k_contact_type_edge;
+ }
+
+ /* P inside region, Q via barycentric coordinates uvw */
+ f32 d = 1.0f/(va+vb+vc),
+ v = vb*d,
+ w = vc*d;
+
+ v3_muladds( tri[0], ab, v, dest );
+ v3_muladds( dest, ac, w, dest );
+
+ return k_contact_type_default;
+}
+
+static void closest_point_elipse( v2f p, v2f e, v2f o )
+{
+ v2f pabs, ei, e2, ve, t;
+
+ v2_abs( p, pabs );
+ v2_div( (v2f){ 1.0f, 1.0f }, e, ei );
+ v2_mul( e, e, e2 );
+ v2_mul( ei, (v2f){ e2[0]-e2[1], e2[1]-e2[0] }, ve );
+
+ v2_fill( t, 0.70710678118654752f );
+
+ for( int i=0; i<3; i++ ){
+ v2f v, u, ud, w;
+
+ v2_mul( ve, t, v ); /* ve*t*t*t */
+ v2_mul( v, t, v );
+ v2_mul( v, t, v );
+
+ v2_sub( pabs, v, u );
+ v2_normalize( u );
+
+ v2_mul( t, e, ud );
+ v2_sub( ud, v, ud );
+
+ v2_muls( u, v2_length( ud ), u );
+
+ v2_add( v, u, w );
+ v2_mul( w, ei, w );
+
+ v2_maxv( (v2f){0.0f,0.0f}, w, t );
+ v2_normalize( t );
+ }
+
+ v2_mul( t, e, o );
+ v2_copysign( o, p );
+}
+
+/*
+ * -----------------------------------------------------------------------------
+ * Section 5.d Raycasts & Spherecasts
+ * -----------------------------------------------------------------------------
+ */
+
+int ray_aabb1( boxf box, v3f co, v3f dir_inv, f32 dist )
+{
+ v3f v0, v1;
+ f32 tmin, tmax;
+
+ v3_sub( box[0], co, v0 );
+ v3_sub( box[1], co, v1 );
+
+ v3_mul( v0, dir_inv, v0 );
+ v3_mul( v1, dir_inv, v1 );
+
+ tmin = vg_minf( v0[0], v1[0] );
+ tmax = vg_maxf( v0[0], v1[0] );
+ tmin = vg_maxf( tmin, vg_minf( v0[1], v1[1] ));
+ tmax = vg_minf( tmax, vg_maxf( v0[1], v1[1] ));
+ tmin = vg_maxf( tmin, vg_minf( v0[2], v1[2] ));
+ tmax = vg_minf( tmax, vg_maxf( v0[2], v1[2] ));
+
+ return (tmax >= tmin) && (tmin <= dist) && (tmax >= 0.0f);
+}
+
+/* Time of intersection with ray vs triangle */
+static int ray_tri( v3f tri[3], v3f co,
+ v3f dir, f32 *dist )
+{
+ f32 const kEpsilon = 0.00001f;