2 /*
3 this code written by Forest Hale, on 2004-10-17, and placed into public domain
4 this implements Quadratic BSpline surfaces as seen in Quake3 by id Software
6 a small rant on misuse of the name 'bezier': many people seem to think that
7 bezier is a generic term for splines, but it is not, it is a term for a
8 specific type of bspline (4 control points, cubic bspline), bsplines are the
9 generalization of the bezier spline to support dimensions other than cubic.
11 example equations for 1-5 control point bsplines being sampled as t=0...1
12 1: flat (0th dimension)
13 o = a
14 2: linear (1st dimension)
15 o = a * (1 - t) + b * t
16 3: quadratic bspline (2nd dimension)
17 o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
18 4: cubic (bezier) bspline (3rd dimension)
19 o = a * (1 - t) * (1 - t) * (1 - t) + 3 * b * (1 - t) * (1 - t) * t + 3 * c * (1 - t) * t * t + d * t * t * t
20 5: quartic bspline (4th dimension)
21 o = a * (1 - t) * (1 - t) * (1 - t) * (1 - t) + 4 * b * (1 - t) * (1 - t) * (1 - t) * t + 6 * c * (1 - t) * (1 - t) * t * t + 4 * d * (1 - t) * t * t * t + e * t * t * t * t
23 arbitrary dimension bspline
24 double factorial(int n)
25 {
26         int i;
27         double f;
28         f = 1;
29         for (i = 1;i < n;i++)
30                 f = f * i;
31         return f;
32 }
33 double bsplinesample(int dimensions, double t, double *param)
34 {
35         double o = 0;
36         for (i = 0;i < dimensions + 1;i++)
37                 o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
38         return o;
39 }
40 */
42 #include <math.h>
43 #include "curves.h"
45 // usage:
46 // to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
47 // Q3PatchSubdivideFloat(3, sizeof(float), outvertices, 5, 5, sizeof(float), patchvertices, 4, 4);
48 void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
49 {
50         int k, l, x, y, component, outputwidth = (patchwidth-1)*tesselationwidth+1;
51         float px, py, *v, a, b, c, *cp, temp;
52         // iterate over the individual 3x3 quadratic spline surfaces one at a time
53         // expanding them to fill the output array (with some overlap to ensure
54         // the edges are filled)
55         for (k = 0;k < patchheight-1;k += 2)
56         {
57                 for (l = 0;l < patchwidth-1;l += 2)
58                 {
59                         // set up control point pointers for quicker lookup later
60                         for (y = 0;y < 3;y++)
61                                 for (x = 0;x < 3;x++)
62                                         cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
63                         // for each row...
64                         for (y = 0;y <= tesselationheight*2;y++)
65                         {
66                                 // calculate control points for this row by collapsing the 3
67                                 // rows of control points to one row using py
68                                 py = (float)y / (float)(tesselationheight*2);
69                                 // calculate quadratic spline weights for py
70                                 a = ((1.0f - py) * (1.0f - py));
71                                 b = ((1.0f - py) * (2.0f * py));
72                                 c = ((       py) * (       py));
73                                 for (component = 0;component < numcomponents;component++)
74                                 {
75                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
76                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
77                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
78                                 }
79                                 // fetch a pointer to the beginning of the output vertex row
80                                 v = (float *)((unsigned char *)outputvertices + ((k * tesselationheight + y) * outputwidth + l * tesselationwidth) * outputstride);
81                                 // for each column of the row...
82                                 for (x = 0;x <= (tesselationwidth*2);x++)
83                                 {
84                                         // calculate point based on the row control points
85                                         px = (float)x / (float)(tesselationwidth*2);
86                                         // calculate quadratic spline weights for px
87                                         // (could be precalculated)
88                                         a = ((1.0f - px) * (1.0f - px));
89                                         b = ((1.0f - px) * (2.0f * px));
90                                         c = ((       px) * (       px));
91                                         for (component = 0;component < numcomponents;component++)
92                                                 v[component] = temp[component] * a + temp[component] * b + temp[component] * c;
93                                         // advance to next output vertex using outputstride
94                                         // (the next vertex may not be directly following this
95                                         // one, as this may be part of a larger structure)
96                                         v = (float *)((unsigned char *)v + outputstride);
97                                 }
98                         }
99                 }
100         }
101 #if 0
102         // enable this if you want results printed out
103         printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
104         for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
105         {
106                 for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
107                 {
108                         printf("(");
109                         for (component = 0;component < numcomponents;component++)
110                                 printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
111                         printf(") ");
112                 }
113                 printf("\n");
114         }
115         printf("}\n");
116 #endif
117 }
119 // returns how much tesselation of each segment is needed to remain under tolerance
120 int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
121 {
122         int c, x, y;
123         const float *patch;
124         float deviation, squareddeviation, bestsquareddeviation;
125         bestsquareddeviation = 0;
126         for (y = 0;y < patchheight;y++)
127         {
128                 for (x = 0;x < patchwidth-1;x += 2)
129                 {
130                         squareddeviation = 0;
131                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
132                         {
133                                 deviation = patch[components] * 0.5f - patch * 0.25f - patch[2*components] * 0.25f;
134                                 squareddeviation += deviation*deviation;
135                         }
136                         if (bestsquareddeviation < squareddeviation)
137                                 bestsquareddeviation = squareddeviation;
138                 }
139         }
140         return (int)floor(log(sqrt(bestsquareddeviation) / tolerance) / log(2)) + 1;
141 }
143 // returns how much tesselation of each segment is needed to remain under tolerance
144 int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
145 {
146         int c, x, y;
147         const float *patch;
148         float deviation, squareddeviation, bestsquareddeviation;
149         bestsquareddeviation = 0;
150         for (y = 0;y < patchheight-1;y += 2)
151         {
152                 for (x = 0;x < patchwidth;x++)
153                 {
154                         squareddeviation = 0;
155                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
156                         {
157                                 deviation = patch[patchwidth*components] * 0.5f - patch * 0.25f - patch[2*patchwidth*components] * 0.25f;
158                                 squareddeviation += deviation*deviation;
159                         }
160                         if (bestsquareddeviation < squareddeviation)
161                                 bestsquareddeviation = squareddeviation;
162                 }
163         }
164         return (int)floor(log(sqrt(bestsquareddeviation) / tolerance) / log(2)) + 1;
165 }
167 // calculates elements for a grid of vertices
168 // (such as those produced by Q3PatchTesselate)
169 // (note: width and height are the actual vertex size, this produces
170 //  (width-1)*(height-1)*2 triangles, 3 elements each)
171 void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
172 {
173         int x, y, row0, row1;
174         for (y = 0;y < height - 1;y++)
175         {
176                 row0 = firstvertex + (y + 0) * width;
177                 row1 = firstvertex + (y + 1) * width;
178                 for (x = 0;x < width - 1;x++)
179                 {
180                         *elements++ = row0;
181                         *elements++ = row1;
182                         *elements++ = row0 + 1;
183                         *elements++ = row1;
184                         *elements++ = row1 + 1;
185                         *elements++ = row0 + 1;
186                         row0++;
187                         row1++;
188                 }
189         }
190 }