git.mutantstargoat.com Git - csgray/blob - src/geom.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <float.h>
   5 #include "geom.h"
   6 #include "matrix.h"
   7
   8 #define EPSILON         1e-6f
   9
  10 static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b);
  11 static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b);
  12 static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b);
  13
  14 /* TODO custom hit allocator */
  15 struct hinterv *alloc_hit(void)
  16 {
  17         struct hinterv *hit = calloc(sizeof *hit, 1);
  18         if(!hit) {
  19                 perror("failed to allocate ray hit node");
  20                 abort();
  21         }
  22         return hit;
  23 }
  24
  25 struct hinterv *alloc_hits(int n)
  26 {
  27         int i;
  28         struct hinterv *list = 0;
  29
  30         for(i=0; i<n; i++) {
  31                 struct hinterv *hit = alloc_hit();
  32                 hit->next = list;
  33                 list = hit;
  34         }
  35         return list;
  36 }
  37
  38
  39 void free_hit(struct hinterv *hit)
  40 {
  41         free(hit);
  42 }
  43
  44 void free_hit_list(struct hinterv *hit)
  45 {
  46         while(hit) {
  47                 struct hinterv *tmp = hit;
  48                 hit = hit->next;
  49                 free_hit(tmp);
  50         }
  51 }
  52
  53 struct hinterv *ray_intersect(struct ray *ray, csg_object *o)
  54 {
  55         switch(o->ob.type) {
  56         case OB_SPHERE:
  57                 return ray_sphere(ray, o);
  58         case OB_CYLINDER:
  59                 return ray_cylinder(ray, o);
  60         case OB_PLANE:
  61                 return ray_plane(ray, o);
  62         case OB_UNION:
  63                 return ray_csg_un(ray, o);
  64         case OB_INTERSECTION:
  65                 return ray_csg_isect(ray, o);
  66         case OB_SUBTRACTION:
  67                 return ray_csg_sub(ray, o);
  68         default:
  69                 break;
  70         }
  71         return 0;
  72 }
  73
  74 struct hinterv *ray_sphere(struct ray *ray, csg_object *o)
  75 {
  76         int i;
  77         float a, b, c, d, sqrt_d, t[2], sq_rad, tmp;
  78         struct hinterv *hit;
  79         struct ray locray = *ray;
  80
  81         if(o->sph.rad == 0.0f) {
  82                 return 0;
  83         }
  84         sq_rad = o->sph.rad * o->sph.rad;
  85
  86         xform_ray(&locray, o->ob.inv_xform);
  87
  88         a = locray.dx * locray.dx + locray.dy * locray.dy + locray.dz * locray.dz;
  89         b = 2.0f * (locray.dx * locray.x + locray.dy * locray.y + locray.dz * locray.z);
  90         c = (locray.x * locray.x + locray.y * locray.y + locray.z * locray.z) - sq_rad;
  91
  92         d = b * b - 4.0f * a * c;
  93         if(d < EPSILON) return 0;
  94
  95         sqrt_d = sqrt(d);
  96         t[0] = (-b + sqrt_d) / (2.0f * a);
  97         t[1] = (-b - sqrt_d) / (2.0f * a);
  98
  99         if(t[0] < EPSILON && t[1] < EPSILON) {
 100                 return 0;
 101         }
 102
 103         if(t[1] < t[0]) {
 104                 tmp = t[0];
 105                 t[0] = t[1];
 106                 t[1] = tmp;
 107         }
 108
 109         hit = alloc_hits(1);
 110         hit->o = o;
 111         for(i=0; i<2; i++) {
 112                 float c[3] = {0, 0, 0};
 113                 float x, y, z;
 114
 115                 mat4_xform3(c, o->ob.xform, c);
 116
 117                 x = ray->x + ray->dx * t[i];
 118                 y = ray->y + ray->dy * t[i];
 119                 z = ray->z + ray->dz * t[i];
 120
 121                 hit->end[i].t = t[i];
 122                 hit->end[i].x = x;
 123                 hit->end[i].y = y;
 124                 hit->end[i].z = z;
 125                 hit->end[i].nx = (x - c[0]) / o->sph.rad;
 126                 hit->end[i].ny = (y - c[1]) / o->sph.rad;
 127                 hit->end[i].nz = (z - c[2]) / o->sph.rad;
 128                 hit->end[i].o = o;
 129         }
 130         return hit;
 131 }
 132
 133 struct hinterv *ray_cylinder(struct ray *ray, csg_object *o)
 134 {
 135         struct ray locray = *ray;
 136
 137         xform_ray(&locray, o->ob.inv_xform);
 138         return 0;       /* TODO */
 139 }
 140
 141 struct hinterv *ray_plane(struct ray *ray, csg_object *o)
 142 {
 143         float vx, vy, vz, ndotv, ndotr, t;
 144         struct hinterv *hit;
 145         struct ray locray = *ray;
 146
 147         xform_ray(&locray, o->ob.inv_xform);
 148
 149         ndotr = o->plane.nx * locray.dx + o->plane.ny * locray.dy + o->plane.nz * locray.dz;
 150         if(fabs(ndotr) < EPSILON) return 0;
 151
 152         vx = o->plane.nx * o->plane.d - locray.x;
 153         vy = o->plane.ny * o->plane.d - locray.y;
 154         vz = o->plane.nz * o->plane.d - locray.z;
 155
 156         ndotv = o->plane.nx * vx + o->plane.ny * vy + o->plane.nz * vz;
 157
 158         t = ndotv / ndotr;
 159         if(t < EPSILON) {
 160                 return 0;
 161         }
 162
 163         hit = alloc_hits(1);
 164         hit->o = hit->end[0].o = hit->end[1].o = o;
 165         hit->end[0].t = t;
 166         hit->end[0].x = ray->x + ray->dx * t;
 167         hit->end[0].y = ray->y + ray->dy * t;
 168         hit->end[0].z = ray->z + ray->dz * t;
 169
 170         hit->end[0].nx = hit->end[1].nx = o->plane.nx;
 171         hit->end[0].ny = hit->end[1].ny = o->plane.ny;
 172         hit->end[0].nz = hit->end[1].nz = o->plane.nz;
 173
 174         hit->end[1].t = FLT_MAX;
 175         hit->end[1].x = ray->x + ray->dx * 10000.0f;
 176         hit->end[1].y = ray->y + ray->dy * 10000.0f;
 177         hit->end[1].z = ray->z + ray->dz * 10000.0f;
 178         return hit;
 179 }
 180
 181 struct hinterv *ray_csg_un(struct ray *ray, csg_object *o)
 182 {
 183         struct hinterv *hita, *hitb, *res;
 184
 185         hita = ray_intersect(ray, o->un.a);
 186         hitb = ray_intersect(ray, o->un.b);
 187
 188         if(!hita) return hitb;
 189         if(!hitb) return hita;
 190
 191         res = interval_union(hita, hitb);
 192         free_hit_list(hita);
 193         free_hit_list(hitb);
 194         return res;
 195 }
 196
 197 struct hinterv *ray_csg_isect(struct ray *ray, csg_object *o)
 198 {
 199         struct hinterv *hita, *hitb, *res;
 200
 201         hita = ray_intersect(ray, o->isect.a);
 202         hitb = ray_intersect(ray, o->isect.b);
 203
 204         if(!hita || !hitb) {
 205                 free_hit_list(hita);
 206                 free_hit_list(hitb);
 207                 return 0;
 208         }
 209
 210         res = interval_isect(hita, hitb);
 211         free_hit_list(hita);
 212         free_hit_list(hitb);
 213         return res;
 214 }
 215
 216 struct hinterv *ray_csg_sub(struct ray *ray, csg_object *o)
 217 {
 218         struct hinterv *hita, *hitb, *res;
 219
 220         hita = ray_intersect(ray, o->un.a);
 221         hitb = ray_intersect(ray, o->un.b);
 222
 223         if(!hita) return 0;
 224         if(!hitb) return hita;
 225
 226         res = interval_sub(hita, hitb);
 227         free_hit_list(hita);
 228         free_hit_list(hitb);
 229         return res;
 230 }
 231
 232
 233 void xform_ray(struct ray *ray, float *mat)
 234 {
 235         float m3x3[16];
 236
 237         mat4_copy(m3x3, mat);
 238         mat4_upper3x3(m3x3);
 239
 240         mat4_xform3(&ray->x, mat, &ray->x);
 241         mat4_xform3(&ray->dx, m3x3, &ray->dx);
 242 }
 243
 244 static void flip_hit(struct hit *hit)
 245 {
 246         hit->nx = -hit->nx;
 247         hit->ny = -hit->ny;
 248         hit->nz = -hit->nz;
 249 }
 250
 251 static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b)
 252 {
 253         struct hinterv *res, *res2;
 254
 255         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 256                 /* disjoint */
 257                 res = alloc_hits(2);
 258                 res2 = res->next;
 259
 260                 if(a->end[0].t < b->end[0].t) {
 261                         *res = *a;
 262                         *res2 = *b;
 263                 } else {
 264                         *res = *b;
 265                         *res2 = *a;
 266                 }
 267                 res->next = res2;
 268                 res2->next = 0;
 269                 return res;
 270         }
 271
 272         res = alloc_hits(1);
 273         res->end[0] = a->end[0].t <= b->end[0].t ? a->end[0] : b->end[0];
 274         res->end[1] = a->end[1].t >= b->end[1].t ? a->end[1] : b->end[1];
 275         return res;
 276 }
 277
 278 static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b)
 279 {
 280         struct hinterv *res;
 281
 282         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 283                 /* disjoint */
 284                 return 0;
 285         }
 286
 287         res = alloc_hits(1);
 288
 289         if(a->end[0].t <= b->end[0].t && a->end[1].t >= b->end[1].t) {
 290                 /* B in A */
 291                 *res = *a;
 292                 res->next = 0;
 293                 return res;
 294         }
 295         if(a->end[0].t > b->end[0].t && a->end[1].t < b->end[1].t) {
 296                 /* A in B */
 297                 *res = *b;
 298                 res->next = 0;
 299                 return res;
 300         }
 301
 302         /* partial overlap */
 303         if(a->end[0].t < b->end[0].t) {
 304                 res->end[0] = b->end[0];
 305                 res->end[1] = a->end[1];
 306         } else {
 307                 res->end[0] = a->end[0];
 308                 res->end[1] = b->end[1];
 309         }
 310         return res;
 311 }
 312
 313 static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b)
 314 {
 315         struct hinterv *res;
 316
 317         if(a->end[0].t >= b->end[0].t && a->end[1].t <= b->end[1].t) {
 318                 /* A in B */
 319                 return 0;
 320         }
 321
 322         if(a->end[0].t < b->end[0].t && a->end[1].t > b->end[1].t) {
 323                 /* B in A */
 324                 res = alloc_hits(2);
 325                 res->end[0] = a->end[0];
 326                 res->end[1] = b->end[0];
 327                 res->next->end[0] = b->end[1];
 328                 res->next->end[1] = a->end[1];
 329                 return res;
 330         }
 331
 332         res = alloc_hits(1);
 333
 334         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 335                 /* disjoint */
 336                 *res = *a;
 337                 res->next = 0;
 338                 return res;
 339         }
 340
 341         /* partial overlap */
 342         if(a->end[0].t <= b->end[0].t) {
 343                 res->end[0] = a->end[0];
 344                 res->end[1] = b->end[0];
 345         } else {
 346                 res->end[0] = b->end[1];
 347                 res->end[1] = a->end[1];
 348         }
 349
 350         flip_hit(res->end + 0);
 351         flip_hit(res->end + 1);
 352         return res;
 353 }