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).
12 void QuadraticSplineSubdivideFloat(int inpoints, int components, const float *in, int instride, float *out, int outstride)
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
22 if (components == 1 && instride == (int)sizeof(float) && outstride == instride)
24 // simple case, single component and no special stride
25 float cpprev0 = 0, cpcurr0 = 0, cpnext0;
27 for (s = 0;s < inpoints - 1;s++)
35 // 50% flattened control point
36 // cp1 = average(cp1, average(cp0, cp2));
37 *out++ = (cpcurr0 + (cpprev0 + cpnext0) * 0.5f) * 0.5f;
41 // copy the control point directly
45 // mid = average(cp0, cp1);
46 *out++ = (cpcurr0 + cpnext0) * 0.5f;
48 // copy the final control point
53 // multiple components or stride is used (complex case)
55 float cpprev[4], cpcurr[4], cpnext[4];
56 // check if there are too many components for the buffers
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);
64 for (c = 0;c < components;c++)
66 (unsigned char *)in += instride;
67 for (s = 0;s < inpoints - 1;s++)
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++)
75 (unsigned char *)in += instride;
76 // the end points are copied as-is
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;
86 // copy the control point directly
87 for (c = 0;c < components;c++)
90 (unsigned char *)out += outstride;
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;
97 // copy the final control point
98 for (c = 0;c < components;c++)
100 //(unsigned char *)out += outstride;
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)
108 int finalwidth, finalheight, xstep, ystep, x, y, c;
111 // error out on various bogus conditions
112 if (xlevel < 0 || ylevel < 0 || xlevel > 16 || ylevel > 16 || cpwidth < 3 || cpheight < 3)
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++)
125 if (xlevel == 1 && ylevel == 0)
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);
131 if (xlevel == 0 && ylevel == 1)
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);
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++)
144 // subdivide in place in the destination buffer
145 while (xstep > 1 || ystep > 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;
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;
164 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
166 int c, x, y, outwidth, outheight, halfstep, xstep, ystep;
167 float prev, curr, next;
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)
181 halfstep = xstep >> 1;
182 for (y = 0;y < outheight;y += ystep)
184 for (c = 0;c < components;c++)
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)
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;
208 halfstep = ystep >> 1;
209 for (x = 0;x < outwidth;x += xstep)
211 for (c = 0;c < components;c++)
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)
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;
233 // flatten control points on X
234 for (y = 0;y < outheight;y += ystep)
236 for (c = 0;c < components;c++)
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)
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;;
251 // flatten control points on Y
252 for (x = 0;x < outwidth;x += xstep)
254 for (c = 0;c < components;c++)
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)
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;;
271 for (y = ystep;y < outheight - ystep;y += ystep)
273 for (c = 0;c < components;c++)
275 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])
278 outp[-halfstep * components] = (prev + curr) * 0.5f;
279 // flatten control point
280 outp[0] = (curr + (prev + next) * 0.5f) * 0.5f;
281 // next midpoint (only needed for end segment)
282 outp[halfstep * components] = (next + curr) * 0.5f;
290 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
292 int outwidth, outheight;
293 outwidth = ((cpwidth - 1) << xlevel) + 1;
294 outheight = ((cpheight - 1) << ylevel) + 1;
295 for (y = 0;y < cpheight;y++)
297 for (x = 0;x < cpwidth;x++)
299 for (c = 0;c < components;c++)
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++)
305 for (sx = 0;sx < expandx;sx++)
307 d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
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
323 float QuadraticSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
327 float deviation, squareddeviation, bestsquareddeviation;
328 bestsquareddeviation = 0;
329 for (y = 0;y < cpheight;y++)
331 for (x = 1;x < cpwidth-1;x++)
333 squareddeviation = 0;
334 for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
336 deviation = cp[0] * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
337 squareddeviation += deviation*deviation;
339 if (bestsquareddeviation < squareddeviation)
340 bestsquareddeviation = squareddeviation;
343 return (float)sqrt(bestsquareddeviation);
346 float QuadraticSplinePatchLargestDeviationOnY(int cpwidth, int cpheight, int components, const float *in)
350 float deviation, squareddeviation, bestsquareddeviation;
351 bestsquareddeviation = 0;
352 for (y = 1;y < cpheight-1;y++)
354 for (x = 0;x < cpwidth;x++)
356 squareddeviation = 0;
357 for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
359 deviation = cp[0] * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
360 squareddeviation += deviation*deviation;
362 if (bestsquareddeviation < squareddeviation)
363 bestsquareddeviation = squareddeviation;
366 return (float)sqrt(bestsquareddeviation);
369 int QuadraticSplinePatchSubdivisionLevelForDeviation(float deviation, float level1tolerance, int levellimit)
372 // count the automatic flatten step which reduces deviation by 50%
374 // count the levels to subdivide to come under the tolerance
375 for (level = 0;level < levellimit && deviation > level1tolerance;level++)
380 int QuadraticSplinePatchSubdivisionLevelOnX(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
382 return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnX(cpwidth, cpheight, components, in), level1tolerance, levellimit);
385 int QuadraticSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
387 return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnY(cpwidth, cpheight, components, in), level1tolerance, levellimit);
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[0] + (in[1] - in[0]) * 2 * t + (in[0] - in[1] + in[2] - in[1]) * t * t;