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 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
13 void QuadraticBSplineSubdivideFloat(int inpoints, int components, const float *in, int instride, float *out, int outstride)
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
23 if (components == 1 && instride == (int)sizeof(float) && outstride == instride)
25 // simple case, single component and no special stride
26 float cpprev0 = 0, cpcurr0 = 0, cpnext0;
28 for (s = 0;s < inpoints - 1;s++)
36 // 50% flattened control point
37 // cp1 = average(cp1, average(cp0, cp2));
38 *out++ = (cpcurr0 + (cpprev0 + cpnext0) * 0.5f) * 0.5f;
42 // copy the control point directly
46 // mid = average(cp0, cp1);
47 *out++ = (cpcurr0 + cpnext0) * 0.5f;
49 // copy the final control point
54 // multiple components or stride is used (complex case)
56 float cpprev[4], cpcurr[4], cpnext[4];
57 // check if there are too many components for the buffers
60 // more components can be handled, but slowly, by calling self multiple times...
61 for (c = 0;c < components;c++, in++, out++)
62 QuadraticBSplineSubdivideFloat(inpoints, 1, in, instride, out, outstride);
65 for (c = 0;c < components;c++)
67 (unsigned char *)in += instride;
68 for (s = 0;s < inpoints - 1;s++)
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++)
76 (unsigned char *)in += instride;
77 // the end points are copied as-is
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;
87 // copy the control point directly
88 for (c = 0;c < components;c++)
91 (unsigned char *)out += outstride;
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;
98 // copy the final control point
99 for (c = 0;c < components;c++)
101 //(unsigned char *)out += outstride;
105 // note: out must have enough room!
106 // (see finalwidth/finalheight calcs below)
107 void QuadraticBSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
109 int finalwidth, finalheight, xstep, ystep, x, y, c;
112 // error out on various bogus conditions
113 if (xlevel < 0 || ylevel < 0 || xlevel > 16 || ylevel > 16 || cpwidth < 3 || cpheight < 3)
118 finalwidth = (cpwidth - 1) * xstep + 1;
119 finalheight = (cpheight - 1) * ystep + 1;
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++)
126 if (xlevel == 1 && ylevel == 0)
128 for (y = 0;y < finalheight;y++)
129 QuadraticBSplineSubdivideFloat(cpwidth, components, in + y * cpwidth * components, sizeof(float) * components, out + y * finalwidth * components, sizeof(float) * components);
132 if (xlevel == 0 && ylevel == 1)
134 for (x = 0;x < finalwidth;x++)
135 QuadraticBSplineSubdivideFloat(cpheight, components, in + x * components, sizeof(float) * cpwidth * components, out + x * components, sizeof(float) * finalwidth * components);
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++)
145 // subdivide in place in the destination buffer
146 while (xstep > 1 || ystep > 1)
151 for (y = 0;y < finalheight;y += ystep)
152 QuadraticBSplineSubdivideFloat(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;
158 for (x = 0;x < finalwidth;x += xstep)
159 QuadraticBSplineSubdivideFloat(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;
165 void QuadraticBSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
167 int c, x, y, outwidth, outheight, halfstep, xstep, ystep;
168 float prev, curr, next;
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)
182 halfstep = xstep >> 1;
183 for (y = 0;y < outheight;y += ystep)
185 for (c = 0;c < components;c++)
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)
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;
209 halfstep = ystep >> 1;
210 for (x = 0;x < outwidth;x += xstep)
212 for (c = 0;c < components;c++)
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)
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;
234 // flatten control points on X
235 for (y = 0;y < outheight;y += ystep)
237 for (c = 0;c < components;c++)
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)
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;
252 // flatten control points on Y
253 for (x = 0;x < outwidth;x += xstep)
255 for (c = 0;c < components;c++)
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)
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;
272 for (y = ystep;y < outheight - ystep;y += ystep)
274 for (c = 0;c < components;c++)
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])
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;
291 void QuadraticBSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
293 int outwidth, outheight;
294 outwidth = ((cpwidth - 1) << xlevel) + 1;
295 outheight = ((cpheight - 1) << ylevel) + 1;
296 for (y = 0;y < cpheight;y++)
298 for (x = 0;x < cpwidth;x++)
300 for (c = 0;c < components;c++)
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++)
306 for (sx = 0;sx < expandx;sx++)
308 d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
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
324 float QuadraticBSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
328 float deviation, squareddeviation, bestsquareddeviation;
329 bestsquareddeviation = 0;
330 for (y = 0;y < cpheight;y++)
332 for (x = 1;x < cpwidth-1;x++)
334 squareddeviation = 0;
335 for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
337 deviation = cp[0] * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
338 squareddeviation += deviation*deviation;
340 if (bestsquareddeviation < squareddeviation)
341 bestsquareddeviation = squareddeviation;
344 return (float)sqrt(bestsquareddeviation);
347 float QuadraticBSplinePatchLargestDeviationOnY(int cpwidth, int cpheight, int components, const float *in)
351 float deviation, squareddeviation, bestsquareddeviation;
352 bestsquareddeviation = 0;
353 for (y = 1;y < cpheight-1;y++)
355 for (x = 0;x < cpwidth;x++)
357 squareddeviation = 0;
358 for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
360 deviation = cp[0] * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
361 squareddeviation += deviation*deviation;
363 if (bestsquareddeviation < squareddeviation)
364 bestsquareddeviation = squareddeviation;
367 return (float)sqrt(bestsquareddeviation);
370 int QuadraticBSplinePatchSubdivisionLevelForDeviation(float deviation, float level1tolerance, int levellimit)
373 // count the automatic flatten step which reduces deviation by 50%
375 // count the levels to subdivide to come under the tolerance
376 for (level = 0;level < levellimit && deviation > level1tolerance;level++)
381 int QuadraticBSplinePatchSubdivisionLevelOnX(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
383 return QuadraticBSplinePatchSubdivisionLevelForDeviation(QuadraticBSplinePatchLargestDeviationOnX(cpwidth, cpheight, components, in), level1tolerance, levellimit);
386 int QuadraticBSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
388 return QuadraticBSplinePatchSubdivisionLevelForDeviation(QuadraticBSplinePatchLargestDeviationOnY(cpwidth, cpheight, components, in), level1tolerance, levellimit);
392 // 1: flat (0th dimension)
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
403 // n: arbitrary dimension bspline
404 double factorial(int n)
409 for (i = 1;i < n;i++)
413 double bsplinesample(int dimensions, double t, double *param)
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);