more documentation on B-Splines
[xonotic/darkplaces.git] / curves.c
1
2 // this code written by Forest Hale, on 2003-08-23, and placed into public domain
3 // this code deals with quadratic splines (minimum of 3 points), the same kind used in Quake3 maps.
4
5 // LordHavoc's rant on misuse of the name 'bezier': many people seem to think that bezier is a generic term for splines, but it is not, it is a term for a specific type of bspline (4 control points, cubic bspline), bsplines are the generalization of the bezier spline to support dimensions other than just cubic.
6 // this implements Quadratic BSpline surfaces
7
8 #include <math.h>
9 #include "curves.h"
10 #include "zone.h"
11
12 #if 0
13 void QuadraticSplineSubdivideFloat(int inpoints, int components, const float *in, int instride, float *out, int outstride)
14 {
15         int s;
16         // the input (control points) is read as a stream of points, and buffered
17         // by the cpprev, cpcurr, and cpnext variables (to allow subdivision in
18         // overlapping memory buffers, even subdivision in-place with pre-spaced
19         // control points in the buffer)
20         // the output (resulting curve) is written as a stream of points
21         // this subdivision is meant to be repeated until the desired flatness
22         // level is reached
23         if (components == 1 && instride == (int)sizeof(float) && outstride == instride)
24         {
25                 // simple case, single component and no special stride
26                 float cpprev0 = 0, cpcurr0 = 0, cpnext0;
27                 cpnext0 = *in++;
28                 for (s = 0;s < inpoints - 1;s++)
29                 {
30                         cpprev0 = cpcurr0;
31                         cpcurr0 = cpnext0;
32                         if (s < inpoints - 1)
33                                 cpnext0 = *in++;
34                         if (s > 0)
35                         {
36                                 // 50% flattened control point
37                                 // cp1 = average(cp1, average(cp0, cp2));
38                                 *out++ = (cpcurr0 + (cpprev0 + cpnext0) * 0.5f) * 0.5f;
39                         }
40                         else
41                         {
42                                 // copy the control point directly
43                                 *out++ = cpcurr0;
44                         }
45                         // midpoint
46                         // mid = average(cp0, cp1);
47                         *out++ = (cpcurr0 + cpnext0) * 0.5f;
48                 }
49                 // copy the final control point
50                 *out++ = cpnext0;
51         }
52         else
53         {
54                 // multiple components or stride is used (complex case)
55                 int c;
56                 float cpprev[4], cpcurr[4], cpnext[4];
57                 // check if there are too many components for the buffers
58                 if (components > 1)
59                 {
60                         // more components can be handled, but slowly, by calling self multiple times...
61                         for (c = 0;c < components;c++, in++, out++)
62                                 QuadraticSplineSubdivideFloat(inpoints, 1, in, instride, out, outstride);
63                         return;
64                 }
65                 for (c = 0;c < components;c++)
66                         cpnext[c] = in[c];
67                 (unsigned char *)in += instride;
68                 for (s = 0;s < inpoints - 1;s++)
69                 {
70                         for (c = 0;c < components;c++)
71                                 cpprev[c] = cpcurr[c];
72                         for (c = 0;c < components;c++)
73                                 cpcurr[c] = cpnext[c];
74                         for (c = 0;c < components;c++)
75                                 cpnext[c] = in[c];
76                         (unsigned char *)in += instride;
77                         // the end points are copied as-is
78                         if (s > 0)
79                         {
80                                 // 50% flattened control point
81                                 // cp1 = average(cp1, average(cp0, cp2));
82                                 for (c = 0;c < components;c++)
83                                         out[c] = (cpcurr[c] + (cpprev[c] + cpnext[c]) * 0.5f) * 0.5f;
84                         }
85                         else
86                         {
87                                 // copy the control point directly
88                                 for (c = 0;c < components;c++)
89                                         out[c] = cpcurr[c];
90                         }
91                         (unsigned char *)out += outstride;
92                         // midpoint
93                         // mid = average(cp0, cp1);
94                         for (c = 0;c < components;c++)
95                                 out[c] = (cpcurr[c] + cpnext[c]) * 0.5f;
96                         (unsigned char *)out += outstride;
97                 }
98                 // copy the final control point
99                 for (c = 0;c < components;c++)
100                         out[c] = cpnext[c];
101                 //(unsigned char *)out += outstride;
102         }
103 }
104
105 // note: out must have enough room!
106 // (see finalwidth/finalheight calcs below)
107 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
108 {
109         int finalwidth, finalheight, xstep, ystep, x, y, c;
110         float *o;
111
112         // error out on various bogus conditions
113         if (xlevel < 0 || ylevel < 0 || xlevel > 16 || ylevel > 16 || cpwidth < 3 || cpheight < 3)
114                 return;
115
116         xstep = 1 << xlevel;
117         ystep = 1 << ylevel;
118         finalwidth = (cpwidth - 1) * xstep + 1;
119         finalheight = (cpheight - 1) * ystep + 1;
120
121         for (y = 0;y < finalheight;y++)
122                 for (x = 0;x < finalwidth;x++)
123                         for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
124                                 o[c] = 0;
125
126         if (xlevel == 1 && ylevel == 0)
127         {
128                 for (y = 0;y < finalheight;y++)
129                         QuadraticSplineSubdivideFloat(cpwidth, components, in + y * cpwidth * components, sizeof(float) * components, out + y * finalwidth * components, sizeof(float) * components);
130                 return;
131         }
132         if (xlevel == 0 && ylevel == 1)
133         {
134                 for (x = 0;x < finalwidth;x++)
135                         QuadraticSplineSubdivideFloat(cpheight, components, in + x * components, sizeof(float) * cpwidth * components, out + x * components, sizeof(float) * finalwidth * components);
136                 return;
137         }
138
139         // copy control points into correct positions in destination buffer
140         for (y = 0;y < finalheight;y += ystep)
141                 for (x = 0;x < finalwidth;x += xstep)
142                         for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
143                                 o[c] = *in++;
144
145         // subdivide in place in the destination buffer
146         while (xstep > 1 || ystep > 1)
147         {
148                 if (xstep > 1)
149                 {
150                         xstep >>= 1;
151                         for (y = 0;y < finalheight;y += ystep)
152                                 QuadraticSplineSubdivideFloat(cpwidth, components, out + y * finalwidth * components, sizeof(float) * xstep * 2 * components, out + y * finalwidth * components, sizeof(float) * xstep * components);
153                         cpwidth = (cpwidth - 1) * 2 + 1;
154                 }
155                 if (ystep > 1)
156                 {
157                         ystep >>= 1;
158                         for (x = 0;x < finalwidth;x += xstep)
159                                 QuadraticSplineSubdivideFloat(cpheight, components, out + x * components, sizeof(float) * ystep * 2 * finalwidth * components, out + x * components, sizeof(float) * ystep * finalwidth * components);
160                         cpheight = (cpheight - 1) * 2 + 1;
161                 }
162         }
163 }
164 #elif 1
165 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
166 {
167         int c, x, y, outwidth, outheight, halfstep, xstep, ystep;
168         float prev, curr, next;
169         xstep = 1 << xlevel;
170         ystep = 1 << ylevel;
171         outwidth = ((cpwidth - 1) * xstep) + 1;
172         outheight = ((cpheight - 1) * ystep) + 1;
173         for (y = 0;y < cpheight;y++)
174                 for (x = 0;x < cpwidth;x++)
175                         for (c = 0;c < components;c++)
176                                 out[(y * ystep * outwidth + x * xstep) * components + c] = in[(y * cpwidth + x) * components + c];
177         while (xstep > 1 || ystep > 1)
178         {
179                 if (xstep >= ystep)
180                 {
181                         // subdivide on X
182                         halfstep = xstep >> 1;
183                         for (y = 0;y < outheight;y += ystep)
184                         {
185                                 for (c = 0;c < components;c++)
186                                 {
187                                         x = xstep;
188                                         // fetch first two control points 
189                                         prev = out[(y * outwidth + (x - xstep)) * components + c];
190                                         curr = out[(y * outwidth + x) * components + c];
191                                         // create first midpoint
192                                         out[(y * outwidth + (x - halfstep)) * components + c] = (curr + prev) * 0.5f;
193                                         for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
194                                         {
195                                                 // fetch next control point
196                                                 next = out[(y * outwidth + (x + xstep)) * components + c];
197                                                 // flatten central control point 
198                                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
199                                                 // create following midpoint
200                                                 out[(y * outwidth + (x + halfstep)) * components + c] = (curr + next) * 0.5f;
201                                         }
202                                 }
203                         }
204                         xstep >>= 1;
205                 }
206                 else
207                 {
208                         // subdivide on Y
209                         halfstep = ystep >> 1;
210                         for (x = 0;x < outwidth;x += xstep)
211                         {
212                                 for (c = 0;c < components;c++)
213                                 {
214                                         y = ystep;
215                                         // fetch first two control points 
216                                         prev = out[((y - ystep) * outwidth + x) * components + c];
217                                         curr = out[(y * outwidth + x) * components + c];
218                                         // create first midpoint
219                                         out[((y - halfstep) * outwidth + x) * components + c] = (curr + prev) * 0.5f;
220                                         for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
221                                         {
222                                                 // fetch next control point
223                                                 next = out[((y + ystep) * outwidth + x) * components + c];
224                                                 // flatten central control point 
225                                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
226                                                 // create following midpoint
227                                                 out[((y + halfstep) * outwidth + x) * components + c] = (curr + next) * 0.5f;
228                                         }
229                                 }
230                         }
231                         ystep >>= 1;
232                 }
233         }
234         // flatten control points on X
235         for (y = 0;y < outheight;y += ystep)
236         {
237                 for (c = 0;c < components;c++)
238                 {
239                         x = xstep;
240                         // fetch first two control points 
241                         prev = out[(y * outwidth + (x - xstep)) * components + c];
242                         curr = out[(y * outwidth + x) * components + c];
243                         for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
244                         {
245                                 // fetch next control point 
246                                 next = out[(y * outwidth + (x + xstep)) * components + c];
247                                 // flatten central control point 
248                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
249                         }
250                 }
251         }
252         // flatten control points on Y
253         for (x = 0;x < outwidth;x += xstep)
254         {
255                 for (c = 0;c < components;c++)
256                 {
257                         y = ystep;
258                         // fetch first two control points 
259                         prev = out[((y - ystep) * outwidth + x) * components + c];
260                         curr = out[(y * outwidth + x) * components + c];
261                         for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
262                         {
263                                 // fetch next control point 
264                                 next = out[((y + ystep) * outwidth + x) * components + c];
265                                 // flatten central control point 
266                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
267                         }
268                 }
269         }
270
271         /*
272         for (y = ystep;y < outheight - ystep;y += ystep)
273         {
274                 for (c = 0;c < components;c++)
275                 {
276                         for (x = xstep, outp = out + (y * outwidth + x) * components + c, prev = outp[-xstep * components], curr = outp[0], next = outp[xstep * components];x < outwidth;x += xstep, outp += ystep * outwidth * components, prev = curr, curr = next, next = outp[xstep * components])
277                         {
278                                 // midpoint
279                                 outp[-halfstep * components] = (prev + curr) * 0.5f;
280                                 // flatten control point
281                                 outp[0] = (curr + (prev + next) * 0.5f) * 0.5f;
282                                 // next midpoint (only needed for end segment)
283                                 outp[halfstep * components] = (next + curr) * 0.5f;
284                         }
285                 }
286         }
287         */
288 }
289 #else
290 // unfinished code
291 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
292 {
293         int outwidth, outheight;
294         outwidth = ((cpwidth - 1) << xlevel) + 1;
295         outheight = ((cpheight - 1) << ylevel) + 1;
296         for (y = 0;y < cpheight;y++)
297         {
298                 for (x = 0;x < cpwidth;x++)
299                 {
300                         for (c = 0;c < components;c++)
301                         {
302                                 inp = in + (y * cpwidth + x) * components + c;
303                                 outp = out + ((y<<ylevel) * outwidth + (x<<xlevel)) * components + c;
304                                 for (sy = 0;sy < expandy;sy++)
305                                 {
306                                         for (sx = 0;sx < expandx;sx++)
307                                         {
308                                                 d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
309                                         }
310                                 }
311                         }
312                 }
313         }
314 }
315 #endif
316
317 /*
318 0.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 1.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 0.00000 deviation: 0.5
319 0.00000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.00000 deviation: 0.125
320 0.00000 ?.????? 0.25000 ?.????? 0.37500 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.37500 ?.????? 0.25000 ?.????? 0.00000 deviation: 0.03125
321 0.00000 0.12500 0.21875 0.31250 0.37500 0.43750 0.46875 0.50000 0.50000 0.50000 0.46875 0.43750 0.37500 0.31250 0.21875 0.12500 0.00000 deviation: not available
322 */
323
324 float QuadraticSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
325 {
326         int c, x, y;
327         const float *cp;
328         float deviation, squareddeviation, bestsquareddeviation;
329         bestsquareddeviation = 0;
330         for (y = 0;y < cpheight;y++)
331         {
332                 for (x = 1;x < cpwidth-1;x++)
333                 {
334                         squareddeviation = 0;
335                         for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
336                         {
337                                 deviation = cp[0] * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
338                                 squareddeviation += deviation*deviation;
339                         }
340                         if (bestsquareddeviation < squareddeviation)
341                                 bestsquareddeviation = squareddeviation;
342                 }
343         }
344         return (float)sqrt(bestsquareddeviation);
345 }
346
347 float QuadraticSplinePatchLargestDeviationOnY(int cpwidth, int cpheight, int components, const float *in)
348 {
349         int c, x, y;
350         const float *cp;
351         float deviation, squareddeviation, bestsquareddeviation;
352         bestsquareddeviation = 0;
353         for (y = 1;y < cpheight-1;y++)
354         {
355                 for (x = 0;x < cpwidth;x++)
356                 {
357                         squareddeviation = 0;
358                         for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
359                         {
360                                 deviation = cp[0] * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
361                                 squareddeviation += deviation*deviation;
362                         }
363                         if (bestsquareddeviation < squareddeviation)
364                                 bestsquareddeviation = squareddeviation;
365                 }
366         }
367         return (float)sqrt(bestsquareddeviation);
368 }
369
370 int QuadraticSplinePatchSubdivisionLevelForDeviation(float deviation, float level1tolerance, int levellimit)
371 {
372         int level;
373         // count the automatic flatten step which reduces deviation by 50%
374         deviation *= 0.5f;
375         // count the levels to subdivide to come under the tolerance
376         for (level = 0;level < levellimit && deviation > level1tolerance;level++)
377                 deviation *= 0.25f;
378         return level;
379 }
380
381 int QuadraticSplinePatchSubdivisionLevelOnX(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
382 {
383         return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnX(cpwidth, cpheight, components, in), level1tolerance, levellimit);
384 }
385
386 int QuadraticSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
387 {
388         return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnY(cpwidth, cpheight, components, in), level1tolerance, levellimit);
389 }
390
391 /*
392         // 1: flat (0th dimension)
393         o = a
394         // 2: linear (1st dimension)
395         o = a * (1 - t) + b * t
396         // 3: quadratic bspline (2nd dimension)
397         o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
398         // 4: cubic (bezier) bspline (3rd dimension)
399         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
400         // 5: quartic bspline (4th dimension)
401         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
402
403         // n: arbitrary dimension bspline
404 double factorial(int n)
405 {
406         int i;
407         double f;
408         f = 1;
409         for (i = 1;i < n;i++)
410                 f = f * i;
411         return f;
412 }
413 double bsplinesample(int dimensions, double t, double *param)
414 {
415         double o = 0;
416         for (i = 0;i < dimensions + 1;i++)
417                 o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
418 }
419 */
420