index a94dc71..cf7569d 100644 (file)
--- a/curves.c
+++ b/curves.c

-// this code written by Forest Hale, on 2003-08-23, and placed into public domain
-// this code deals with quadratic splines (minimum of 3 points), the same kind used in Quake3 maps.
+/*
+this code written by Forest Hale, on 2004-10-17, and placed into public domain
+this implements Quadratic BSpline surfaces as seen in Quake3 by id Software
+
+a small 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 cubic.
+
+example equations for 1-5 control point bsplines being sampled as t=0...1
+1: flat (0th dimension)
+o = a
+2: linear (1st dimension)
+o = a * (1 - t) + b * t
+o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
+4: cubic (bezier) bspline (3rd dimension)
+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
+5: quartic bspline (4th dimension)
+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
+
+arbitrary dimension bspline
+double factorial(int n)
+{
+       int i;
+       double f;
+       f = 1;
+       for (i = 1;i < n;i++)
+               f = f * i;
+       return f;
+}
+double bsplinesample(int dimensions, double t, double *param)
+{
+       double o = 0;
+       for (i = 0;i < dimensions + 1;i++)
+               o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
+       return o;
+}
+*/

-// 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).
+#include "quakedef.h"
+#include "mathlib.h"

#include <math.h>
#include "curves.h"
-#include "zone.h"

