X-Git-Url: http://git.mutantstargoat.com/user/nuclear/?p=csgray;a=blobdiff_plain;f=src%2Fgeom.c;h=1e6110cedb4216da22feb3d4649a1b7e81fb466d;hp=1961aed1ee871ef53a1c5fedbc639eb8095ee853;hb=HEAD;hpb=7abe0c51fc87e90674427c5ac5848f394af4575a

diff --git a/src/geom.c b/src/geom.c
index 1961aed..1e6110c 100644
--- a/src/geom.c
+++ b/src/geom.c
@@ -1,15 +1,21 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <float.h>
+#include <assert.h>
 #include "geom.h"
 #include "matrix.h"
 
 #define EPSILON		1e-6f
 
+static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b);
+static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b);
+static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b);
+
 /* TODO custom hit allocator */
-struct hit *alloc_hit(void)
+struct hinterv *alloc_hit(void)
 {
-	struct hit *hit = calloc(sizeof *hit, 1);
+	struct hinterv *hit = calloc(sizeof *hit, 1);
 	if(!hit) {
 		perror("failed to allocate ray hit node");
 		abort();
@@ -17,13 +23,13 @@ struct hit *alloc_hit(void)
 	return hit;
 }
 
-struct hit *alloc_hits(int n)
+struct hinterv *alloc_hits(int n)
 {
 	int i;
-	struct hit *list = 0;
+	struct hinterv *list = 0;
 
 	for(i=0; i<n; i++) {
-		struct hit *hit = alloc_hit();
+		struct hinterv *hit = alloc_hit();
 		hit->next = list;
 		list = hit;
 	}
@@ -31,21 +37,21 @@ struct hit *alloc_hits(int n)
 }
 
 
-void free_hit(struct hit *hit)
+void free_hit(struct hinterv *hit)
 {
 	free(hit);
 }
 
-void free_hit_list(struct hit *hit)
+void free_hit_list(struct hinterv *hit)
 {
 	while(hit) {
-		struct hit *tmp = hit;
+		struct hinterv *tmp = hit;
 		hit = hit->next;
 		free_hit(tmp);
 	}
 }
 
-struct hit *ray_intersect(struct ray *ray, csg_object *o)
+struct hinterv *ray_intersect(struct ray *ray, csg_object *o)
 {
 	switch(o->ob.type) {
 	case OB_SPHERE:
@@ -54,6 +60,8 @@ struct hit *ray_intersect(struct ray *ray, csg_object *o)
 		return ray_cylinder(ray, o);
 	case OB_PLANE:
 		return ray_plane(ray, o);
+	case OB_BOX:
+		return ray_box(ray, o);
 	case OB_UNION:
 		return ray_csg_un(ray, o);
 	case OB_INTERSECTION:
@@ -66,11 +74,11 @@ struct hit *ray_intersect(struct ray *ray, csg_object *o)
 	return 0;
 }
 
-struct hit *ray_sphere(struct ray *ray, csg_object *o)
+struct hinterv *ray_sphere(struct ray *ray, csg_object *o)
 {
 	int i;
 	float a, b, c, d, sqrt_d, t[2], sq_rad, tmp;
-	struct hit *hit, *hitlist;
+	struct hinterv *hit;
 	struct ray locray = *ray;
 
 	if(o->sph.rad == 0.0f) {
@@ -101,37 +109,172 @@ struct hit *ray_sphere(struct ray *ray, csg_object *o)
 		t[1] = tmp;
 	}
 
-	hitlist = hit = alloc_hits(2);
+	hit = alloc_hits(1);
+	hit->o = o;
 	for(i=0; i<2; i++) {
 		float c[3] = {0, 0, 0};
+		float x, y, z;
+
 		mat4_xform3(c, o->ob.xform, c);
 
-		hit->t = t[i];
-		hit->x = ray->x + ray->dx * t[i];
-		hit->y = ray->y + ray->dy * t[i];
-		hit->z = ray->z + ray->dz * t[i];
-		hit->nx = (hit->x - c[0]) / o->sph.rad;
-		hit->ny = (hit->y - c[1]) / o->sph.rad;
-		hit->nz = (hit->z - c[2]) / o->sph.rad;
-		hit->o = o;
+		x = ray->x + ray->dx * t[i];
+		y = ray->y + ray->dy * t[i];
+		z = ray->z + ray->dz * t[i];
+
+		hit->end[i].t = t[i];
+		hit->end[i].x = x;
+		hit->end[i].y = y;
+		hit->end[i].z = z;
+		hit->end[i].nx = (x - c[0]) / o->sph.rad;
+		hit->end[i].ny = (y - c[1]) / o->sph.rad;
+		hit->end[i].nz = (z - c[2]) / o->sph.rad;
+		hit->end[i].o = o;
+	}
+	return hit;
+}
 
-		hit = hit->next;
+static int ray_cylcap(struct ray *ray, float y, float rad, float *tres)
+{
+	float ndotr, ndotv, vy, t;
+	float ny = y > 0.0f ? 1.0f : -1.0f;
+	float x, z, lensq;
+
+	ndotr = ny * ray->dy;
+	if(fabs(ndotr) < EPSILON) return 0;
+
+	vy = y - ray->y;
+
+	ndotv = ny * vy;
+
+	t = ndotv / ndotr;
+
+	x = ray->x + ray->dx * t;
+	z = ray->z + ray->dz * t;
+	lensq = x * x + z * z;
+
+	if(lensq <= rad * rad) {
+		*tres = t;
+		return 1;
 	}
-	return hitlist;
+	return 0;
 }
 
-struct hit *ray_cylinder(struct ray *ray, csg_object *o)
+struct hinterv *ray_cylinder(struct ray *ray, csg_object *o)
 {
+	int i, out[2] = {0}, t_is_cap[2] = {0};
+	float a, b, c, d, sqrt_d, t[2], sq_rad, tmp, y[2], hh, cap_t;
+	struct hinterv *hit;
 	struct ray locray = *ray;
 
+	if(o->cyl.rad == 0.0f || o->cyl.height == 0.0f) {
+		return 0;
+	}
+	sq_rad = o->cyl.rad * o->cyl.rad;
+	hh = o->cyl.height / 2.0f;
+
 	xform_ray(&locray, o->ob.inv_xform);
-	return 0;	/* TODO */
+
+	a = locray.dx * locray.dx + locray.dz * locray.dz;
+	b = 2.0f * (locray.dx * locray.x + locray.dz * locray.z);
+	c = locray.x * locray.x + locray.z * locray.z - sq_rad;
+
+	d = b * b - 4.0f * a * c;
+	if(d < EPSILON) return 0;
+
+	sqrt_d = sqrt(d);
+	t[0] = (-b + sqrt_d) / (2.0f * a);
+	t[1] = (-b - sqrt_d) / (2.0f * a);
+
+	if(t[0] < EPSILON && t[1] < EPSILON) {
+		return 0;
+	}
+	if(t[1] < t[0]) {
+		tmp = t[0];
+		t[0] = t[1];
+		t[1] = tmp;
+	}
+
+	y[0] = locray.y + locray.dy * t[0];
+	y[1] = locray.y + locray.dy * t[1];
+
+	if(y[0] < -hh || y[0] > hh) {
+		out[0] = 1;
+	}
+	if(y[1] < -hh || y[1] > hh) {
+		out[1] = 1;
+	}
+
+	if(out[0]) {
+		t[0] = t[1];
+	}
+	if(out[1]) {
+		t[1] = t[0];
+	}
+
+	if(ray_cylcap(ray, hh, o->cyl.rad, &cap_t)) {
+		if(cap_t < t[0]) {
+			t[0] = cap_t;
+			t_is_cap[0] = 1;
+			out[0] = 0;
+		}
+		if(cap_t > t[1]) {
+			t[1] = cap_t;
+			t_is_cap[1] = 1;
+			out[1] = 0;
+		}
+	}
+	if(ray_cylcap(ray, -hh, o->cyl.rad, &cap_t)) {
+		if(cap_t < t[0]) {
+			t[0] = cap_t;
+			t_is_cap[0] = -1;
+			out[0] = 0;
+		}
+		if(cap_t > t[1]) {
+			t[1] = cap_t;
+			t_is_cap[1] = -1;
+			out[1] = 0;
+		}
+	}
+
+	if(out[0] && out[1]) {
+		return 0;
+	}
+
+	hit = alloc_hits(1);
+	hit->o = o;
+	for(i=0; i<2; i++) {
+		float c[3] = {0, 0, 0};
+		float x, y, z;
+
+		x = ray->x + ray->dx * t[i];
+		y = ray->y + ray->dy * t[i];
+		z = ray->z + ray->dz * t[i];
+
+		if(t_is_cap[i]) {
+			hit->end[i].nx = hit->end[i].nz = 0.0f;
+			hit->end[i].ny = t_is_cap[i] > 0 ? 1.0f : -1.0f;
+		} else {
+			c[1] = locray.y + locray.dy * t[i];
+			mat4_xform3(c, o->ob.xform, c);
+
+			hit->end[i].nx = (x - c[0]) / o->cyl.rad;
+			hit->end[i].ny = (y - c[1]) / o->cyl.rad;
+			hit->end[i].nz = (z - c[2]) / o->cyl.rad;
+		}
+
+		hit->end[i].t = t[i];
+		hit->end[i].x = x;
+		hit->end[i].y = y;
+		hit->end[i].z = z;
+		hit->end[i].o = o;
+	}
+	return hit;
 }
 
-struct hit *ray_plane(struct ray *ray, csg_object *o)
+struct hinterv *ray_plane(struct ray *ray, csg_object *o)
 {
 	float vx, vy, vz, ndotv, ndotr, t;
-	struct hit *hit = 0;
+	struct hinterv *hit;
 	struct ray locray = *ray;
 
 	xform_ray(&locray, o->ob.inv_xform);
@@ -146,24 +289,113 @@ struct hit *ray_plane(struct ray *ray, csg_object *o)
 	ndotv = o->plane.nx * vx + o->plane.ny * vy + o->plane.nz * vz;
 
 	t = ndotv / ndotr;
+	if(t < EPSILON) {
+		return 0;
+	}
+
+	hit = alloc_hits(1);
+	hit->o = hit->end[0].o = hit->end[1].o = o;
+	hit->end[0].t = t;
+	hit->end[0].x = ray->x + ray->dx * t;
+	hit->end[0].y = ray->y + ray->dy * t;
+	hit->end[0].z = ray->z + ray->dz * t;
+
+	hit->end[0].nx = hit->end[1].nx = o->plane.nx;
+	hit->end[0].ny = hit->end[1].ny = o->plane.ny;
+	hit->end[0].nz = hit->end[1].nz = o->plane.nz;
+
+	hit->end[1].t = FLT_MAX;
+	hit->end[1].x = ray->x + ray->dx * 10000.0f;
+	hit->end[1].y = ray->y + ray->dy * 10000.0f;
+	hit->end[1].z = ray->z + ray->dz * 10000.0f;
+	return hit;
+}
+
+#define BEXT(x)	((x) * 0.49999)
+
+struct hinterv *ray_box(struct ray *ray, csg_object *o)
+{
+	int i, sign[3];
+	float param[2][3];
+	float inv_dir[3];
+	float tmin, tmax, tymin, tymax, tzmin, tzmax;
+	struct hinterv *hit;
+	struct ray locray = *ray;
+	float dirmat[16];
+
+	xform_ray(&locray, o->ob.inv_xform);
+
+	for(i=0; i<3; i++) {
+		float sz = *(&o->box.xsz + i);
+		param[0][i] = -0.5 * sz;
+		param[1][i] = 0.5 * sz;
+
+		inv_dir[i] = 1.0f / *(&locray.dx + i);
+		sign[i] = inv_dir[i] < 0;
+	}
+
+	tmin = (param[sign[0]][0] - locray.x) * inv_dir[0];
+	tmax = (param[1 - sign[0]][0] - locray.x) * inv_dir[0];
+	tymin = (param[sign[1]][1] - locray.y) * inv_dir[1];
+	tymax = (param[1 - sign[1]][1] - locray.y) * inv_dir[1];
 
-	if(t > EPSILON) {
-		hit = alloc_hits(1);
-		hit->t = t;
-		hit->x = ray->x + ray->dx * t;
-		hit->y = ray->y + ray->dy * t;
-		hit->z = ray->z + ray->dz * t;
-		hit->nx = o->plane.nx;
-		hit->ny = o->plane.ny;
-		hit->nz = o->plane.nz;
-		hit->o = o;
+	if(tmin > tymax || tymin > tmax) {
+		return 0;
+	}
+	if(tymin > tmin) {
+		tmin = tymin;
+	}
+	if(tymax < tmax) {
+		tmax = tymax;
+	}
+
+	tzmin = (param[sign[2]][2] - locray.z) * inv_dir[2];
+	tzmax = (param[1 - sign[2]][2] - locray.z) * inv_dir[2];
+
+	if(tmin > tzmax || tzmin > tmax) {
+		return 0;
+	}
+	if(tzmin > tmin) {
+		tmin = tzmin;
+	}
+	if(tzmax < tmax) {
+		tmax = tzmax;
+	}
+
+	mat4_copy(dirmat, o->ob.xform);
+	mat4_upper3x3(dirmat);
+
+	hit = alloc_hits(1);
+	hit->o = o;
+	for(i=0; i<2; i++) {
+		float n[3] = {0};
+		float t = i == 0 ? tmin : tmax;
+
+		float x = (locray.x + locray.dx * t) / o->box.xsz;
+		float y = (locray.y + locray.dy * t) / o->box.ysz;
+		float z = (locray.z + locray.dz * t) / o->box.zsz;
+
+		if(fabs(x) > fabs(y) && fabs(x) > fabs(z)) {
+			n[0] = x > 0.0f ? 1.0f : -1.0f;
+		} else if(fabs(y) > fabs(z)) {
+			n[1] = y > 0.0f ? 1.0f : -1.0f;
+		} else {
+			n[2] = z > 0.0f ? 1.0f : -1.0f;
+		}
+
+		hit->end[i].o = o;
+		hit->end[i].t = t;
+		hit->end[i].x = ray->x + ray->dx * t;
+		hit->end[i].y = ray->y + ray->dy * t;
+		hit->end[i].z = ray->z + ray->dz * t;
+		mat4_xform3(&hit->end[i].nx, dirmat, n);
 	}
 	return hit;
 }
 
-struct hit *ray_csg_un(struct ray *ray, csg_object *o)
+struct hinterv *ray_csg_un(struct ray *ray, csg_object *o)
 {
-	struct hit *hita, *hitb;
+	struct hinterv *hita, *hitb, *res;
 
 	hita = ray_intersect(ray, o->un.a);
 	hitb = ray_intersect(ray, o->un.b);
@@ -171,73 +403,166 @@ struct hit *ray_csg_un(struct ray *ray, csg_object *o)
 	if(!hita) return hitb;
 	if(!hitb) return hita;
 
-	if(hita->t < hitb->t) {
+	res = interval_union(hita, hitb);
+	free_hit_list(hita);
+	free_hit_list(hitb);
+	return res;
+}
+
+struct hinterv *ray_csg_isect(struct ray *ray, csg_object *o)
+{
+	struct hinterv *hita, *hitb, *res;
+
+	hita = ray_intersect(ray, o->isect.a);
+	hitb = ray_intersect(ray, o->isect.b);
+
+	if(!hita || !hitb) {
+		free_hit_list(hita);
 		free_hit_list(hitb);
-		return hita;
+		return 0;
 	}
+
+	res = interval_isect(hita, hitb);
 	free_hit_list(hita);
-	return hitb;
+	free_hit_list(hitb);
+	return res;
 }
 
-struct hit *ray_csg_isect(struct ray *ray, csg_object *o)
+struct hinterv *ray_csg_sub(struct ray *ray, csg_object *o)
 {
-	return 0;
+	struct hinterv *hita, *hitb, *res;
+
+	hita = ray_intersect(ray, o->un.a);
+	hitb = ray_intersect(ray, o->un.b);
+
+	if(!hita) return 0;
+	if(!hitb) return hita;
+
+	res = interval_sub(hita, hitb);
+	free_hit_list(hita);
+	free_hit_list(hitb);
+	return res;
+}
+
+
+void xform_ray(struct ray *ray, float *mat)
+{
+	float m3x3[16];
+
+	mat4_copy(m3x3, mat);
+	mat4_upper3x3(m3x3);
+
+	mat4_xform3(&ray->x, mat, &ray->x);
+	mat4_xform3(&ray->dx, m3x3, &ray->dx);
 }
 
-static struct hit *first(struct hit *hlist)
+static void flip_hit(struct hit *hit)
 {
-	return hlist;
+	hit->nx = -hit->nx;
+	hit->ny = -hit->ny;
+	hit->nz = -hit->nz;
 }
 
-static struct hit *last(struct hit *hlist)
+static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b)
 {
-	while(hlist->next) {
-		hlist = hlist->next;
+	struct hinterv *res, *res2;
+
+	if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
+		/* disjoint */
+		res = alloc_hits(2);
+		res2 = res->next;
+
+		if(a->end[0].t < b->end[0].t) {
+			*res = *a;
+			*res2 = *b;
+		} else {
+			*res = *b;
+			*res2 = *a;
+		}
+		res->next = res2;
+		res2->next = 0;
+		return res;
 	}
-	return hlist;
+
+	res = alloc_hits(1);
+	res->end[0] = a->end[0].t <= b->end[0].t ? a->end[0] : b->end[0];
+	res->end[1] = a->end[1].t >= b->end[1].t ? a->end[1] : b->end[1];
+	return res;
 }
 
-struct hit *ray_csg_sub(struct ray *ray, csg_object *o)
+static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b)
 {
-	struct hit *hita, *hitb, *lasthit, tmp;
+	struct hinterv *res;
 
-	hita = ray_intersect(ray, o->sub.a);
-	hitb = ray_intersect(ray, o->sub.b);
+	if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
+		/* disjoint */
+		return 0;
+	}
 
-	if(!hita) return 0;
-	if(!hitb) return hita;
+	res = alloc_hits(1);
 
-	if(first(hita)->t < first(hitb)->t) {
-		free_hit_list(hitb);
-		return hita;
+	if(a->end[0].t <= b->end[0].t && a->end[1].t >= b->end[1].t) {
+		/* B in A */
+		*res = *a;
+		res->next = 0;
+		return res;
+	}
+	if(a->end[0].t > b->end[0].t && a->end[1].t < b->end[1].t) {
+		/* A in B */
+		*res = *b;
+		res->next = 0;
+		return res;
 	}
-	if(first(hita)->t < (lasthit = last(hitb))->t) {
-		/* overlapping */
-		tmp = *hitb;
-		*hitb = *lasthit;
 
-		free_hit_list(tmp.next);
+	/* partial overlap */
+	if(a->end[0].t < b->end[0].t) {
+		res->end[0] = b->end[0];
+		res->end[1] = a->end[1];
+	} else {
+		res->end[0] = a->end[0];
+		res->end[1] = b->end[1];
+	}
+	return res;
+}
 
-		hitb->nx = -hitb->nx;
-		hitb->ny = -hitb->ny;
-		hitb->nz = -hitb->nz;
+static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b)
+{
+	struct hinterv *res;
 
-		free_hit_list(hita);
-		return hitb;
+	if(a->end[0].t >= b->end[0].t && a->end[1].t <= b->end[1].t) {
+		/* A in B */
+		return 0;
 	}
 
-	free_hit_list(hitb);
-	return hita;
-}
+	if(a->end[0].t < b->end[0].t && a->end[1].t > b->end[1].t) {
+		/* B in A */
+		res = alloc_hits(2);
+		res->end[0] = a->end[0];
+		res->end[1] = b->end[0];
+		res->next->end[0] = b->end[1];
+		res->next->end[1] = a->end[1];
+		return res;
+	}
 
+	res = alloc_hits(1);
 
-void xform_ray(struct ray *ray, float *mat)
-{
-	float m3x3[16];
+	if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
+		/* disjoint */
+		*res = *a;
+		res->next = 0;
+		return res;
+	}
 
-	mat4_copy(m3x3, mat);
-	mat4_upper3x3(m3x3);
+	/* partial overlap */
+	if(a->end[0].t <= b->end[0].t) {
+		res->end[0] = a->end[0];
+		res->end[1] = b->end[0];
+	} else {
+		res->end[0] = b->end[1];
+		res->end[1] = a->end[1];
+	}
 
-	mat4_xform3(&ray->x, mat, &ray->x);
-	mat4_xform3(&ray->dx, m3x3, &ray->dx);
+	flip_hit(res->end + 0);
+	flip_hit(res->end + 1);
+	return res;
 }