curves.c

   1
   2 /*
   3 this code written by Forest Hale, on 2004-10-17, and placed into public domain
   4 this implements Quadratic BSpline surfaces as seen in Quake3 by id Software
   5
   6 a small rant on misuse of the name 'bezier': many people seem to think that
   7 bezier is a generic term for splines, but it is not, it is a term for a
   8 specific type of bspline (4 control points, cubic bspline), bsplines are the
   9 generalization of the bezier spline to support dimensions other than cubic.
  10
  11 example equations for 1-5 control point bsplines being sampled as t=0...1
  12 1: flat (0th dimension)
  13 o = a
  14 2: linear (1st dimension)
  15 o = a * (1 - t) + b * t
  16 3: quadratic bspline (2nd dimension)
  17 o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
  18 4: cubic (bezier) bspline (3rd dimension)
  19 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
  20 5: quartic bspline (4th dimension)
  21 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
  22
  23 arbitrary dimension bspline
  24 double factorial(int n)
  25 {
  26         int i;
  27         double f;
  28         f = 1;
  29         for (i = 1;i < n;i++)
  30                 f = f * i;
  31         return f;
  32 }
  33 double bsplinesample(int dimensions, double t, double *param)
  34 {
  35         double o = 0;
  36         for (i = 0;i < dimensions + 1;i++)
  37                 o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
  38         return o;
  39 }
  40 */
  41
  42 #include "quakedef.h"
  43 #include "mathlib.h"
  44
  45 #include <math.h>
  46 #include "curves.h"
  47
  48 // Calculate number of resulting vertex rows/columns by given patch size and tesselation factor
  49 // tess=0 means that we reduce detalization of base 3x3 patches by removing middle row and column of vertices
  50 // "DimForTess" is "DIMension FOR TESSelation factor"
  51 // 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)
  52 int Q3PatchDimForTess(int size, int tess)
  53 {
  54         if (tess > 0)
  55                 return (size - 1) * tess + 1;
  56         else if (tess == 0)
  57                 return (size - 1) / 2 + 1;
  58         else
  59                 return 0; // Maybe warn about wrong tess here?
  60 }
  61
  62 // usage:
  63 // to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
  64 // Q3PatchSubdivideFloat(3, sizeof(float[3]), outvertices, 5, 5, sizeof(float[3]), patchvertices, 4, 4);
  65 void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
  66 {
  67         int k, l, x, y, component, outputwidth = Q3PatchDimForTess(patchwidth, tesselationwidth);
  68         float px, py, *v, a, b, c, *cp[3][3], temp[3][64];
  69         int xmax = max(1, 2*tesselationwidth);
  70         int ymax = max(1, 2*tesselationheight);
  71
  72         // iterate over the individual 3x3 quadratic spline surfaces one at a time
  73         // expanding them to fill the output array (with some overlap to ensure
  74         // the edges are filled)
  75         for (k = 0;k < patchheight-1;k += 2)
  76         {
  77                 for (l = 0;l < patchwidth-1;l += 2)
  78                 {
  79                         // set up control point pointers for quicker lookup later
  80                         for (y = 0;y < 3;y++)
  81                                 for (x = 0;x < 3;x++)
  82                                         cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
  83                         // for each row...
  84                         for (y = 0;y <= ymax;y++)
  85                         {
  86                                 // calculate control points for this row by collapsing the 3
  87                                 // rows of control points to one row using py
  88                                 py = (float)y / (float)ymax;
  89                                 // calculate quadratic spline weights for py
  90                                 a = ((1.0f - py) * (1.0f - py));
  91                                 b = ((1.0f - py) * (2.0f * py));
  92                                 c = ((       py) * (       py));
  93                                 for (component = 0;component < numcomponents;component++)
  94                                 {
  95                                         temp[0][component] = cp[0][0][component] * a + cp[1][0][component] * b + cp[2][0][component] * c;
  96                                         temp[1][component] = cp[0][1][component] * a + cp[1][1][component] * b + cp[2][1][component] * c;
  97                                         temp[2][component] = cp[0][2][component] * a + cp[1][2][component] * b + cp[2][2][component] * c;
  98                                 }
  99                                 // fetch a pointer to the beginning of the output vertex row
 100                                 v = (float *)((unsigned char *)outputvertices + ((k * ymax / 2 + y) * outputwidth + l * xmax / 2) * outputstride);
 101                                 // for each column of the row...
 102                                 for (x = 0;x <= xmax;x++)
 103                                 {
 104                                         // calculate point based on the row control points
 105                                         px = (float)x / (float)xmax;
 106                                         // calculate quadratic spline weights for px
 107                                         // (could be precalculated)
 108                                         a = ((1.0f - px) * (1.0f - px));
 109                                         b = ((1.0f - px) * (2.0f * px));
 110                                         c = ((       px) * (       px));
 111                                         for (component = 0;component < numcomponents;component++)
 112                                                 v[component] = temp[0][component] * a + temp[1][component] * b + temp[2][component] * c;
 113                                         // advance to next output vertex using outputstride
 114                                         // (the next vertex may not be directly following this
 115                                         // one, as this may be part of a larger structure)
 116                                         v = (float *)((unsigned char *)v + outputstride);
 117                                 }
 118                         }
 119                 }
 120         }
 121 #if 0
 122         // enable this if you want results printed out
 123         printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
 124         for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
 125         {
 126                 for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
 127                 {
 128                         printf("(");
 129                         for (component = 0;component < numcomponents;component++)
 130                                 printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
 131                         printf(") ");
 132                 }
 133                 printf("\n");
 134         }
 135         printf("}\n");
 136 #endif
 137 }
 138
 139 static int Q3PatchTesselation(float largestsquared2xcurvearea, float tolerance)
 140 {
 141         float f;
 142         // f is actually a squared 2x curve area... so the formula had to be adjusted to give roughly the same subdivisions
 143         f = pow(largestsquared2xcurvearea / 64.0, 0.25) / tolerance;
 144         //if(f < 0.25) // VERY flat patches
 145         if(f < 0.0001) // TOTALLY flat patches
 146                 return 0;
 147         else if(f < 2)
 148                 return 1;
 149         else
 150                 return (int) floor(log(f) / log(2)) + 1;
 151                 // this is always at least 2
 152                 // maps [0.25..0.5[ to -1 (actually, 1 is returned)
 153                 // maps [0.5..1[ to 0 (actually, 1 is returned)
 154                 // maps [1..2[ to 1
 155                 // maps [2..4[ to 2
 156                 // maps [4..8[ to 4
 157 }
 158
 159 float Squared2xCurveArea(const float *a, const float *control, const float *b, int components)
 160 {
 161 #if 0
 162         // mimicing the old behaviour with the new code...
 163
 164         float deviation;
 165         float quartercurvearea = 0;
 166         int c;
 167         for (c = 0;c < components;c++)
 168         {
 169                 deviation = control[c] * 0.5f - a[c] * 0.25f - b[c] * 0.25f;
 170                 quartercurvearea += deviation*deviation;
 171         }
 172
 173         // But as the new code now works on the squared 2x curve area, let's scale the value
 174         return quartercurvearea * quartercurvearea * 64.0;
 175
 176 #else
 177         // ideally, we'd like the area between the spline a->control->b and the line a->b.
 178         // but as this is hard to calculate, let's calculate an upper bound of it:
 179         // the area of the triangle a->control->b->a.
 180         //
 181         // why are we going for the spline area anyway?
 182         // we know that:
 183         //
 184         //   the area between the spline and the line a->b is a measure of the
 185         //   error of approximation of the spline by the line.
 186         //
 187         //   also, on circle-like or parabola-like curves, you easily get that the
 188         //   double amount of line approximation segments reduces the error to its quarter
 189         //
 190         // so...
 191         //
 192         // let's calculate the area! but we have to avoid the cross product, as
 193         // components is not necessarily 3
 194         //
 195         // the area of a triangle spanned by vectors a and b is
 196         //
 197         // 0.5 * |a| |b| sin gamma
 198         //
 199         // now, cos gamma is
 200         //
 201         // a.b / (|a| |b|)
 202         //
 203         // so the area is
 204         //
 205         // 0.5 * sqrt(|a|^2 |b|^2 - (a.b)^2)
 206         int c;
 207         float aa = 0, bb = 0, ab = 0;
 208         for (c = 0;c < components;c++)
 209         {
 210                 float xa = a[c] - control[c];
 211                 float xb = b[c] - control[c];
 212                 aa += xa * xa;
 213                 ab += xa * xb;
 214                 bb += xb * xb;
 215         }
 216         // area is 0.5 * sqrt(aa*bb - ab*ab)
 217         // 2x area is sqrt(aa*bb - ab*ab)
 218         return aa * bb - ab * ab;
 219 #endif
 220 }
 221
 222 // returns how much tesselation of each segment is needed to remain under tolerance
 223 int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
 224 {
 225         int x, y;
 226         const float *patch;
 227         float squared2xcurvearea, largestsquared2xcurvearea;
 228         largestsquared2xcurvearea = 0;
 229         for (y = 0;y < patchheight;y++)
 230         {
 231                 for (x = 0;x < patchwidth-1;x += 2)
 232                 {
 233                         patch = in + ((y * patchwidth) + x) * components;
 234                         squared2xcurvearea = Squared2xCurveArea(&patch[0], &patch[components], &patch[2*components], components);
 235                         if (largestsquared2xcurvearea < squared2xcurvearea)
 236                                 largestsquared2xcurvearea = squared2xcurvearea;
 237                 }
 238         }
 239         return Q3PatchTesselation(largestsquared2xcurvearea, tolerance);
 240 }
 241
 242 // returns how much tesselation of each segment is needed to remain under tolerance
 243 int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
 244 {
 245         int x, y;
 246         const float *patch;
 247         float squared2xcurvearea, largestsquared2xcurvearea;
 248         largestsquared2xcurvearea = 0;
 249         for (y = 0;y < patchheight-1;y += 2)
 250         {
 251                 for (x = 0;x < patchwidth;x++)
 252                 {
 253                         patch = in + ((y * patchwidth) + x) * components;
 254                         squared2xcurvearea = Squared2xCurveArea(&patch[0], &patch[patchwidth*components], &patch[2*patchwidth*components], components);
 255                         if (largestsquared2xcurvearea < squared2xcurvearea)
 256                                 largestsquared2xcurvearea = squared2xcurvearea;
 257                 }
 258         }
 259         return Q3PatchTesselation(largestsquared2xcurvearea, tolerance);
 260 }
 261
 262 // Find an equal vertex in array. Check only vertices with odd X and Y
 263 static int FindEqualOddVertexInArray(int numcomponents, float *vertex, float *vertices, int width, int height)
 264 {
 265         int x, y, j;
 266         for (y=0; y<height; y+=2)
 267         {
 268                 for (x=0; x<width; x+=2)
 269                 {
 270                         qboolean found = true;
 271                         for (j=0; j<numcomponents; j++)
 272                                 if (fabs(*(vertex+j) - *(vertices+j)) > 0.05)
 273                                 // div0: this is notably smaller than the smallest radiant grid
 274                                 // but large enough so we don't need to get scared of roundoff
 275                                 // errors
 276                                 {
 277                                         found = false;
 278                                         break;
 279                                 }
 280                         if(found)
 281                                 return y*width+x;
 282                         vertices += numcomponents*2;
 283                 }
 284                 vertices += numcomponents*(width-1);
 285         }
 286         return -1;
 287 }
 288
 289 #define SIDE_INVALID -1
 290 #define SIDE_X 0
 291 #define SIDE_Y 1
 292
 293 static int GetSide(int p1, int p2, int width, int height, int *pointdist)
 294 {
 295         int x1 = p1 % width, y1 = p1 / width;
 296         int x2 = p2 % width, y2 = p2 / width;
 297         if (p1 < 0 || p2 < 0)
 298                 return SIDE_INVALID;
 299         if (x1 == x2)
 300         {
 301                 if (y1 != y2)
 302                 {
 303                         *pointdist = abs(y2 - y1);
 304                         return SIDE_Y;
 305                 }
 306                 else
 307                         return SIDE_INVALID;
 308         }
 309         else if (y1 == y2)
 310         {
 311                 *pointdist = abs(x2 - x1);
 312                 return SIDE_X;
 313         }
 314         else
 315                 return SIDE_INVALID;
 316 }
 317
 318 // Increase tesselation of one of two touching patches to make a seamless connection between them
 319 // Returns 0 in case if patches were not modified, otherwise 1
 320 int Q3PatchAdjustTesselation(int numcomponents, patchinfo_t *patch1, float *patchvertices1, patchinfo_t *patch2, float *patchvertices2)
 321 {
 322         // what we are doing here is:
 323         //   we take for each corner of one patch
 324         //   and check if the other patch contains that corner
 325         //   once we have a pair of such matches
 326
 327         struct {int id1,id2;} commonverts[8];
 328         int i, j, k, side1, side2, *tess1, *tess2;
 329         int dist1, dist2;
 330         qboolean modified = false;
 331
 332         // Potential paired vertices (corners of the first patch)
 333         commonverts[0].id1 = 0;
 334         commonverts[1].id1 = patch1->xsize-1;
 335         commonverts[2].id1 = patch1->xsize*(patch1->ysize-1);
 336         commonverts[3].id1 = patch1->xsize*patch1->ysize-1;
 337         for (i=0;i<4;++i)
 338                 commonverts[i].id2 = FindEqualOddVertexInArray(numcomponents, patchvertices1+numcomponents*commonverts[i].id1, patchvertices2, patch2->xsize, patch2->ysize);
 339
 340         // Corners of the second patch
 341         commonverts[4].id2 = 0;
 342         commonverts[5].id2 = patch2->xsize-1;
 343         commonverts[6].id2 = patch2->xsize*(patch2->ysize-1);
 344         commonverts[7].id2 = patch2->xsize*patch2->ysize-1;
 345         for (i=4;i<8;++i)
 346                 commonverts[i].id1 = FindEqualOddVertexInArray(numcomponents, patchvertices2+numcomponents*commonverts[i].id2, patchvertices1, patch1->xsize, patch1->ysize);
 347
 348         for (i=0;i<8;++i)
 349                 for (j=i+1;j<8;++j)
 350                 {
 351                         side1 = GetSide(commonverts[i].id1,commonverts[j].id1,patch1->xsize,patch1->ysize,&dist1);
 352                         side2 = GetSide(commonverts[i].id2,commonverts[j].id2,patch2->xsize,patch2->ysize,&dist2);
 353
 354                         if (side1 == SIDE_INVALID || side2 == SIDE_INVALID)
 355                                 continue;
 356
 357                         if(dist1 != dist2)
 358                         {
 359                                 // no patch welding if the resolutions mismatch
 360                                 continue;
 361                         }
 362
 363                         // Update every lod level
 364                         for (k=0;k<PATCH_LODS_NUM;++k)
 365                         {
 366                                 tess1 = side1 == SIDE_X ? &patch1->lods[k].xtess : &patch1->lods[k].ytess;
 367                                 tess2 = side2 == SIDE_X ? &patch2->lods[k].xtess : &patch2->lods[k].ytess;
 368                                 if (*tess1 != *tess2)
 369                                 {
 370                                         if (*tess1 < *tess2)
 371                                                 *tess1 = *tess2;
 372                                         else
 373                                                 *tess2 = *tess1;
 374                                         modified = true;
 375                                 }
 376                         }
 377                 }
 378
 379         return modified;
 380 }
 381
 382 #undef SIDE_INVALID
 383 #undef SIDE_X
 384 #undef SIDE_Y
 385
 386 // calculates elements for a grid of vertices
 387 // (such as those produced by Q3PatchTesselate)
 388 // (note: width and height are the actual vertex size, this produces
 389 // (width-1)*(height-1)*2 triangles, 3 elements each)
 390 void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
 391 {
 392         int x, y, row0, row1;
 393         for (y = 0;y < height - 1;y++)
 394         {
 395                 row0 = firstvertex + (y + 0) * width;
 396                 row1 = firstvertex + (y + 1) * width;
 397                 for (x = 0;x < width - 1;x++)
 398                 {
 399                         *elements++ = row0;
 400                         *elements++ = row1;
 401                         *elements++ = row0 + 1;
 402                         *elements++ = row1;
 403                         *elements++ = row1 + 1;
 404                         *elements++ = row0 + 1;
 405                         row0++;
 406                         row1++;
 407                 }
 408         }
 409 }
 410