-#if 0
-void QuadraticSplineSubdivideFloat(int inpoints, int components, const float *in, int instride, float *out, int outstride)
+// Calculate number of resulting vertex rows/columns by given patch size and tesselation factor
+// tess=0 means that we reduce detalization of base 3x3 patches by removing middle row and column of vertices
+// "DimForTess" is "DIMension FOR TESSelation factor"
+// NB: tess=0 actually means that tess must be 0.5, but obviously it can't because it is of int type. (so "a*tess"-like code is replaced by "a/2" if tess=0)
+int Q3PatchDimForTess(int size, int tess)
+{
+       if (tess > 0)
+               return (size - 1) * tess + 1;
+       else if (tess == 0)
+               return (size - 1) / 2 + 1;
+       else
+               return 0; // Maybe warn about wrong tess here?
+}
+
+// usage:
+// to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
+// Q3PatchSubdivideFloat(3, sizeof(float), outvertices, 5, 5, sizeof(float), patchvertices, 4, 4);
+void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
{
-       int s;
-       // the input (control points) is read as a stream of points, and buffered
-       // by the cpprev, cpcurr, and cpnext variables (to allow subdivision in
-       // overlapping memory buffers, even subdivision in-place with pre-spaced
-       // control points in the buffer)
-       // the output (resulting curve) is written as a stream of points
-       // this subdivision is meant to be repeated until the desired flatness
-       // level is reached
-       if (components == 1 && instride == (int)sizeof(float) && outstride == instride)
+       int k, l, x, y, component, outputwidth = Q3PatchDimForTess(patchwidth, tesselationwidth);
+       float px, py, *v, a, b, c, *cp, temp;
+       int xmax = max(1, 2*tesselationwidth);
+       int ymax = max(1, 2*tesselationheight);
+
+       // iterate over the individual 3x3 quadratic spline surfaces one at a time
+       // expanding them to fill the output array (with some overlap to ensure
+       // the edges are filled)
+       for (k = 0;k < patchheight-1;k += 2)
{
-               // simple case, single component and no special stride
-               float cpprev0 = 0, cpcurr0 = 0, cpnext0;
-               cpnext0 = *in++;
-               for (s = 0;s < inpoints - 1;s++)
+               for (l = 0;l < patchwidth-1;l += 2)
{
-                       cpprev0 = cpcurr0;
-                       cpcurr0 = cpnext0;
-                       if (s < inpoints - 1)
-                               cpnext0 = *in++;
-                       if (s > 0)
-                       {
-                               // 50% flattened control point
-                               // cp1 = average(cp1, average(cp0, cp2));
-                               *out++ = (cpcurr0 + (cpprev0 + cpnext0) * 0.5f) * 0.5f;
-                       }
-                       else
+                       // set up control point pointers for quicker lookup later
+                       for (y = 0;y < 3;y++)
+                               for (x = 0;x < 3;x++)
+                                       cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
+                       // for each row...
+                       for (y = 0;y <= ymax;y++)
{
-                               // copy the control point directly
-                               *out++ = cpcurr0;
+                               // calculate control points for this row by collapsing the 3
+                               // rows of control points to one row using py
+                               py = (float)y / (float)ymax;
+                               // calculate quadratic spline weights for py
+                               a = ((1.0f - py) * (1.0f - py));
+                               b = ((1.0f - py) * (2.0f * py));
+                               c = ((       py) * (       py));
+                               for (component = 0;component < numcomponents;component++)
+                               {
+                                       temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
+                                       temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
+                                       temp[component] = cp[component] * a + cp[component] * b + cp[component] * c;
+                               }
+                               // fetch a pointer to the beginning of the output vertex row
+                               v = (float *)((unsigned char *)outputvertices + ((k * ymax / 2 + y) * outputwidth + l * xmax / 2) * outputstride);
+                               // for each column of the row...
+                               for (x = 0;x <= xmax;x++)
+                               {
+                                       // calculate point based on the row control points
+                                       px = (float)x / (float)xmax;
+                                       // calculate quadratic spline weights for px
+                                       // (could be precalculated)
+                                       a = ((1.0f - px) * (1.0f - px));
+                                       b = ((1.0f - px) * (2.0f * px));
+                                       c = ((       px) * (       px));
+                                       for (component = 0;component < numcomponents;component++)
+                                               v[component] = temp[component] * a + temp[component] * b + temp[component] * c;
+                                       // advance to next output vertex using outputstride
+                                       // (the next vertex may not be directly following this
+                                       // one, as this may be part of a larger structure)
+                                       v = (float *)((unsigned char *)v + outputstride);
+                               }
}
-                       // midpoint
-                       // mid = average(cp0, cp1);
-                       *out++ = (cpcurr0 + cpnext0) * 0.5f;
}
-               // copy the final control point
-               *out++ = cpnext0;
}
-       else
+#if 0
+       // enable this if you want results printed out
+       printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
+       for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
{
-               // multiple components or stride is used (complex case)
-               int c;
-               float cpprev, cpcurr, cpnext;
-               // check if there are too many components for the buffers
-               if (components > 1)
-               {
-                       // more components can be handled, but slowly, by calling self multiple times...
-                       for (c = 0;c < components;c++, in++, out++)
-                               QuadraticSplineSubdivideFloat(inpoints, 1, in, instride, out, outstride);
-                       return;
-               }
-               for (c = 0;c < components;c++)
-                       cpnext[c] = in[c];
-               (unsigned char *)in += instride;
-               for (s = 0;s < inpoints - 1;s++)
+               for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
{
-                       for (c = 0;c < components;c++)
-                               cpprev[c] = cpcurr[c];
-                       for (c = 0;c < components;c++)
-                               cpcurr[c] = cpnext[c];
-                       for (c = 0;c < components;c++)
-                               cpnext[c] = in[c];
-                       (unsigned char *)in += instride;
-                       // the end points are copied as-is
-                       if (s > 0)
-                       {
-                               // 50% flattened control point
-                               // cp1 = average(cp1, average(cp0, cp2));
-                               for (c = 0;c < components;c++)
-                                       out[c] = (cpcurr[c] + (cpprev[c] + cpnext[c]) * 0.5f) * 0.5f;
-                       }
-                       else
-                       {
-                               // copy the control point directly
-                               for (c = 0;c < components;c++)
-                                       out[c] = cpcurr[c];
-                       }
-                       (unsigned char *)out += outstride;
-                       // midpoint
-                       // mid = average(cp0, cp1);
-                       for (c = 0;c < components;c++)
-                               out[c] = (cpcurr[c] + cpnext[c]) * 0.5f;
-                       (unsigned char *)out += outstride;
+                       printf("(");
+                       for (component = 0;component < numcomponents;component++)
+                               printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
+                       printf(") ");
}
-               // copy the final control point
-               for (c = 0;c < components;c++)
-                       out[c] = cpnext[c];
-               //(unsigned char *)out += outstride;
+               printf("\n");
}
+       printf("}\n");
+#endif
}

