+#elif 1
+void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
+{
+ int c, x, y, outwidth, outheight, halfstep, xstep, ystep;
+ float prev, curr, next;
+ xstep = 1 << xlevel;
+ ystep = 1 << ylevel;
+ outwidth = ((cpwidth - 1) * xstep) + 1;
+ outheight = ((cpheight - 1) * ystep) + 1;
+ for (y = 0;y < cpheight;y++)
+ for (x = 0;x < cpwidth;x++)
+ for (c = 0;c < components;c++)
+ out[(y * ystep * outwidth + x * xstep) * components + c] = in[(y * cpwidth + x) * components + c];
+ while (xstep > 1 || ystep > 1)
+ {
+ if (xstep >= ystep)
+ {
+ // subdivide on X
+ halfstep = xstep >> 1;
+ for (y = 0;y < outheight;y += ystep)
+ {
+ for (c = 0;c < components;c++)
+ {
+ x = xstep;
+ // fetch first two control points
+ prev = out[(y * outwidth + (x - xstep)) * components + c];
+ curr = out[(y * outwidth + x) * components + c];
+ // create first midpoint
+ out[(y * outwidth + (x - halfstep)) * components + c] = (curr + prev) * 0.5f;
+ for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
+ {
+ // fetch next control point
+ next = out[(y * outwidth + (x + xstep)) * components + c];
+ // flatten central control point
+ out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
+ // create following midpoint
+ out[(y * outwidth + (x + halfstep)) * components + c] = (curr + next) * 0.5f;
+ }
+ }
+ }
+ xstep >>= 1;
+ }
+ else
+ {
+ // subdivide on Y
+ halfstep = ystep >> 1;
+ for (x = 0;x < outwidth;x += xstep)
+ {
+ for (c = 0;c < components;c++)
+ {
+ y = ystep;
+ // fetch first two control points
+ prev = out[((y - ystep) * outwidth + x) * components + c];
+ curr = out[(y * outwidth + x) * components + c];
+ // create first midpoint
+ out[((y - halfstep) * outwidth + x) * components + c] = (curr + prev) * 0.5f;
+ for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
+ {
+ // fetch next control point
+ next = out[((y + ystep) * outwidth + x) * components + c];
+ // flatten central control point
+ out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;;
+ // create following midpoint
+ out[((y + halfstep) * outwidth + x) * components + c] = (curr + next) * 0.5f;
+ }
+ }
+ }
+ ystep >>= 1;
+ }
+ }
+ // flatten control points on X
+ for (y = 0;y < outheight;y += ystep)
+ {
+ for (c = 0;c < components;c++)
+ {
+ x = xstep;
+ // fetch first two control points
+ prev = out[(y * outwidth + (x - xstep)) * components + c];
+ curr = out[(y * outwidth + x) * components + c];
+ for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
+ {
+ // fetch next control point
+ next = out[(y * outwidth + (x + xstep)) * components + c];
+ // flatten central control point
+ out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;;
+ }
+ }
+ }
+ // flatten control points on Y
+ for (x = 0;x < outwidth;x += xstep)
+ {
+ for (c = 0;c < components;c++)
+ {
+ y = ystep;
+ // fetch first two control points
+ prev = out[((y - ystep) * outwidth + x) * components + c];
+ curr = out[(y * outwidth + x) * components + c];
+ for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
+ {
+ // fetch next control point
+ next = out[((y + ystep) * outwidth + x) * components + c];
+ // flatten central control point
+ out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;;
+ }
+ }
+ }
+
+ /*
+ for (y = ystep;y < outheight - ystep;y += ystep)
+ {
+ for (c = 0;c < components;c++)
+ {
+ 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])
+ {
+ // midpoint
+ outp[-halfstep * components] = (prev + curr) * 0.5f;
+ // flatten control point
+ outp[0] = (curr + (prev + next) * 0.5f) * 0.5f;
+ // next midpoint (only needed for end segment)
+ outp[halfstep * components] = (next + curr) * 0.5f;
+ }
+ }
+ }
+ */
+}
+#else
+// unfinished code
+void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
+{
+ int outwidth, outheight;
+ outwidth = ((cpwidth - 1) << xlevel) + 1;
+ outheight = ((cpheight - 1) << ylevel) + 1;
+ for (y = 0;y < cpheight;y++)
+ {
+ for (x = 0;x < cpwidth;x++)
+ {
+ for (c = 0;c < components;c++)
+ {
+ inp = in + (y * cpwidth + x) * components + c;
+ outp = out + ((y<<ylevel) * outwidth + (x<<xlevel)) * components + c;
+ for (sy = 0;sy < expandy;sy++)
+ {
+ for (sx = 0;sx < expandx;sx++)
+ {
+ d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
+ }
+ }
+ }
+ }
+ }
+}
+#endif
+
+/*
+0.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 1.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 0.00000 deviation: 0.5
+0.00000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.00000 deviation: 0.125
+0.00000 ?.????? 0.25000 ?.????? 0.37500 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.37500 ?.????? 0.25000 ?.????? 0.00000 deviation: 0.03125
+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
+*/
+
+float QuadraticSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
+{
+ int c, x, y;
+ const float *cp;
+ float deviation, squareddeviation, bestsquareddeviation;
+ bestsquareddeviation = 0;
+ for (y = 0;y < cpheight;y++)
+ {
+ for (x = 1;x < cpwidth-1;x++)
+ {
+ squareddeviation = 0;
+ for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
+ {
+ deviation = cp[0] * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
+ squareddeviation += deviation*deviation;
+ }
+ if (bestsquareddeviation < squareddeviation)
+ bestsquareddeviation = squareddeviation;
+ }
+ }
+ return (float)sqrt(bestsquareddeviation);
+}
+
+float QuadraticSplinePatchLargestDeviationOnY(int cpwidth, int cpheight, int components, const float *in)
+{
+ int c, x, y;
+ const float *cp;
+ float deviation, squareddeviation, bestsquareddeviation;
+ bestsquareddeviation = 0;
+ for (y = 1;y < cpheight-1;y++)
+ {
+ for (x = 0;x < cpwidth;x++)
+ {
+ squareddeviation = 0;
+ for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
+ {
+ deviation = cp[0] * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
+ squareddeviation += deviation*deviation;
+ }
+ if (bestsquareddeviation < squareddeviation)
+ bestsquareddeviation = squareddeviation;
+ }
+ }
+ return (float)sqrt(bestsquareddeviation);
+}
+
+int QuadraticSplinePatchSubdivisionLevelForDeviation(float deviation, float level1tolerance, int levellimit)
+{
+ int level;
+ // count the automatic flatten step which reduces deviation by 50%
+ deviation *= 0.5f;
+ // count the levels to subdivide to come under the tolerance
+ for (level = 0;level < levellimit && deviation > level1tolerance;level++)
+ deviation *= 0.25f;
+ return level;
+}
+
+int QuadraticSplinePatchSubdivisionLevelOnX(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
+{
+ return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnX(cpwidth, cpheight, components, in), level1tolerance, levellimit);
+}
+
+int QuadraticSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
+{
+ return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnY(cpwidth, cpheight, components, in), level1tolerance, levellimit);
+}