1 /*
2 * quadric.c - This file contains the functions for dealing with quadrics.
3 *
4 * $Id: quadric.c,v 1.24 2007/02/04 05:38:08 johns Exp $
5 */
6
7 #include "machine.h"
8 #include "types.h"
9 #include "macros.h"
10 #include "quadric.h"
11 #include "vector.h"
12 #include "intersect.h"
13 #include "util.h"
14
quadric_bbox(void * obj,vector * min,vector * max)15 int quadric_bbox(void * obj, vector * min, vector * max) {
16 return 0;
17 }
18
19 static object_methods quadric_methods = {
20 (void (*)(const void *, void *))(quadric_intersect),
21 (void (*)(const void *, const void *, const void *, void *))(quadric_normal),
22 quadric_bbox,
23 free
24 };
25
newquadric(void)26 quadric * newquadric(void) {
27 quadric * q;
28
29 q=(quadric *) malloc(sizeof(quadric));
30 memset(q, 0, sizeof(quadric));
31 q->ctr.x=0.0;
32 q->ctr.y=0.0;
33 q->ctr.z=0.0;
34 q->methods = &quadric_methods;
35
36 return q;
37 }
38
quadric_intersect(const quadric * q,ray * ry)39 void quadric_intersect(const quadric * q, ray * ry) {
40 flt Aq, Bq, Cq;
41 flt t1, t2;
42 flt disc;
43 vector rd;
44 vector ro;
45
46 rd=ry->d;
47 VNorm(&rd);
48
49 ro.x = ry->o.x - q->ctr.x;
50 ro.y = ry->o.y - q->ctr.y;
51 ro.z = ry->o.z - q->ctr.z;
52
53
54 Aq = (q->mat.a*(rd.x * rd.x)) +
55 (2.0 * q->mat.b * rd.x * rd.y) +
56 (2.0 * q->mat.c * rd.x * rd.z) +
57 (q->mat.e * (rd.y * rd.y)) +
58 (2.0 * q->mat.f * rd.y * rd.z) +
59 (q->mat.h * (rd.z * rd.z));
60
61 Bq = 2.0 * (
62 (q->mat.a * ro.x * rd.x) +
63 (q->mat.b * ((ro.x * rd.y) + (rd.x * ro.y))) +
64 (q->mat.c * ((ro.x * rd.z) + (rd.x * ro.z))) +
65 (q->mat.d * rd.x) +
66 (q->mat.e * ro.y * rd.y) +
67 (q->mat.f * ((ro.y * rd.z) + (rd.y * ro.z))) +
68 (q->mat.g * rd.y) +
69 (q->mat.h * ro.z * rd.z) +
70 (q->mat.i * rd.z)
71 );
72
73 Cq = (q->mat.a * (ro.x * ro.x)) +
74 (2.0 * q->mat.b * ro.x * ro.y) +
75 (2.0 * q->mat.c * ro.x * ro.z) +
76 (2.0 * q->mat.d * ro.x) +
77 (q->mat.e * (ro.y * ro.y)) +
78 (2.0 * q->mat.f * ro.y * ro.z) +
79 (2.0 * q->mat.g * ro.y) +
80 (q->mat.h * (ro.z * ro.z)) +
81 (2.0 * q->mat.i * ro.z) +
82 q->mat.j;
83
84 if (Aq == 0.0) {
85 t1 = - Cq / Bq;
86 ry->add_intersection(t1, (object *) q, ry);
87 }
88 else {
89 disc=(Bq*Bq - 4.0 * Aq * Cq);
90 if (disc > 0.0) {
91 disc=sqrt(disc);
92 t1 = (-Bq + disc) / (2.0 * Aq);
93 t2 = (-Bq - disc) / (2.0 * Aq);
94 ry->add_intersection(t1, (object *) q, ry);
95 ry->add_intersection(t2, (object *) q, ry);
96 }
97 }
98 }
99
quadric_normal(const quadric * q,const vector * pnt,const ray * incident,vector * N)100 void quadric_normal(const quadric * q, const vector * pnt, const ray * incident, vector * N) {
101 flt invlen;
102
103 N->x = (q->mat.a*(pnt->x - q->ctr.x) +
104 q->mat.b*(pnt->y - q->ctr.y) +
105 q->mat.c*(pnt->z - q->ctr.z) + q->mat.d);
106
107 N->y = (q->mat.b*(pnt->x - q->ctr.x) +
108 q->mat.e*(pnt->y - q->ctr.y) +
109 q->mat.f*(pnt->z - q->ctr.z) + q->mat.g);
110
111 N->z = (q->mat.c*(pnt->x - q->ctr.x) +
112 q->mat.f*(pnt->y - q->ctr.y) +
113 q->mat.h*(pnt->z - q->ctr.z) + q->mat.i);
114
115 invlen = 1.0 / sqrt(N->x*N->x + N->y*N->y + N->z*N->z);
116 N->x *= invlen;
117 N->y *= invlen;
118 N->z *= invlen;
119
120 /* Flip surface normal to point toward the viewer if necessary */
121 if (VDot(N, &(incident->d)) > 0.0) {
122 N->x=-N->x;
123 N->y=-N->y;
124 N->z=-N->z;
125 }
126 }
127
128
129