-// note: out must have enough room!
-// (see finalwidth/finalheight calcs below)
-void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
+static int Q3PatchTesselation(float largestsquared3xcurvearea, float tolerance)
{
-       int finalwidth, finalheight, xstep, ystep, x, y, c;
-       float *o;
-
-       // error out on various bogus conditions
-       if (xlevel < 0 || ylevel < 0 || xlevel > 16 || ylevel > 16 || cpwidth < 3 || cpheight < 3)
-               return;
-
-       xstep = 1 << xlevel;
-       ystep = 1 << ylevel;
-       finalwidth = (cpwidth - 1) * xstep + 1;
-       finalheight = (cpheight - 1) * ystep + 1;
+       float f;
+       // f is actually a squared 2x curve area... so the formula had to be adjusted to give roughly the same subdivisions
+       f = pow(largestsquared3xcurvearea / 64.0f, 0.25f) / tolerance;
+       //if(f < 0.25) // VERY flat patches
+       if(f < 0.0001) // TOTALLY flat patches
+               return 0;
+       else if(f < 2)
+               return 1;
+       else
+               return (int) floor(log(f) / log(2.0f)) + 1;
+               // this is always at least 2
+               // maps [0.25..0.5[ to -1 (actually, 1 is returned)
+               // maps [0.5..1[ to 0 (actually, 1 is returned)
+               // maps [1..2[ to 1
+               // maps [2..4[ to 2
+               // maps [4..8[ to 4
+}

-       for (y = 0;y < finalheight;y++)
-               for (x = 0;x < finalwidth;x++)
-                       for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
-                               o[c] = 0;
+static float Squared3xCurveArea(const float *a, const float *control, const float *b, int components)
+{
+#if 0
+       // mimicing the old behaviour with the new code...

-       if (xlevel == 1 && ylevel == 0)
+       float deviation;
+       float quartercurvearea = 0;
+       int c;
+       for (c = 0;c < components;c++)
{
-               for (y = 0;y < finalheight;y++)
-                       QuadraticSplineSubdivideFloat(cpwidth, components, in + y * cpwidth * components, sizeof(float) * components, out + y * finalwidth * components, sizeof(float) * components);
-               return;
+               deviation = control[c] * 0.5f - a[c] * 0.25f - b[c] * 0.25f;
+               quartercurvearea += deviation*deviation;
}
-       if (xlevel == 0 && ylevel == 1)
+
+       // But as the new code now works on the squared 2x curve area, let's scale the value
+       return quartercurvearea * quartercurvearea * 64.0;
+
+#else
+       // ideally, we'd like the area between the spline a->control->b and the line a->b.
+       // but as this is hard to calculate, let's calculate an upper bound of it:
+       // the area of the triangle a->control->b->a.
+       //
+       // one can prove that the area of a quadratic spline = 2/3 * the area of
+       // the triangle of its control points!
+       // to do it, first prove it for the spline through (0,0), (1,1), (2,0)
+       // (which is a parabola) and then note that moving the control point
+       // left/right is just shearing and keeps the area of both the spline and
+       // the triangle invariant.
+       //
+       // why are we going for the spline area anyway?
+       // we know that:
+       //
+       //   the area between the spline and the line a->b is a measure of the
+       //   error of approximation of the spline by the line.
+       //
+       //   also, on circle-like or parabola-like curves, you easily get that the
+       //   double amount of line approximation segments reduces the error to its quarter
+       //   (also, easy to prove for splines by doing it for one specific one, and using
+       //   affine transforms to get all other splines)
+       //
+       // so...
+       //
+       // let's calculate the area! but we have to avoid the cross product, as
+       // components is not necessarily 3
+       //
+       // the area of a triangle spanned by vectors a and b is
+       //
+       // 0.5 * |a| |b| sin gamma
+       //
+       // now, cos gamma is
+       //
+       // a.b / (|a| |b|)
+       //
+       // so the area is
+       //
+       // 0.5 * sqrt(|a|^2 |b|^2 - (a.b)^2)
+       int c;
+       float aa = 0, bb = 0, ab = 0;
+       for (c = 0;c < components;c++)
{
-               for (x = 0;x < finalwidth;x++)
-                       QuadraticSplineSubdivideFloat(cpheight, components, in + x * components, sizeof(float) * cpwidth * components, out + x * components, sizeof(float) * finalwidth * components);
-               return;
+               float xa = a[c] - control[c];
+               float xb = b[c] - control[c];
+               aa += xa * xa;
+               ab += xa * xb;
+               bb += xb * xb;
}
+       // area is 0.5 * sqrt(aa*bb - ab*ab)
+       // 2x TRIANGLE area is sqrt(aa*bb - ab*ab)
+       // 3x CURVE area is sqrt(aa*bb - ab*ab)
+       return aa * bb - ab * ab;
+#endif
+}

