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 <assert.h>
   6 #include "geom.h"
   7 #include "matrix.h"
   8
   9 #define EPSILON         1e-6f
  10
  11 static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b);
  12 static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b);
  13 static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b);
  14
  15 /* TODO custom hit allocator */
  16 struct hinterv *alloc_hit(void)
  17 {
  18         struct hinterv *hit = calloc(sizeof *hit, 1);
  19         if(!hit) {
  20                 perror("failed to allocate ray hit node");
  21                 abort();
  22         }
  23         return hit;
  24 }
  25
  26 struct hinterv *alloc_hits(int n)
  27 {
  28         int i;
  29         struct hinterv *list = 0;
  30
  31         for(i=0; i<n; i++) {
  32                 struct hinterv *hit = alloc_hit();
  33                 hit->next = list;
  34                 list = hit;
  35         }
  36         return list;
  37 }
  38
  39
  40 void free_hit(struct hinterv *hit)
  41 {
  42         free(hit);
  43 }
  44
  45 void free_hit_list(struct hinterv *hit)
  46 {
  47         while(hit) {
  48                 struct hinterv *tmp = hit;
  49                 hit = hit->next;
  50                 free_hit(tmp);
  51         }
  52 }
  53
  54 struct hinterv *ray_intersect(struct ray *ray, csg_object *o)
  55 {
  56         switch(o->ob.type) {
  57         case OB_SPHERE:
  58                 return ray_sphere(ray, o);
  59         case OB_CYLINDER:
  60                 return ray_cylinder(ray, o);
  61         case OB_PLANE:
  62                 return ray_plane(ray, o);
  63         case OB_BOX:
  64                 return ray_box(ray, o);
  65         case OB_UNION:
  66                 return ray_csg_un(ray, o);
  67         case OB_INTERSECTION:
  68                 return ray_csg_isect(ray, o);
  69         case OB_SUBTRACTION:
  70                 return ray_csg_sub(ray, o);
  71         default:
  72                 break;
  73         }
  74         return 0;
  75 }
  76
  77 struct hinterv *ray_sphere(struct ray *ray, csg_object *o)
  78 {
  79         int i;
  80         float a, b, c, d, sqrt_d, t[2], sq_rad, tmp;
  81         struct hinterv *hit;
  82         struct ray locray = *ray;
  83
  84         if(o->sph.rad == 0.0f) {
  85                 return 0;
  86         }
  87         sq_rad = o->sph.rad * o->sph.rad;
  88
  89         xform_ray(&locray, o->ob.inv_xform);
  90
  91         a = locray.dx * locray.dx + locray.dy * locray.dy + locray.dz * locray.dz;
  92         b = 2.0f * (locray.dx * locray.x + locray.dy * locray.y + locray.dz * locray.z);
  93         c = (locray.x * locray.x + locray.y * locray.y + locray.z * locray.z) - sq_rad;
  94
  95         d = b * b - 4.0f * a * c;
  96         if(d < EPSILON) return 0;
  97
  98         sqrt_d = sqrt(d);
  99         t[0] = (-b + sqrt_d) / (2.0f * a);
 100         t[1] = (-b - sqrt_d) / (2.0f * a);
 101
 102         if(t[0] < EPSILON && t[1] < EPSILON) {
 103                 return 0;
 104         }
 105
 106         if(t[1] < t[0]) {
 107                 tmp = t[0];
 108                 t[0] = t[1];
 109                 t[1] = tmp;
 110         }
 111
 112         hit = alloc_hits(1);
 113         hit->o = o;
 114         for(i=0; i<2; i++) {
 115                 float c[3] = {0, 0, 0};
 116                 float x, y, z;
 117
 118                 mat4_xform3(c, o->ob.xform, c);
 119
 120                 x = ray->x + ray->dx * t[i];
 121                 y = ray->y + ray->dy * t[i];
 122                 z = ray->z + ray->dz * t[i];
 123
 124                 hit->end[i].t = t[i];
 125                 hit->end[i].x = x;
 126                 hit->end[i].y = y;
 127                 hit->end[i].z = z;
 128                 hit->end[i].nx = (x - c[0]) / o->sph.rad;
 129                 hit->end[i].ny = (y - c[1]) / o->sph.rad;
 130                 hit->end[i].nz = (z - c[2]) / o->sph.rad;
 131                 hit->end[i].o = o;
 132         }
 133         return hit;
 134 }
 135
 136 static int ray_cylcap(struct ray *ray, float y, float rad, float *tres)
 137 {
 138         float ndotr, ndotv, vy, t;
 139         float ny = y > 0.0f ? 1.0f : -1.0f;
 140         float x, z, lensq;
 141
 142         ndotr = ny * ray->dy;
 143         if(fabs(ndotr) < EPSILON) return 0;
 144
 145         vy = y - ray->y;
 146
 147         ndotv = ny * vy;
 148
 149         t = ndotv / ndotr;
 150
 151         x = ray->x + ray->dx * t;
 152         z = ray->z + ray->dz * t;
 153         lensq = x * x + z * z;
 154
 155         if(lensq <= rad * rad) {
 156                 *tres = t;
 157                 return 1;
 158         }
 159         return 0;
 160 }
 161
 162 struct hinterv *ray_cylinder(struct ray *ray, csg_object *o)
 163 {
 164         int i, out[2] = {0}, t_is_cap[2] = {0};
 165         float a, b, c, d, sqrt_d, t[2], sq_rad, tmp, y[2], hh, cap_t;
 166         struct hinterv *hit;
 167         struct ray locray = *ray;
 168
 169         if(o->cyl.rad == 0.0f || o->cyl.height == 0.0f) {
 170                 return 0;
 171         }
 172         sq_rad = o->cyl.rad * o->cyl.rad;
 173         hh = o->cyl.height / 2.0f;
 174
 175         xform_ray(&locray, o->ob.inv_xform);
 176
 177         a = locray.dx * locray.dx + locray.dz * locray.dz;
 178         b = 2.0f * (locray.dx * locray.x + locray.dz * locray.z);
 179         c = locray.x * locray.x + locray.z * locray.z - sq_rad;
 180
 181         d = b * b - 4.0f * a * c;
 182         if(d < EPSILON) return 0;
 183
 184         sqrt_d = sqrt(d);
 185         t[0] = (-b + sqrt_d) / (2.0f * a);
 186         t[1] = (-b - sqrt_d) / (2.0f * a);
 187
 188         if(t[0] < EPSILON && t[1] < EPSILON) {
 189                 return 0;
 190         }
 191         if(t[1] < t[0]) {
 192                 tmp = t[0];
 193                 t[0] = t[1];
 194                 t[1] = tmp;
 195         }
 196
 197         y[0] = locray.y + locray.dy * t[0];
 198         y[1] = locray.y + locray.dy * t[1];
 199
 200         if(y[0] < -hh || y[0] > hh) {
 201                 out[0] = 1;
 202         }
 203         if(y[1] < -hh || y[1] > hh) {
 204                 out[1] = 1;
 205         }
 206
 207         if(out[0]) {
 208                 t[0] = t[1];
 209         }
 210         if(out[1]) {
 211                 t[1] = t[0];
 212         }
 213
 214         if(ray_cylcap(ray, hh, o->cyl.rad, &cap_t)) {
 215                 if(cap_t < t[0]) {
 216                         t[0] = cap_t;
 217                         t_is_cap[0] = 1;
 218                         out[0] = 0;
 219                 }
 220                 if(cap_t > t[1]) {
 221                         t[1] = cap_t;
 222                         t_is_cap[1] = 1;
 223                         out[1] = 0;
 224                 }
 225         }
 226         if(ray_cylcap(ray, -hh, o->cyl.rad, &cap_t)) {
 227                 if(cap_t < t[0]) {
 228                         t[0] = cap_t;
 229                         t_is_cap[0] = -1;
 230                         out[0] = 0;
 231                 }
 232                 if(cap_t > t[1]) {
 233                         t[1] = cap_t;
 234                         t_is_cap[1] = -1;
 235                         out[1] = 0;
 236                 }
 237         }
 238
 239         if(out[0] && out[1]) {
 240                 return 0;
 241         }
 242
 243         hit = alloc_hits(1);
 244         hit->o = o;
 245         for(i=0; i<2; i++) {
 246                 float c[3] = {0, 0, 0};
 247                 float x, y, z;
 248
 249                 x = ray->x + ray->dx * t[i];
 250                 y = ray->y + ray->dy * t[i];
 251                 z = ray->z + ray->dz * t[i];
 252
 253                 if(t_is_cap[i]) {
 254                         hit->end[i].nx = hit->end[i].nz = 0.0f;
 255                         hit->end[i].ny = t_is_cap[i] > 0 ? 1.0f : -1.0f;
 256                 } else {
 257                         c[1] = locray.y + locray.dy * t[i];
 258                         mat4_xform3(c, o->ob.xform, c);
 259
 260                         hit->end[i].nx = (x - c[0]) / o->cyl.rad;
 261                         hit->end[i].ny = (y - c[1]) / o->cyl.rad;
 262                         hit->end[i].nz = (z - c[2]) / o->cyl.rad;
 263                 }
 264
 265                 hit->end[i].t = t[i];
 266                 hit->end[i].x = x;
 267                 hit->end[i].y = y;
 268                 hit->end[i].z = z;
 269                 hit->end[i].o = o;
 270         }
 271         return hit;
 272 }
 273
 274 struct hinterv *ray_plane(struct ray *ray, csg_object *o)
 275 {
 276         float vx, vy, vz, ndotv, ndotr, t;
 277         struct hinterv *hit;
 278         struct ray locray = *ray;
 279
 280         xform_ray(&locray, o->ob.inv_xform);
 281
 282         ndotr = o->plane.nx * locray.dx + o->plane.ny * locray.dy + o->plane.nz * locray.dz;
 283         if(fabs(ndotr) < EPSILON) return 0;
 284
 285         vx = o->plane.nx * o->plane.d - locray.x;
 286         vy = o->plane.ny * o->plane.d - locray.y;
 287         vz = o->plane.nz * o->plane.d - locray.z;
 288
 289         ndotv = o->plane.nx * vx + o->plane.ny * vy + o->plane.nz * vz;
 290
 291         t = ndotv / ndotr;
 292         if(t < EPSILON) {
 293                 return 0;
 294         }
 295
 296         hit = alloc_hits(1);
 297         hit->o = hit->end[0].o = hit->end[1].o = o;
 298         hit->end[0].t = t;
 299         hit->end[0].x = ray->x + ray->dx * t;
 300         hit->end[0].y = ray->y + ray->dy * t;
 301         hit->end[0].z = ray->z + ray->dz * t;
 302
 303         hit->end[0].nx = hit->end[1].nx = o->plane.nx;
 304         hit->end[0].ny = hit->end[1].ny = o->plane.ny;
 305         hit->end[0].nz = hit->end[1].nz = o->plane.nz;
 306
 307         hit->end[1].t = FLT_MAX;
 308         hit->end[1].x = ray->x + ray->dx * 10000.0f;
 309         hit->end[1].y = ray->y + ray->dy * 10000.0f;
 310         hit->end[1].z = ray->z + ray->dz * 10000.0f;
 311         return hit;
 312 }
 313
 314 #define BEXT(x) ((x) * 0.49999)
 315
 316 struct hinterv *ray_box(struct ray *ray, csg_object *o)
 317 {
 318         int i, sign[3];
 319         float param[2][3];
 320         float inv_dir[3];
 321         float tmin, tmax, tymin, tymax, tzmin, tzmax;
 322         struct hinterv *hit;
 323         struct ray locray = *ray;
 324         float dirmat[16];
 325
 326         xform_ray(&locray, o->ob.inv_xform);
 327
 328         for(i=0; i<3; i++) {
 329                 float sz = *(&o->box.xsz + i);
 330                 param[0][i] = -0.5 * sz;
 331                 param[1][i] = 0.5 * sz;
 332
 333                 inv_dir[i] = 1.0f / *(&locray.dx + i);
 334                 sign[i] = inv_dir[i] < 0;
 335         }
 336
 337         tmin = (param[sign[0]][0] - locray.x) * inv_dir[0];
 338         tmax = (param[1 - sign[0]][0] - locray.x) * inv_dir[0];
 339         tymin = (param[sign[1]][1] - locray.y) * inv_dir[1];
 340         tymax = (param[1 - sign[1]][1] - locray.y) * inv_dir[1];
 341
 342         if(tmin > tymax || tymin > tmax) {
 343                 return 0;
 344         }
 345         if(tymin > tmin) {
 346                 tmin = tymin;
 347         }
 348         if(tymax < tmax) {
 349                 tmax = tymax;
 350         }
 351
 352         tzmin = (param[sign[2]][2] - locray.z) * inv_dir[2];
 353         tzmax = (param[1 - sign[2]][2] - locray.z) * inv_dir[2];
 354
 355         if(tmin > tzmax || tzmin > tmax) {
 356                 return 0;
 357         }
 358         if(tzmin > tmin) {
 359                 tmin = tzmin;
 360         }
 361         if(tzmax < tmax) {
 362                 tmax = tzmax;
 363         }
 364
 365         mat4_copy(dirmat, o->ob.xform);
 366         mat4_upper3x3(dirmat);
 367
 368         hit = alloc_hits(1);
 369         hit->o = o;
 370         for(i=0; i<2; i++) {
 371                 float n[3] = {0};
 372                 float t = i == 0 ? tmin : tmax;
 373
 374                 float x = (locray.x + locray.dx * t) / o->box.xsz;
 375                 float y = (locray.y + locray.dy * t) / o->box.ysz;
 376                 float z = (locray.z + locray.dz * t) / o->box.zsz;
 377
 378                 if(fabs(x) > fabs(y) && fabs(x) > fabs(z)) {
 379                         n[0] = x > 0.0f ? 1.0f : -1.0f;
 380                 } else if(fabs(y) > fabs(z)) {
 381                         n[1] = y > 0.0f ? 1.0f : -1.0f;
 382                 } else {
 383                         n[2] = z > 0.0f ? 1.0f : -1.0f;
 384                 }
 385
 386                 hit->end[i].o = o;
 387                 hit->end[i].t = t;
 388                 hit->end[i].x = ray->x + ray->dx * t;
 389                 hit->end[i].y = ray->y + ray->dy * t;
 390                 hit->end[i].z = ray->z + ray->dz * t;
 391                 mat4_xform3(&hit->end[i].nx, dirmat, n);
 392         }
 393         return hit;
 394 }
 395
 396 struct hinterv *ray_csg_un(struct ray *ray, csg_object *o)
 397 {
 398         struct hinterv *hita, *hitb, *res;
 399
 400         hita = ray_intersect(ray, o->un.a);
 401         hitb = ray_intersect(ray, o->un.b);
 402
 403         if(!hita) return hitb;
 404         if(!hitb) return hita;
 405
 406         res = interval_union(hita, hitb);
 407         free_hit_list(hita);
 408         free_hit_list(hitb);
 409         return res;
 410 }
 411
 412 struct hinterv *ray_csg_isect(struct ray *ray, csg_object *o)
 413 {
 414         struct hinterv *hita, *hitb, *res;
 415
 416         hita = ray_intersect(ray, o->isect.a);
 417         hitb = ray_intersect(ray, o->isect.b);
 418
 419         if(!hita || !hitb) {
 420                 free_hit_list(hita);
 421                 free_hit_list(hitb);
 422                 return 0;
 423         }
 424
 425         res = interval_isect(hita, hitb);
 426         free_hit_list(hita);
 427         free_hit_list(hitb);
 428         return res;
 429 }
 430
 431 struct hinterv *ray_csg_sub(struct ray *ray, csg_object *o)
 432 {
 433         struct hinterv *hita, *hitb, *res;
 434
 435         hita = ray_intersect(ray, o->un.a);
 436         hitb = ray_intersect(ray, o->un.b);
 437
 438         if(!hita) return 0;
 439         if(!hitb) return hita;
 440
 441         res = interval_sub(hita, hitb);
 442         free_hit_list(hita);
 443         free_hit_list(hitb);
 444         return res;
 445 }
 446
 447
 448 void xform_ray(struct ray *ray, float *mat)
 449 {
 450         float m3x3[16];
 451
 452         mat4_copy(m3x3, mat);
 453         mat4_upper3x3(m3x3);
 454
 455         mat4_xform3(&ray->x, mat, &ray->x);
 456         mat4_xform3(&ray->dx, m3x3, &ray->dx);
 457 }
 458
 459 static void flip_hit(struct hit *hit)
 460 {
 461         hit->nx = -hit->nx;
 462         hit->ny = -hit->ny;
 463         hit->nz = -hit->nz;
 464 }
 465
 466 static struct hinterv *interval_union(struct hinterv *a, struct hinterv *b)
 467 {
 468         struct hinterv *res, *res2;
 469
 470         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 471                 /* disjoint */
 472                 res = alloc_hits(2);
 473                 res2 = res->next;
 474
 475                 if(a->end[0].t < b->end[0].t) {
 476                         *res = *a;
 477                         *res2 = *b;
 478                 } else {
 479                         *res = *b;
 480                         *res2 = *a;
 481                 }
 482                 res->next = res2;
 483                 res2->next = 0;
 484                 return res;
 485         }
 486
 487         res = alloc_hits(1);
 488         res->end[0] = a->end[0].t <= b->end[0].t ? a->end[0] : b->end[0];
 489         res->end[1] = a->end[1].t >= b->end[1].t ? a->end[1] : b->end[1];
 490         return res;
 491 }
 492
 493 static struct hinterv *interval_isect(struct hinterv *a, struct hinterv *b)
 494 {
 495         struct hinterv *res;
 496
 497         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 498                 /* disjoint */
 499                 return 0;
 500         }
 501
 502         res = alloc_hits(1);
 503
 504         if(a->end[0].t <= b->end[0].t && a->end[1].t >= b->end[1].t) {
 505                 /* B in A */
 506                 *res = *a;
 507                 res->next = 0;
 508                 return res;
 509         }
 510         if(a->end[0].t > b->end[0].t && a->end[1].t < b->end[1].t) {
 511                 /* A in B */
 512                 *res = *b;
 513                 res->next = 0;
 514                 return res;
 515         }
 516
 517         /* partial overlap */
 518         if(a->end[0].t < b->end[0].t) {
 519                 res->end[0] = b->end[0];
 520                 res->end[1] = a->end[1];
 521         } else {
 522                 res->end[0] = a->end[0];
 523                 res->end[1] = b->end[1];
 524         }
 525         return res;
 526 }
 527
 528 static struct hinterv *interval_sub(struct hinterv *a, struct hinterv *b)
 529 {
 530         struct hinterv *res;
 531
 532         if(a->end[0].t >= b->end[0].t && a->end[1].t <= b->end[1].t) {
 533                 /* A in B */
 534                 return 0;
 535         }
 536
 537         if(a->end[0].t < b->end[0].t && a->end[1].t > b->end[1].t) {
 538                 /* B in A */
 539                 res = alloc_hits(2);
 540                 res->end[0] = a->end[0];
 541                 res->end[1] = b->end[0];
 542                 res->next->end[0] = b->end[1];
 543                 res->next->end[1] = a->end[1];
 544                 return res;
 545         }
 546
 547         res = alloc_hits(1);
 548
 549         if(a->end[0].t > b->end[1].t || a->end[1].t < b->end[0].t) {
 550                 /* disjoint */
 551                 *res = *a;
 552                 res->next = 0;
 553                 return res;
 554         }
 555
 556         /* partial overlap */
 557         if(a->end[0].t <= b->end[0].t) {
 558                 res->end[0] = a->end[0];
 559                 res->end[1] = b->end[0];
 560         } else {
 561                 res->end[0] = b->end[1];
 562                 res->end[1] = a->end[1];
 563         }
 564
 565         flip_hit(res->end + 0);
 566         flip_hit(res->end + 1);
 567         return res;
 568 }