Experimental new patch welding code by "someone", fixing the problem with seams when...
[divverent/darkplaces.git] / curves.c
1
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
5
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.
10
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
22
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 */
41
42 #include "quakedef.h"
43 #include "mathlib.h"
44
45 #include <math.h>
46 #include "curves.h"
47
48 // Calculate number of resulting vertex rows/columns by given patch size and tesselation factor
49 // tess=0 means that we reduce detalization of base 3x3 patches by removing middle row and column of vertices
50 // "DimForTess" is "DIMension FOR TESSelation factor"
51 // 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)
52 int Q3PatchDimForTess(int size, int tess)
53 {
54         if (tess > 0)
55                 return (size - 1) * tess + 1;
56         else if (tess == 0)
57                 return (size - 1) / 2 + 1;
58         else
59                 return 0; // Maybe warn about wrong tess here?
60 }
61
62 // usage:
63 // to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
64 // Q3PatchSubdivideFloat(3, sizeof(float[3]), outvertices, 5, 5, sizeof(float[3]), patchvertices, 4, 4);
65 void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
66 {
67         int k, l, x, y, component, outputwidth = Q3PatchDimForTess(patchwidth, tesselationwidth);
68         float px, py, *v, a, b, c, *cp[3][3], temp[3][64];
69         int xmax = max(1, 2*tesselationwidth);
70         int ymax = max(1, 2*tesselationheight);
71         
72         // iterate over the individual 3x3 quadratic spline surfaces one at a time
73         // expanding them to fill the output array (with some overlap to ensure
74         // the edges are filled)
75         for (k = 0;k < patchheight-1;k += 2)
76         {
77                 for (l = 0;l < patchwidth-1;l += 2)
78                 {
79                         // set up control point pointers for quicker lookup later
80                         for (y = 0;y < 3;y++)
81                                 for (x = 0;x < 3;x++)
82                                         cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
83                         // for each row...
84                         for (y = 0;y <= ymax;y++)
85                         {
86                                 // calculate control points for this row by collapsing the 3
87                                 // rows of control points to one row using py
88                                 py = (float)y / (float)ymax;
89                                 // calculate quadratic spline weights for py
90                                 a = ((1.0f - py) * (1.0f - py));
91                                 b = ((1.0f - py) * (2.0f * py));
92                                 c = ((       py) * (       py));
93                                 for (component = 0;component < numcomponents;component++)
94                                 {
95                                         temp[0][component] = cp[0][0][component] * a + cp[1][0][component] * b + cp[2][0][component] * c;
96                                         temp[1][component] = cp[0][1][component] * a + cp[1][1][component] * b + cp[2][1][component] * c;
97                                         temp[2][component] = cp[0][2][component] * a + cp[1][2][component] * b + cp[2][2][component] * c;
98                                 }
99                                 // fetch a pointer to the beginning of the output vertex row
100                                 v = (float *)((unsigned char *)outputvertices + ((k * ymax / 2 + y) * outputwidth + l * xmax / 2) * outputstride);
101                                 // for each column of the row...
102                                 for (x = 0;x <= xmax;x++)
103                                 {
104                                         // calculate point based on the row control points
105                                         px = (float)x / (float)xmax;
106                                         // calculate quadratic spline weights for px
107                                         // (could be precalculated)
108                                         a = ((1.0f - px) * (1.0f - px));
109                                         b = ((1.0f - px) * (2.0f * px));
110                                         c = ((       px) * (       px));
111                                         for (component = 0;component < numcomponents;component++)
112                                                 v[component] = temp[0][component] * a + temp[1][component] * b + temp[2][component] * c;
113                                         // advance to next output vertex using outputstride
114                                         // (the next vertex may not be directly following this
115                                         // one, as this may be part of a larger structure)
116                                         v = (float *)((unsigned char *)v + outputstride);
117                                 }
118                         }
119                 }
120         }
121 #if 0
122         // enable this if you want results printed out
123         printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
124         for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
125         {
126                 for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
127                 {
128                         printf("(");
129                         for (component = 0;component < numcomponents;component++)
130                                 printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
131                         printf(") ");
132                 }
133                 printf("\n");
134         }
135         printf("}\n");
136 #endif
137 }
138
139 static int Q3PatchTesselation(float bestsquareddeviation, float tolerance)
140 {
141         float f;
142         f = sqrt(bestsquareddeviation) / tolerance;
143         //if(f < 0.25) // REALLY 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)) + 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 }
157
158 // returns how much tesselation of each segment is needed to remain under tolerance
159 int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
160 {
161         int c, x, y;
162         const float *patch;
163         float deviation, squareddeviation, bestsquareddeviation;
164         bestsquareddeviation = 0;
165         for (y = 0;y < patchheight;y++)
166         {
167                 for (x = 0;x < patchwidth-1;x += 2)
168                 {
169                         squareddeviation = 0;
170                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
171                         {
172                                 deviation = patch[components] * 0.5f - patch[0] * 0.25f - patch[2*components] * 0.25f;
173                                 squareddeviation += deviation*deviation;
174                         }
175                         if (bestsquareddeviation < squareddeviation)
176                                 bestsquareddeviation = squareddeviation;
177                 }
178         }
179         return Q3PatchTesselation(bestsquareddeviation, tolerance);
180 }
181
182 // returns how much tesselation of each segment is needed to remain under tolerance
183 int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
184 {
185         int c, x, y;
186         const float *patch;
187         float deviation, squareddeviation, bestsquareddeviation;
188         bestsquareddeviation = 0;
189         for (y = 0;y < patchheight-1;y += 2)
190         {
191                 for (x = 0;x < patchwidth;x++)
192                 {
193                         squareddeviation = 0;
194                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
195                         {
196                                 deviation = patch[patchwidth*components] * 0.5f - patch[0] * 0.25f - patch[2*patchwidth*components] * 0.25f;
197                                 squareddeviation += deviation*deviation;
198                         }
199                         if (bestsquareddeviation < squareddeviation)
200                                 bestsquareddeviation = squareddeviation;
201                 }
202         }
203         return Q3PatchTesselation(bestsquareddeviation, tolerance);
204 }
205
206 // Find an equal vertex in array. Check only vertices with odd X and Y
207 static int FindEqualOddVertexInArray(int numcomponents, float *vertex, float *vertices, int width, int height)
208 {
209         int x, y, j;
210         for (y=0; y<height; y+=2)
211         {
212                 for (x=0; x<width; x+=2)
213                 {
214                         qboolean found = true;
215                         for (j=0; j<numcomponents; j++)
216                                 if (fabs(*(vertex+j) - *(vertices+j)) > 0.05)
217                                 // div0: this is notably smaller than the smallest radiant grid
218                                 // but large enough so we don't need to get scared of roundoff
219                                 // errors
220                                 {
221                                         found = false;
222                                         break;
223                                 }
224                         if(found)
225                                 return y*width+x;
226                         vertices += numcomponents*2;
227                 }
228                 vertices += numcomponents*(width-1);
229         }
230         return -1;
231 }
232
233 #define SIDE_INVALID -1
234 #define SIDE_X 0
235 #define SIDE_Y 1
236
237 static int GetSide(int p1, int p2, int width, int height, int *pointdist)
238 {
239         int x1 = p1 % width, y1 = p1 / width;
240         int x2 = p2 % width, y2 = p2 / width;
241         if (p1 < 0 || p2 < 0)
242                 return SIDE_INVALID;
243         if (x1 == x2)
244         {
245                 if (y1 != y2)
246                 {
247                         *pointdist = abs(y2 - y1);
248                         return SIDE_Y;
249                 }
250                 else
251                         return SIDE_INVALID;
252         }
253         else if (y1 == y2)
254         {
255                 *pointdist = abs(x2 - x1);
256                 return SIDE_X;
257         }
258         else
259                 return SIDE_INVALID;
260 }
261
262 // Increase tesselation of one of two touching patches to make a seamless connection between them
263 // Returns 0 in case if patches were not modified, otherwise 1
264 int Q3PatchAdjustTesselation(int numcomponents, patchinfo_t *patch1, float *patchvertices1, patchinfo_t *patch2, float *patchvertices2)
265 {
266         // what we are doing here is:
267         //   we take for each corner of one patch
268         //   and check if the other patch contains that corner
269         //   once we have a pair of such matches
270
271         struct {int id1,id2;} commonverts[8];
272         int i, j, k, side1, side2, *tess1, *tess2;
273         int dist1, dist2;
274         qboolean modified = false;
275
276         // Potential paired vertices (corners of the first patch)
277         commonverts[0].id1 = 0;
278         commonverts[1].id1 = patch1->xsize-1;
279         commonverts[2].id1 = patch1->xsize*(patch1->ysize-1);
280         commonverts[3].id1 = patch1->xsize*patch1->ysize-1;
281         for (i=0;i<4;++i)
282                 commonverts[i].id2 = FindEqualOddVertexInArray(numcomponents, patchvertices1+numcomponents*commonverts[i].id1, patchvertices2, patch2->xsize, patch2->ysize);
283
284         // Corners of the second patch
285         commonverts[4].id2 = 0;
286         commonverts[5].id2 = patch2->xsize-1;
287         commonverts[6].id2 = patch2->xsize*(patch2->ysize-1);
288         commonverts[7].id2 = patch2->xsize*patch2->ysize-1;
289         for (i=4;i<8;++i)
290                 commonverts[i].id1 = FindEqualOddVertexInArray(numcomponents, patchvertices2+numcomponents*commonverts[i].id2, patchvertices1, patch1->xsize, patch1->ysize);
291
292         for (i=0;i<8;++i)
293                 for (j=i+1;j<8;++j)
294                 {
295                         side1 = GetSide(commonverts[i].id1,commonverts[j].id1,patch1->xsize,patch1->ysize,&dist1);
296                         side2 = GetSide(commonverts[i].id2,commonverts[j].id2,patch2->xsize,patch2->ysize,&dist2);
297
298                         if (side1 == SIDE_INVALID || side2 == SIDE_INVALID)
299                                 continue;
300
301                         if(dist1 != dist2)
302                         {
303                                 // no patch welding if the resolutions mismatch
304                                 continue;
305                         }
306
307                         // Update every lod level
308                         for (k=0;k<PATCH_LODS_NUM;++k)
309                         {
310                                 tess1 = side1 == SIDE_X ? &patch1->lods[k].xtess : &patch1->lods[k].ytess;
311                                 tess2 = side2 == SIDE_X ? &patch2->lods[k].xtess : &patch2->lods[k].ytess;
312                                 if (*tess1 != *tess2)
313                                 {
314                                         if (*tess1 < *tess2)
315                                                 *tess1 = *tess2;
316                                         else
317                                                 *tess2 = *tess1;
318                                         modified = true;
319                                 }
320                         }
321                 }
322
323         return modified;
324 }
325
326 #undef SIDE_INVALID 
327 #undef SIDE_X
328 #undef SIDE_Y
329
330 // calculates elements for a grid of vertices
331 // (such as those produced by Q3PatchTesselate)
332 // (note: width and height are the actual vertex size, this produces
333 // (width-1)*(height-1)*2 triangles, 3 elements each)
334 void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
335 {
336         int x, y, row0, row1;
337         for (y = 0;y < height - 1;y++)
338         {
339                 row0 = firstvertex + (y + 0) * width;
340                 row1 = firstvertex + (y + 1) * width;
341                 for (x = 0;x < width - 1;x++)
342                 {
343                         *elements++ = row0;
344                         *elements++ = row1;
345                         *elements++ = row0 + 1;
346                         *elements++ = row1;
347                         *elements++ = row1 + 1;
348                         *elements++ = row0 + 1;
349                         row0++;
350                         row1++;
351                 }
352         }
353 }
354