-       // copy control points into correct positions in destination buffer
-       for (y = 0;y < finalheight;y += ystep)
-               for (x = 0;x < finalwidth;x += xstep)
-                       for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
-                               o[c] = *in++;
-
-       // subdivide in place in the destination buffer
-       while (xstep > 1 || ystep > 1)
+// returns how much tesselation of each segment is needed to remain under tolerance
+int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
+{
+       int x, y;
+       const float *patch;
+       float squared3xcurvearea, largestsquared3xcurvearea;
+       largestsquared3xcurvearea = 0;
+       for (y = 0;y < patchheight;y++)
{
-               if (xstep > 1)
+               for (x = 0;x < patchwidth-1;x += 2)
{
-                       xstep >>= 1;
-                       for (y = 0;y < finalheight;y += ystep)
-                               QuadraticSplineSubdivideFloat(cpwidth, components, out + y * finalwidth * components, sizeof(float) * xstep * 2 * components, out + y * finalwidth * components, sizeof(float) * xstep * components);
-                       cpwidth = (cpwidth - 1) * 2 + 1;
+                       patch = in + ((y * patchwidth) + x) * components;
+                       squared3xcurvearea = Squared3xCurveArea(&patch, &patch[components], &patch[2*components], components);
+                       if (largestsquared3xcurvearea < squared3xcurvearea)
+                               largestsquared3xcurvearea = squared3xcurvearea;
}
-               if (ystep > 1)
+       }
+       return Q3PatchTesselation(largestsquared3xcurvearea, tolerance);
+}
+
+// returns how much tesselation of each segment is needed to remain under tolerance
+int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
+{
+       int x, y;
+       const float *patch;
+       float squared3xcurvearea, largestsquared3xcurvearea;
+       largestsquared3xcurvearea = 0;
+       for (y = 0;y < patchheight-1;y += 2)
+       {
+               for (x = 0;x < patchwidth;x++)
{
-                       ystep >>= 1;
-                       for (x = 0;x < finalwidth;x += xstep)
-                               QuadraticSplineSubdivideFloat(cpheight, components, out + x * components, sizeof(float) * ystep * 2 * finalwidth * components, out + x * components, sizeof(float) * ystep * finalwidth * components);
-                       cpheight = (cpheight - 1) * 2 + 1;
+                       patch = in + ((y * patchwidth) + x) * components;
+                       squared3xcurvearea = Squared3xCurveArea(&patch, &patch[patchwidth*components], &patch[2*patchwidth*components], components);
+                       if (largestsquared3xcurvearea < squared3xcurvearea)
+                               largestsquared3xcurvearea = squared3xcurvearea;
}
}
+       return Q3PatchTesselation(largestsquared3xcurvearea, tolerance);
}
-#elif 1
-void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
+
+// Find an equal vertex in array. Check only vertices with odd X and Y
+static int FindEqualOddVertexInArray(int numcomponents, float *vertex, float *vertices, int width, int height)
{
-       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)
+       int x, y, j;
+       for (y=0; y<height; y+=2)
{
-               if (xstep >= ystep)
+               for (x=0; x<width; x+=2)
{
-                       // subdivide on X
-                       halfstep = xstep >> 1;
-                       for (y = 0;y < outheight;y += ystep)
-                       {
-                               for (c = 0;c < components;c++)
+                       qboolean found = true;
+                       for (j=0; j<numcomponents; j++)
+                               if (fabs(*(vertex+j) - *(vertices+j)) > 0.05)
+                               // div0: this is notably smaller than the smallest radiant grid
+                               // but large enough so we don't need to get scared of roundoff
+                               // errors
{
-                                       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;
-                                       }
+                                       found = false;
+                                       break;
}
-                       }
-                       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;
+                       if(found)
+                               return y*width+x;
+                       vertices += numcomponents*2;
}
+               vertices += numcomponents*(width-1);
}
-       // flatten control points on X
-       for (y = 0;y < outheight;y += ystep)
+       return -1;
+}
+
+#define SIDE_INVALID -1
+#define SIDE_X 0
+#define SIDE_Y 1
+
+static int GetSide(int p1, int p2, int width, int height, int *pointdist)
+{
+       int x1 = p1 % width, y1 = p1 / width;
+       int x2 = p2 % width, y2 = p2 / width;
+       if (p1 < 0 || p2 < 0)
+               return SIDE_INVALID;
+       if (x1 == x2)
{
-               for (c = 0;c < components;c++)
+               if (y1 != y2)
{
-                       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;;
-                       }
+                       *pointdist = abs(y2 - y1);
+                       return SIDE_Y;
}
+               else
+                       return SIDE_INVALID;
}
-       // flatten control points on Y
-       for (x = 0;x < outwidth;x += xstep)
+       else if (y1 == y2)
{
-               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;;
-                       }
-               }
+               *pointdist = abs(x2 - x1);
+               return SIDE_X;
}
+       else
+               return SIDE_INVALID;
+}

