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
  44 #include <math.h>
  45 #include "curves.h"
  46
  47 // usage:
  48 // to expand a 5x5 patch to 21x21 vertices (4x4 tesselation), one might use this call:
  49 // Q3PatchSubdivideFloat(3, sizeof(float[3]), outvertices, 5, 5, sizeof(float[3]), patchvertices, 4, 4);
  50 void Q3PatchTesselateFloat(int numcomponents, int outputstride, float *outputvertices, int patchwidth, int patchheight, int inputstride, float *patchvertices, int tesselationwidth, int tesselationheight)
  51 {
  52         int k, l, x, y, component, outputwidth = (patchwidth-1)*tesselationwidth+1;
  53         float px, py, *v, a, b, c, *cp[3][3], temp[3][64];
  54         // iterate over the individual 3x3 quadratic spline surfaces one at a time
  55         // expanding them to fill the output array (with some overlap to ensure
  56         // the edges are filled)
  57         for (k = 0;k < patchheight-1;k += 2)
  58         {
  59                 for (l = 0;l < patchwidth-1;l += 2)
  60                 {
  61                         // set up control point pointers for quicker lookup later
  62                         for (y = 0;y < 3;y++)
  63                                 for (x = 0;x < 3;x++)
  64                                         cp[y][x] = (float *)((unsigned char *)patchvertices + ((k+y)*patchwidth+(l+x)) * inputstride);
  65                         // for each row...
  66                         for (y = 0;y <= tesselationheight*2;y++)
  67                         {
  68                                 // calculate control points for this row by collapsing the 3
  69                                 // rows of control points to one row using py
  70                                 py = (float)y / (float)(tesselationheight*2);
  71                                 // calculate quadratic spline weights for py
  72                                 a = ((1.0f - py) * (1.0f - py));
  73                                 b = ((1.0f - py) * (2.0f * py));
  74                                 c = ((       py) * (       py));
  75                                 for (component = 0;component < numcomponents;component++)
  76                                 {
  77                                         temp[0][component] = cp[0][0][component] * a + cp[1][0][component] * b + cp[2][0][component] * c;
  78                                         temp[1][component] = cp[0][1][component] * a + cp[1][1][component] * b + cp[2][1][component] * c;
  79                                         temp[2][component] = cp[0][2][component] * a + cp[1][2][component] * b + cp[2][2][component] * c;
  80                                 }
  81                                 // fetch a pointer to the beginning of the output vertex row
  82                                 v = (float *)((unsigned char *)outputvertices + ((k * tesselationheight + y) * outputwidth + l * tesselationwidth) * outputstride);
  83                                 // for each column of the row...
  84                                 for (x = 0;x <= (tesselationwidth*2);x++)
  85                                 {
  86                                         // calculate point based on the row control points
  87                                         px = (float)x / (float)(tesselationwidth*2);
  88                                         // calculate quadratic spline weights for px
  89                                         // (could be precalculated)
  90                                         a = ((1.0f - px) * (1.0f - px));
  91                                         b = ((1.0f - px) * (2.0f * px));
  92                                         c = ((       px) * (       px));
  93                                         for (component = 0;component < numcomponents;component++)
  94                                                 v[component] = temp[0][component] * a + temp[1][component] * b + temp[2][component] * c;
  95                                         // advance to next output vertex using outputstride
  96                                         // (the next vertex may not be directly following this
  97                                         // one, as this may be part of a larger structure)
  98                                         v = (float *)((unsigned char *)v + outputstride);
  99                                 }
 100                         }
 101                 }
 102         }
 103 #if 0
 104         // enable this if you want results printed out
 105         printf("vertices[%i][%i] =\n{\n", (patchheight-1)*tesselationheight+1, (patchwidth-1)*tesselationwidth+1);
 106         for (y = 0;y < (patchheight-1)*tesselationheight+1;y++)
 107         {
 108                 for (x = 0;x < (patchwidth-1)*tesselationwidth+1;x++)
 109                 {
 110                         printf("(");
 111                         for (component = 0;component < numcomponents;component++)
 112                                 printf("%f ", outputvertices[(y*((patchwidth-1)*tesselationwidth+1)+x)*numcomponents+component]);
 113                         printf(") ");
 114                 }
 115                 printf("\n");
 116         }
 117         printf("}\n");
 118 #endif
 119 }
 120
 121 // returns how much tesselation of each segment is needed to remain under tolerance
 122 int Q3PatchTesselationOnX(int patchwidth, int patchheight, int components, const float *in, float tolerance)
 123 {
 124         int c, x, y;
 125         const float *patch;
 126         float deviation, squareddeviation, bestsquareddeviation;
 127         bestsquareddeviation = 0;
 128         for (y = 0;y < patchheight;y++)
 129         {
 130                 for (x = 0;x < patchwidth-1;x += 2)
 131                 {
 132                         squareddeviation = 0;
 133                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
 134                         {
 135                                 deviation = patch[components] * 0.5f - patch[0] * 0.25f - patch[2*components] * 0.25f;
 136                                 squareddeviation += deviation*deviation;
 137                         }
 138                         if (bestsquareddeviation < squareddeviation)
 139                                 bestsquareddeviation = squareddeviation;
 140                 }
 141         }
 142         return (int)floor(log(sqrt(bestsquareddeviation) / tolerance) / log(2)) + 1;
 143 }
 144
 145 // returns how much tesselation of each segment is needed to remain under tolerance
 146 int Q3PatchTesselationOnY(int patchwidth, int patchheight, int components, const float *in, float tolerance)
 147 {
 148         int c, x, y;
 149         const float *patch;
 150         float deviation, squareddeviation, bestsquareddeviation;
 151         bestsquareddeviation = 0;
 152         for (y = 0;y < patchheight-1;y += 2)
 153         {
 154                 for (x = 0;x < patchwidth;x++)
 155                 {
 156                         squareddeviation = 0;
 157                         for (c = 0, patch = in + ((y * patchwidth) + x) * components;c < components;c++, patch++)
 158                         {
 159                                 deviation = patch[patchwidth*components] * 0.5f - patch[0] * 0.25f - patch[2*patchwidth*components] * 0.25f;
 160                                 squareddeviation += deviation*deviation;
 161                         }
 162                         if (bestsquareddeviation < squareddeviation)
 163                                 bestsquareddeviation = squareddeviation;
 164                 }
 165         }
 166         return (int)floor(log(sqrt(bestsquareddeviation) / tolerance) / log(2)) + 1;
 167 }
 168
 169 // calculates elements for a grid of vertices
 170 // (such as those produced by Q3PatchTesselate)
 171 // (note: width and height are the actual vertex size, this produces
 172 //  (width-1)*(height-1)*2 triangles, 3 elements each)
 173 void Q3PatchTriangleElements(int *elements, int width, int height, int firstvertex)
 174 {
 175         int x, y, row0, row1;
 176         for (y = 0;y < height - 1;y++)
 177         {
 178                 row0 = firstvertex + (y + 0) * width;
 179                 row1 = firstvertex + (y + 1) * width;
 180                 for (x = 0;x < width - 1;x++)
 181                 {
 182                         *elements++ = row0;
 183                         *elements++ = row1;
 184                         *elements++ = row0 + 1;
 185                         *elements++ = row1;
 186                         *elements++ = row1 + 1;
 187                         *elements++ = row0 + 1;
 188                         row0++;
 189                         row1++;
 190                 }
 191         }
 192 }
 193