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