-       /*
-       for (y = ystep;y < outheight - ystep;y += ystep)
-       {
-               for (c = 0;c < components;c++)
+// Increase tesselation of one of two touching patches to make a seamless connection between them
+// Returns 0 in case if patches were not modified, otherwise 1
+int Q3PatchAdjustTesselation(int numcomponents, patchinfo_t *patch1, float *patchvertices1, patchinfo_t *patch2, float *patchvertices2)
+{
+       // what we are doing here is:
+       //   we take for each corner of one patch
+       //   and check if the other patch contains that corner
+       //   once we have a pair of such matches
+
+       struct {int id1,id2;} commonverts;
+       int i, j, k, side1, side2, *tess1, *tess2;
+       int dist1 = 0, dist2 = 0;
+       qboolean modified = false;
+
+       // Potential paired vertices (corners of the first patch)
+       commonverts.id1 = 0;
+       commonverts.id1 = patch1->xsize-1;
+       commonverts.id1 = patch1->xsize*(patch1->ysize-1);
+       commonverts.id1 = patch1->xsize*patch1->ysize-1;
+       for (i=0;i<4;++i)
+               commonverts[i].id2 = FindEqualOddVertexInArray(numcomponents, patchvertices1+numcomponents*commonverts[i].id1, patchvertices2, patch2->xsize, patch2->ysize);
+
+       // Corners of the second patch
+       commonverts.id2 = 0;
+       commonverts.id2 = patch2->xsize-1;
+       commonverts.id2 = patch2->xsize*(patch2->ysize-1);
+       commonverts.id2 = patch2->xsize*patch2->ysize-1;
+       for (i=4;i<8;++i)
+               commonverts[i].id1 = FindEqualOddVertexInArray(numcomponents, patchvertices2+numcomponents*commonverts[i].id2, patchvertices1, patch1->xsize, patch1->ysize);
+
+       for (i=0;i<8;++i)
+               for (j=i+1;j<8;++j)
{
-                       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])
+                       side1 = GetSide(commonverts[i].id1,commonverts[j].id1,patch1->xsize,patch1->ysize,&dist1);
+                       side2 = GetSide(commonverts[i].id2,commonverts[j].id2,patch2->xsize,patch2->ysize,&dist2);
+
+                       if (side1 == SIDE_INVALID || side2 == SIDE_INVALID)
+                               continue;
+
+                       if(dist1 != dist2)
{
-                               // midpoint
-                               outp[-halfstep * components] = (prev + curr) * 0.5f;
-                               // flatten control point
-                               outp = (curr + (prev + next) * 0.5f) * 0.5f;
-                               // next midpoint (only needed for end segment)
-                               outp[halfstep * components] = (next + curr) * 0.5f;
+                               // no patch welding if the resolutions mismatch
+                               continue;
}
-               }
-       }
-       */
-}
-#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++)
+
+                       // Update every lod level
+                       for (k=0;k<PATCH_LODS_NUM;++k)
{
-                               inp = in + (y * cpwidth + x) * components + c;
-                               outp = out + ((y<<ylevel) * outwidth + (x<<xlevel)) * components + c;
-                               for (sy = 0;sy < expandy;sy++)
+                               tess1 = side1 == SIDE_X ? &patch1->lods[k].xtess : &patch1->lods[k].ytess;
+                               tess2 = side2 == SIDE_X ? &patch2->lods[k].xtess : &patch2->lods[k].ytess;
+                               if (*tess1 != *tess2)
{
-                                       for (sx = 0;sx < expandx;sx++)
-                                       {
-                                               d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
-                                       }
+                                       if (*tess1 < *tess2)
+                                               *tess1 = *tess2;
+                                       else
+                                               *tess2 = *tess1;
+                                       modified = true;
}
}
}
-       }
+
+       return modified;
}
-#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
-*/
+#undef SIDE_INVALID
+#undef SIDE_X
+#undef SIDE_Y

