1 /*
2 this code written by Ashley Rose Hale (LadyHavoc), on 2004-10-17, and placed into public domain
3 this implements Quadratic BSpline surfaces as seen in Quake3 by id Software
5 a small rant on misuse of the name 'bezier': many people seem to think that
6 bezier is a generic term for splines, but it is not, it is a term for a
7 specific type of bspline (4 control points, cubic bspline), bsplines are the
8 generalization of the bezier spline to support dimensions other than cubic.
10 example equations for 1-5 control point bsplines being sampled as t=0...1
11 1: flat (0th dimension)
12 o = a
13 2: linear (1st dimension)
14 o = a * (1 - t) + b * t
15 3: quadratic bspline (2nd dimension)
16 o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
17 4: cubic (bezier) bspline (3rd dimension)
18 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
19 5: quartic bspline (4th dimension)
20 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
22 arbitrary dimension bspline
23 double factorial(int n)
24 {
25         int i;
26         double f;
27         f = 1;
28         for (i = 1;i < n;i++)
29                 f = f * i;
30         return f;
31 }
32 double bsplinesample(int dimensions, double t, double *param)
33 {
34         double o = 0;
35         for (i = 0;i < dimensions + 1;i++)
36                 o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
37         return o;
38 }
39 */
41 #include "quakedef.h"
42 #include "mathlib.h"
44 #include <math.h>
45 #include "curves.h"
47 // Calculate number of resulting vertex rows/columns by given patch size and tesselation factor
48 // tess=0 means that we reduce detalization of base 3x3 patches by removing middle row and column of vertices
49 // "DimForTess" is "DIMension FOR TESSelation factor"
50 // NB: tess=0 actually means that tess must be 0.5, but obviously it can't because it is of int type. (so "a*tess"-like code is replaced by "a/2" if tess=0)
51 int Q3PatchDimForTess(int size, int tess)
52 {
53         if (tess > 0)
54                 return (size - 1) * tess + 1;
55         else if (tess == 0)
56                 return (size - 1) / 2 + 1;
57         else
58                 return 0; // Maybe warn about wrong tess here?
59 }
61 // usage:
62 // to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
63 // Q3PatchSubdivideFloat(3, sizeof(float), outvertices, 5, 5, sizeof(float), patchvertices, 4, 4);
64 void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
65 {
66         int k, l, x, y, component, outputwidth = Q3PatchDimForTess(patchwidth, tesselationwidth);
67         float px, py, *v, a, b, c, *cp, temp;
68         int xmax = max(1, 2*tesselationwidth);
69         int ymax = max(1, 2*tesselationheight);
71         // iterate over the individual 3x3 quadratic spline surfaces one at a time
72         // expanding them to fill the output array (with some overlap to ensure
73         // the edges are filled)
74         for (k = 0;k < patchheight-1;k += 2)
75         {
76                 for (l = 0;l < patchwidth-1;l += 2)
77                 {
78                         // set up control point pointers for quicker lookup later
79                         for (y = 0;y < 3;y++)
80                                 for (x = 0;x < 3;x++)
81                                         cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
82                         // for each row...
83                         for (y = 0;y <= ymax;y++)
84                         {
85                                 // calculate control points for this row by collapsing the 3
86                                 // rows of control points to one row using py
87                                 py = (float)y / (float)ymax;
88                                 // calculate quadratic spline weights for py
89                                 a = ((1.0f - py) * (1.0f - py));
90                                 b = ((1.0f - py) * (2.0f * py));
91                                 c = ((       py) * (       py));
92                                 for (component = 0;component < numcomponents;component++)
93                                 {
94                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
95                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
96                                         temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
97                                 }
98                                 // fetch a pointer to the beginning of the output vertex row
99                                 v = (float *)((unsigned char *)outputvertices + ((k * ymax / 2 + y) * outputwidth + l * xmax / 2) * outputstride);
100                                 // for each column of the row...
101                                 for (x = 0;x <= xmax;x++)
102                                 {
103                                         // calculate point based on the row control points
104                                         px = (float)x / (float)xmax;
105                                         // calculate quadratic spline weights for px
106                                         // (could be precalculated)
107                                         a = ((1.0f - px) * (1.0f - px));
108                                         b = ((1.0f - px) * (2.0f * px));
109                                         c = ((       px) * (       px));
110                                         for (component = 0;component < numcomponents;component++)
111                                                 v[component] = temp[component] * a + temp[component] * b + temp[component] * c;
112                                         // advance to next output vertex using outputstride
113                                         // (the next vertex may not be directly following this
114                                         // one, as this may be part of a larger structure)
115                                         v = (float *)((unsigned char *)v + outputstride);
116                                 }
117                         }
118                 }
119         }
120 #if 0
121         // enable this if you want results printed out
122         printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
123         for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
124         {
125                 for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
126                 {
127                         printf("(");
128                         for (component = 0;component < numcomponents;component++)
129                                 printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
130                         printf(") ");
131                 }
132                 printf("\n");
133         }
134         printf("}\n");
135 #endif
136 }
138 static int Q3PatchTesselation(float largestsquared3xcurvearea, float tolerance)
139 {
140         float f;
141         // f is actually a squared 2x curve area... so the formula had to be adjusted to give roughly the same subdivisions
142         f = pow(largestsquared3xcurvearea / 64.0f, 0.25f) / tolerance;
143         //if(f < 0.25) // VERY flat patches
144         if(f < 0.0001) // TOTALLY flat patches
145                 return 0;
146         else if(f < 2)
147                 return 1;
148         else
149                 return (int) floor(log(f) / log(2.0f)) + 1;
150                 // this is always at least 2
151                 // maps [0.25..0.5[ to -1 (actually, 1 is returned)
152                 // maps [0.5..1[ to 0 (actually, 1 is returned)
153                 // maps [1..2[ to 1
154                 // maps [2..4[ to 2
155                 // maps [4..8[ to 4
156 }
158 static float Squared3xCurveArea(const float *a, const float *control, const float *b, int components)
159 {
160 #if 0
161         // mimicing the old behaviour with the new code...
163         float deviation;
164         float quartercurvearea = 0;
165         int c;
166         for (c = 0;c < components;c++)
167         {
168                 deviation = control[c] * 0.5f - a[c] * 0.25f - b[c] * 0.25f;
169                 quartercurvearea += deviation*deviation;
170         }
172         // But as the new code now works on the squared 2x curve area, let's scale the value
173         return quartercurvearea * quartercurvearea * 64.0;
175 #else
176         // ideally, we'd like the area between the spline a->control->b and the line a->b.
177         // but as this is hard to calculate, let's calculate an upper bound of it:
178         // the area of the triangle a->control->b->a.
179         //
180         // one can prove that the area of a quadratic spline = 2/3 * the area of
181         // the triangle of its control points!
182         // to do it, first prove it for the spline through (0,0), (1,1), (2,0)
183         // (which is a parabola) and then note that moving the control point
184         // left/right is just shearing and keeps the area of both the spline and
185         // the triangle invariant.
186         //
187         // why are we going for the spline area anyway?
188         // we know that:
189         //
190         //   the area between the spline and the line a->b is a measure of the
191         //   error of approximation of the spline by the line.
192         //
193         //   also, on circle-like or parabola-like curves, you easily get that the
194         //   double amount of line approximation segments reduces the error to its quarter
195         //   (also, easy to prove for splines by doing it for one specific one, and using
196         //   affine transforms to get all other splines)
197         //
198         // so...
199         //
200         // let's calculate the area! but we have to avoid the cross product, as
201         // components is not necessarily 3
202         //
203         // the area of a triangle spanned by vectors a and b is
204         //
205         // 0.5 * |a| |b| sin gamma
206         //
207         // now, cos gamma is
208         //
209         // a.b / (|a| |b|)
210         //
211         // so the area is
212         //
213         // 0.5 * sqrt(|a|^2 |b|^2 - (a.b)^2)
214         int c;
215         float aa = 0, bb = 0, ab = 0;
216         for (c = 0;c < components;c++)
217         {
218                 float xa = a[c] - control[c];
219                 float xb = b[c] - control[c];
220                 aa += xa * xa;
221                 ab += xa * xb;
222                 bb += xb * xb;
223         }
224         // area is 0.5 * sqrt(aa*bb - ab*ab)
225         // 2x TRIANGLE area is sqrt(aa*bb - ab*ab)
226         // 3x CURVE area is sqrt(aa*bb - ab*ab)
227         return aa * bb - ab * ab;
228 #endif
229 }
231 // returns how much tesselation of each segment is needed to remain under tolerance
232 int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
233 {
234         int x, y;
235         const float *patch;
236         float squared3xcurvearea, largestsquared3xcurvearea;
237         largestsquared3xcurvearea = 0;
238         for (y = 0;y < patchheight;y++)
239         {
240                 for (x = 0;x < patchwidth-1;x += 2)
241                 {
242                         patch = in + ((y * patchwidth) + x) * components;
243                         squared3xcurvearea = Squared3xCurveArea(&patch, &patch[components], &patch[2*components], components);
244                         if (largestsquared3xcurvearea < squared3xcurvearea)
245                                 largestsquared3xcurvearea = squared3xcurvearea;
246                 }
247         }
248         return Q3PatchTesselation(largestsquared3xcurvearea, tolerance);
249 }
251 // returns how much tesselation of each segment is needed to remain under tolerance
252 int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
253 {
254         int x, y;
255         const float *patch;
256         float squared3xcurvearea, largestsquared3xcurvearea;
257         largestsquared3xcurvearea = 0;
258         for (y = 0;y < patchheight-1;y += 2)
259         {
260                 for (x = 0;x < patchwidth;x++)
261                 {
262                         patch = in + ((y * patchwidth) + x) * components;
263                         squared3xcurvearea = Squared3xCurveArea(&patch, &patch[patchwidth*components], &patch[2*patchwidth*components], components);
264                         if (largestsquared3xcurvearea < squared3xcurvearea)
265                                 largestsquared3xcurvearea = squared3xcurvearea;
266                 }
267         }
268         return Q3PatchTesselation(largestsquared3xcurvearea, tolerance);
269 }
271 // Find an equal vertex in array. Check only vertices with odd X and Y
272 static int FindEqualOddVertexInArray(int numcomponents, float *vertex, float *vertices, int width, int height)
273 {
274         int x, y, j;
275         for (y=0; y<height; y+=2)
276         {
277                 for (x=0; x<width; x+=2)
278                 {
279                         qbool found = true;
280                         for (j=0; j<numcomponents; j++)
281                                 if (fabs(*(vertex+j) - *(vertices+j)) > 0.05)
282                                 // div0: this is notably smaller than the smallest radiant grid
283                                 // but large enough so we don't need to get scared of roundoff
284                                 // errors
285                                 {
286                                         found = false;
287                                         break;
288                                 }
289                         if(found)
290                                 return y*width+x;
291                         vertices += numcomponents*2;
292                 }
293                 vertices += numcomponents*(width-1);
294         }
295         return -1;
296 }
298 #define SIDE_INVALID -1
299 #define SIDE_X 0
300 #define SIDE_Y 1
302 static int GetSide(int p1, int p2, int width, int height, int *pointdist)
303 {
304         int x1 = p1 % width, y1 = p1 / width;
305         int x2 = p2 % width, y2 = p2 / width;
306         if (p1 < 0 || p2 < 0)
307                 return SIDE_INVALID;
308         if (x1 == x2)
309         {
310                 if (y1 != y2)
311                 {
312                         *pointdist = abs(y2 - y1);
313                         return SIDE_Y;
314                 }
315                 else
316                         return SIDE_INVALID;
317         }
318         else if (y1 == y2)
319         {
320                 *pointdist = abs(x2 - x1);
321                 return SIDE_X;
322         }
323         else
324                 return SIDE_INVALID;
325 }
327 // Increase tesselation of one of two touching patches to make a seamless connection between them
328 // Returns 0 in case if patches were not modified, otherwise 1
329 int Q3PatchAdjustTesselation(int numcomponents, patchinfo_t *patch1, float *patchvertices1, patchinfo_t *patch2, float *patchvertices2)
330 {
331         // what we are doing here is:
332         //   we take for each corner of one patch
333         //   and check if the other patch contains that corner
334         //   once we have a pair of such matches
336         struct {int id1,id2;} commonverts;
337         int i, j, k, side1, side2, *tess1, *tess2;
338         int dist1 = 0, dist2 = 0;
339         qbool modified = false;
341         // Potential paired vertices (corners of the first patch)
342         commonverts.id1 = 0;
343         commonverts.id1 = patch1->xsize-1;
344         commonverts.id1 = patch1->xsize*(patch1->ysize-1);
345         commonverts.id1 = patch1->xsize*patch1->ysize-1;
346         for (i=0;i<4;++i)
347                 commonverts[i].id2 = FindEqualOddVertexInArray(numcomponents, patchvertices1+numcomponents*commonverts[i].id1, patchvertices2, patch2->xsize, patch2->ysize);
349         // Corners of the second patch
350         commonverts.id2 = 0;
351         commonverts.id2 = patch2->xsize-1;
352         commonverts.id2 = patch2->xsize*(patch2->ysize-1);
353         commonverts.id2 = patch2->xsize*patch2->ysize-1;
354         for (i=4;i<8;++i)
355                 commonverts[i].id1 = FindEqualOddVertexInArray(numcomponents, patchvertices2+numcomponents*commonverts[i].id2, patchvertices1, patch1->xsize, patch1->ysize);
357         for (i=0;i<8;++i)
358                 for (j=i+1;j<8;++j)
359                 {
360                         side1 = GetSide(commonverts[i].id1,commonverts[j].id1,patch1->xsize,patch1->ysize,&dist1);
361                         side2 = GetSide(commonverts[i].id2,commonverts[j].id2,patch2->xsize,patch2->ysize,&dist2);
363                         if (side1 == SIDE_INVALID || side2 == SIDE_INVALID)
364                                 continue;
366                         if(dist1 != dist2)
367                         {
368                                 // no patch welding if the resolutions mismatch
369                                 continue;
370                         }
372                         // Update every lod level
373                         for (k=0;k<PATCH_LODS_NUM;++k)
374                         {
375                                 tess1 = side1 == SIDE_X ? &patch1->lods[k].xtess : &patch1->lods[k].ytess;
376                                 tess2 = side2 == SIDE_X ? &patch2->lods[k].xtess : &patch2->lods[k].ytess;
377                                 if (*tess1 != *tess2)
378                                 {
379                                         if (*tess1 < *tess2)
380                                                 *tess1 = *tess2;
381                                         else
382                                                 *tess2 = *tess1;
383                                         modified = true;
384                                 }
385                         }
386                 }
388         return modified;
389 }
391 #undef SIDE_INVALID
392 #undef SIDE_X
393 #undef SIDE_Y
395 // calculates elements for a grid of vertices
396 // (such as those produced by Q3PatchTesselate)
397 // (note: width and height are the actual vertex size, this produces
398 // (width-1)*(height-1)*2 triangles, 3 elements each)
399 void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
400 {
401         int x, y, row0, row1;
402         for (y = 0;y < height - 1;y++)
403         {
404                 if(y % 2)
405                 {
406                         // swap the triangle order in odd rows as optimization for collision stride
407                         row0 = firstvertex + (y + 0) * width + width - 2;
408                         row1 = firstvertex + (y + 1) * width + width - 2;
409                         for (x = 0;x < width - 1;x++)
410                         {
411                                 *elements++ = row1;
412                                 *elements++ = row1 + 1;
413                                 *elements++ = row0 + 1;
414                                 *elements++ = row0;
415                                 *elements++ = row1;
416                                 *elements++ = row0 + 1;
417                                 row0--;
418                                 row1--;
419                         }
420                 }
421                 else
422                 {
423                         row0 = firstvertex + (y + 0) * width;
424                         row1 = firstvertex + (y + 1) * width;
425                         for (x = 0;x < width - 1;x++)
426                         {
427                                 *elements++ = row0;
428                                 *elements++ = row1;
429                                 *elements++ = row0 + 1;
430                                 *elements++ = row1;
431                                 *elements++ = row1 + 1;
432                                 *elements++ = row0 + 1;
433                                 row0++;
434                                 row1++;
435                         }
436                 }
437         }
438 }