-float QuadraticSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
+// calculates elements for a grid of vertices
+// (such as those produced by Q3PatchTesselate)
+// (note: width and height are the actual vertex size, this produces
+// (width-1)*(height-1)*2 triangles, 3 elements each)
+void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
{
-       int c, x, y;
-       const float *cp;
-       float deviation, squareddeviation, bestsquareddeviation;
-       bestsquareddeviation = 0;
-       for (y = 0;y < cpheight;y++)
+       int x, y, row0, row1;
+       for (y = 0;y < height - 1;y++)
{
-               for (x = 1;x < cpwidth-1;x++)
+               if(y % 2)
{
-                       squareddeviation = 0;
-                       for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
+                       // swap the triangle order in odd rows as optimization for collision stride
+                       row0 = firstvertex + (y + 0) * width + width - 2;
+                       row1 = firstvertex + (y + 1) * width + width - 2;
+                       for (x = 0;x < width - 1;x++)
{
-                               deviation = cp * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
-                               squareddeviation += deviation*deviation;
+                               *elements++ = row1;
+                               *elements++ = row1 + 1;
+                               *elements++ = row0 + 1;
+                               *elements++ = row0;
+                               *elements++ = row1;
+                               *elements++ = row0 + 1;
+                               row0--;
+                               row1--;
}
-                       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++)
+               else
{
-                       squareddeviation = 0;
-                       for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
+                       row0 = firstvertex + (y + 0) * width;
+                       row1 = firstvertex + (y + 1) * width;
+                       for (x = 0;x < width - 1;x++)
{
-                               deviation = cp * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
-                               squareddeviation += deviation*deviation;
+                               *elements++ = row0;
+                               *elements++ = row1;
+                               *elements++ = row0 + 1;
+                               *elements++ = row1;
+                               *elements++ = row1 + 1;
+                               *elements++ = row0 + 1;
+                               row0++;
+                               row1++;
}
-                       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)
-{
-}
-
-int QuadraticSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
-{
-}
-
-/*
-       d = a * (1 - 2 * t + t * t) + b * (2 * t - 2 * t * t) + c * t * t;
-       d = a * (1 + t * t + -2 * t) + b * (2 * t + -2 * t * t) + c * t * t;
-       d = a * 1 + a * t * t + a * -2 * t + b * 2 * t + b * -2 * t * t + c * t * t;
-       d = a * 1 + (a * t + a * -2) * t + (b * 2 + b * -2 * t) * t + (c * t) * t;
-       d = a + ((a * t + a * -2) + (b * 2 + b * -2 * t) + (c * t)) * t;
-       d = a + (a * (t - 2) + b * 2 + b * -2 * t + c * t) * t;
-       d = a + (a * (t - 2) + b * 2 + (b * -2 + c) * t) * t;
-       d = a + (a * (t - 2) + b * 2 + (c + b * -2) * t) * t;
-       d = a + a * (t - 2) * t + b * 2 * t + (c + b * -2) * t * t;
-       d = a * (1 + (t - 2) * t) + b * 2 * t + (c + b * -2) * t * t;
-       d = a * (1 + (t - 2) * t) + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a * 1 + a * (t - 2) * t + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a * 1 + a * t * t + a * -2 * t + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a * (1 - 2 * t + t * t) + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a * (1 - 2 * t) + a * t * t + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a + a * -2 * t + a * t * t + b * 2 * t + c * t * t + b * -2 * t * t;
-       d = a + a * -2 * t + a * t * t + b * 2 * t + b * -2 * t * t + c * t * t;
-       d = a + a * -2 * t + a * t * t + b * 2 * t + b * -2 * t * t + c * t * t;
-       d = a + a * -2 * t + b * 2 * t + b * -2 * t * t + a * t * t + c * t * t;
-       d = a + a * -2 * t + b * 2 * t + (a + c + b * -2) * t * t;
-       d = a + (a * -2 + b * 2) * t + (a + c + b * -2) * t * t;
-       d = a + ((a * -2 + b * 2) + (a + c + b * -2) * t) * t;
-       d = a + ((b + b - a - a) + (a + c - b - b) * t) * t;
-       d = a + (b + b - a - a) * t + (a + c - b - b) * t * t;
-       d = a + (b - a) * 2 * t + (a + c - b * 2) * t * t;
-       d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
-
-       d = in + (in - in) * 2 * t + (in - in + in - in) * t * t;
-*/
-