curves.c

   1
   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.
   4
   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 bspline (4 control points, cubic bspline), bsplines are the generalization of the bezier spline to support dimensions other than just cubic.
   6 // this implements Quadratic BSpline surfaces
   7
   8 #include <math.h>
   9 #include "curves.h"
  10 #include "zone.h"
  11
  12 #if 0
  13 void QuadraticSplineSubdivideFloat(int inpoints, int components, const float *in, int instride, float *out, int outstride)
  14 {
  15         int s;
  16         // the input (control points) is read as a stream of points, and buffered
  17         // by the cpprev, cpcurr, and cpnext variables (to allow subdivision in
  18         // overlapping memory buffers, even subdivision in-place with pre-spaced
  19         // control points in the buffer)
  20         // the output (resulting curve) is written as a stream of points
  21         // this subdivision is meant to be repeated until the desired flatness
  22         // level is reached
  23         if (components == 1 && instride == (int)sizeof(float) && outstride == instride)
  24         {
  25                 // simple case, single component and no special stride
  26                 float cpprev0 = 0, cpcurr0 = 0, cpnext0;
  27                 cpnext0 = *in++;
  28                 for (s = 0;s < inpoints - 1;s++)
  29                 {
  30                         cpprev0 = cpcurr0;
  31                         cpcurr0 = cpnext0;
  32                         if (s < inpoints - 1)
  33                                 cpnext0 = *in++;
  34                         if (s > 0)
  35                         {
  36                                 // 50% flattened control point
  37                                 // cp1 = average(cp1, average(cp0, cp2));
  38                                 *out++ = (cpcurr0 + (cpprev0 + cpnext0) * 0.5f) * 0.5f;
  39                         }
  40                         else
  41                         {
  42                                 // copy the control point directly
  43                                 *out++ = cpcurr0;
  44                         }
  45                         // midpoint
  46                         // mid = average(cp0, cp1);
  47                         *out++ = (cpcurr0 + cpnext0) * 0.5f;
  48                 }
  49                 // copy the final control point
  50                 *out++ = cpnext0;
  51         }
  52         else
  53         {
  54                 // multiple components or stride is used (complex case)
  55                 int c;
  56                 float cpprev[4], cpcurr[4], cpnext[4];
  57                 // check if there are too many components for the buffers
  58                 if (components > 1)
  59                 {
  60                         // more components can be handled, but slowly, by calling self multiple times...
  61                         for (c = 0;c < components;c++, in++, out++)
  62                                 QuadraticSplineSubdivideFloat(inpoints, 1, in, instride, out, outstride);
  63                         return;
  64                 }
  65                 for (c = 0;c < components;c++)
  66                         cpnext[c] = in[c];
  67                 (unsigned char *)in += instride;
  68                 for (s = 0;s < inpoints - 1;s++)
  69                 {
  70                         for (c = 0;c < components;c++)
  71                                 cpprev[c] = cpcurr[c];
  72                         for (c = 0;c < components;c++)
  73                                 cpcurr[c] = cpnext[c];
  74                         for (c = 0;c < components;c++)
  75                                 cpnext[c] = in[c];
  76                         (unsigned char *)in += instride;
  77                         // the end points are copied as-is
  78                         if (s > 0)
  79                         {
  80                                 // 50% flattened control point
  81                                 // cp1 = average(cp1, average(cp0, cp2));
  82                                 for (c = 0;c < components;c++)
  83                                         out[c] = (cpcurr[c] + (cpprev[c] + cpnext[c]) * 0.5f) * 0.5f;
  84                         }
  85                         else
  86                         {
  87                                 // copy the control point directly
  88                                 for (c = 0;c < components;c++)
  89                                         out[c] = cpcurr[c];
  90                         }
  91                         (unsigned char *)out += outstride;
  92                         // midpoint
  93                         // mid = average(cp0, cp1);
  94                         for (c = 0;c < components;c++)
  95                                 out[c] = (cpcurr[c] + cpnext[c]) * 0.5f;
  96                         (unsigned char *)out += outstride;
  97                 }
  98                 // copy the final control point
  99                 for (c = 0;c < components;c++)
 100                         out[c] = cpnext[c];
 101                 //(unsigned char *)out += outstride;
 102         }
 103 }
 104
 105 // note: out must have enough room!
 106 // (see finalwidth/finalheight calcs below)
 107 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
 108 {
 109         int finalwidth, finalheight, xstep, ystep, x, y, c;
 110         float *o;
 111
 112         // error out on various bogus conditions
 113         if (xlevel < 0 || ylevel < 0 || xlevel > 16 || ylevel > 16 || cpwidth < 3 || cpheight < 3)
 114                 return;
 115
 116         xstep = 1 << xlevel;
 117         ystep = 1 << ylevel;
 118         finalwidth = (cpwidth - 1) * xstep + 1;
 119         finalheight = (cpheight - 1) * ystep + 1;
 120
 121         for (y = 0;y < finalheight;y++)
 122                 for (x = 0;x < finalwidth;x++)
 123                         for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
 124                                 o[c] = 0;
 125
 126         if (xlevel == 1 && ylevel == 0)
 127         {
 128                 for (y = 0;y < finalheight;y++)
 129                         QuadraticSplineSubdivideFloat(cpwidth, components, in + y * cpwidth * components, sizeof(float) * components, out + y * finalwidth * components, sizeof(float) * components);
 130                 return;
 131         }
 132         if (xlevel == 0 && ylevel == 1)
 133         {
 134                 for (x = 0;x < finalwidth;x++)
 135                         QuadraticSplineSubdivideFloat(cpheight, components, in + x * components, sizeof(float) * cpwidth * components, out + x * components, sizeof(float) * finalwidth * components);
 136                 return;
 137         }
 138
 139         // copy control points into correct positions in destination buffer
 140         for (y = 0;y < finalheight;y += ystep)
 141                 for (x = 0;x < finalwidth;x += xstep)
 142                         for (c = 0, o = out + (y * finalwidth + x) * components;c < components;c++)
 143                                 o[c] = *in++;
 144
 145         // subdivide in place in the destination buffer
 146         while (xstep > 1 || ystep > 1)
 147         {
 148                 if (xstep > 1)
 149                 {
 150                         xstep >>= 1;
 151                         for (y = 0;y < finalheight;y += ystep)
 152                                 QuadraticSplineSubdivideFloat(cpwidth, components, out + y * finalwidth * components, sizeof(float) * xstep * 2 * components, out + y * finalwidth * components, sizeof(float) * xstep * components);
 153                         cpwidth = (cpwidth - 1) * 2 + 1;
 154                 }
 155                 if (ystep > 1)
 156                 {
 157                         ystep >>= 1;
 158                         for (x = 0;x < finalwidth;x += xstep)
 159                                 QuadraticSplineSubdivideFloat(cpheight, components, out + x * components, sizeof(float) * ystep * 2 * finalwidth * components, out + x * components, sizeof(float) * ystep * finalwidth * components);
 160                         cpheight = (cpheight - 1) * 2 + 1;
 161                 }
 162         }
 163 }
 164 #elif 1
 165 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
 166 {
 167         int c, x, y, outwidth, outheight, halfstep, xstep, ystep;
 168         float prev, curr, next;
 169         xstep = 1 << xlevel;
 170         ystep = 1 << ylevel;
 171         outwidth = ((cpwidth - 1) * xstep) + 1;
 172         outheight = ((cpheight - 1) * ystep) + 1;
 173         for (y = 0;y < cpheight;y++)
 174                 for (x = 0;x < cpwidth;x++)
 175                         for (c = 0;c < components;c++)
 176                                 out[(y * ystep * outwidth + x * xstep) * components + c] = in[(y * cpwidth + x) * components + c];
 177         while (xstep > 1 || ystep > 1)
 178         {
 179                 if (xstep >= ystep)
 180                 {
 181                         // subdivide on X
 182                         halfstep = xstep >> 1;
 183                         for (y = 0;y < outheight;y += ystep)
 184                         {
 185                                 for (c = 0;c < components;c++)
 186                                 {
 187                                         x = xstep;
 188                                         // fetch first two control points
 189                                         prev = out[(y * outwidth + (x - xstep)) * components + c];
 190                                         curr = out[(y * outwidth + x) * components + c];
 191                                         // create first midpoint
 192                                         out[(y * outwidth + (x - halfstep)) * components + c] = (curr + prev) * 0.5f;
 193                                         for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
 194                                         {
 195                                                 // fetch next control point
 196                                                 next = out[(y * outwidth + (x + xstep)) * components + c];
 197                                                 // flatten central control point
 198                                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
 199                                                 // create following midpoint
 200                                                 out[(y * outwidth + (x + halfstep)) * components + c] = (curr + next) * 0.5f;
 201                                         }
 202                                 }
 203                         }
 204                         xstep >>= 1;
 205                 }
 206                 else
 207                 {
 208                         // subdivide on Y
 209                         halfstep = ystep >> 1;
 210                         for (x = 0;x < outwidth;x += xstep)
 211                         {
 212                                 for (c = 0;c < components;c++)
 213                                 {
 214                                         y = ystep;
 215                                         // fetch first two control points
 216                                         prev = out[((y - ystep) * outwidth + x) * components + c];
 217                                         curr = out[(y * outwidth + x) * components + c];
 218                                         // create first midpoint
 219                                         out[((y - halfstep) * outwidth + x) * components + c] = (curr + prev) * 0.5f;
 220                                         for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
 221                                         {
 222                                                 // fetch next control point
 223                                                 next = out[((y + ystep) * outwidth + x) * components + c];
 224                                                 // flatten central control point
 225                                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
 226                                                 // create following midpoint
 227                                                 out[((y + halfstep) * outwidth + x) * components + c] = (curr + next) * 0.5f;
 228                                         }
 229                                 }
 230                         }
 231                         ystep >>= 1;
 232                 }
 233         }
 234         // flatten control points on X
 235         for (y = 0;y < outheight;y += ystep)
 236         {
 237                 for (c = 0;c < components;c++)
 238                 {
 239                         x = xstep;
 240                         // fetch first two control points
 241                         prev = out[(y * outwidth + (x - xstep)) * components + c];
 242                         curr = out[(y * outwidth + x) * components + c];
 243                         for (;x < outwidth - xstep;x += xstep, prev = curr, curr = next)
 244                         {
 245                                 // fetch next control point
 246                                 next = out[(y * outwidth + (x + xstep)) * components + c];
 247                                 // flatten central control point
 248                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
 249                         }
 250                 }
 251         }
 252         // flatten control points on Y
 253         for (x = 0;x < outwidth;x += xstep)
 254         {
 255                 for (c = 0;c < components;c++)
 256                 {
 257                         y = ystep;
 258                         // fetch first two control points
 259                         prev = out[((y - ystep) * outwidth + x) * components + c];
 260                         curr = out[(y * outwidth + x) * components + c];
 261                         for (;y < outheight - ystep;y += ystep, prev = curr, curr = next)
 262                         {
 263                                 // fetch next control point
 264                                 next = out[((y + ystep) * outwidth + x) * components + c];
 265                                 // flatten central control point
 266                                 out[(y * outwidth + x) * components + c] = (curr + (prev + next) * 0.5f) * 0.5f;
 267                         }
 268                 }
 269         }
 270
 271         /*
 272         for (y = ystep;y < outheight - ystep;y += ystep)
 273         {
 274                 for (c = 0;c < components;c++)
 275                 {
 276                         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])
 277                         {
 278                                 // midpoint
 279                                 outp[-halfstep * components] = (prev + curr) * 0.5f;
 280                                 // flatten control point
 281                                 outp[0] = (curr + (prev + next) * 0.5f) * 0.5f;
 282                                 // next midpoint (only needed for end segment)
 283                                 outp[halfstep * components] = (next + curr) * 0.5f;
 284                         }
 285                 }
 286         }
 287         */
 288 }
 289 #else
 290 // unfinished code
 291 void QuadraticSplinePatchSubdivideFloatBuffer(int cpwidth, int cpheight, int xlevel, int ylevel, int components, const float *in, float *out)
 292 {
 293         int outwidth, outheight;
 294         outwidth = ((cpwidth - 1) << xlevel) + 1;
 295         outheight = ((cpheight - 1) << ylevel) + 1;
 296         for (y = 0;y < cpheight;y++)
 297         {
 298                 for (x = 0;x < cpwidth;x++)
 299                 {
 300                         for (c = 0;c < components;c++)
 301                         {
 302                                 inp = in + (y * cpwidth + x) * components + c;
 303                                 outp = out + ((y<<ylevel) * outwidth + (x<<xlevel)) * components + c;
 304                                 for (sy = 0;sy < expandy;sy++)
 305                                 {
 306                                         for (sx = 0;sx < expandx;sx++)
 307                                         {
 308                                                 d = a + (b - a) * 2 * t + (a - b + c - b) * t * t;
 309                                         }
 310                                 }
 311                         }
 312                 }
 313         }
 314 }
 315 #endif
 316
 317 /*
 318 0.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 1.00000 ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? ?.????? 0.00000 deviation: 0.5
 319 0.00000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.50000 ?.????? ?.????? ?.????? 0.00000 deviation: 0.125
 320 0.00000 ?.????? 0.25000 ?.????? 0.37500 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.50000 ?.????? 0.37500 ?.????? 0.25000 ?.????? 0.00000 deviation: 0.03125
 321 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
 322 */
 323
 324 float QuadraticSplinePatchLargestDeviationOnX(int cpwidth, int cpheight, int components, const float *in)
 325 {
 326         int c, x, y;
 327         const float *cp;
 328         float deviation, squareddeviation, bestsquareddeviation;
 329         bestsquareddeviation = 0;
 330         for (y = 0;y < cpheight;y++)
 331         {
 332                 for (x = 1;x < cpwidth-1;x++)
 333                 {
 334                         squareddeviation = 0;
 335                         for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
 336                         {
 337                                 deviation = cp[0] * 0.5f - cp[-components] * 0.25f - cp[components] * 0.25f;
 338                                 squareddeviation += deviation*deviation;
 339                         }
 340                         if (bestsquareddeviation < squareddeviation)
 341                                 bestsquareddeviation = squareddeviation;
 342                 }
 343         }
 344         return (float)sqrt(bestsquareddeviation);
 345 }
 346
 347 float QuadraticSplinePatchLargestDeviationOnY(int cpwidth, int cpheight, int components, const float *in)
 348 {
 349         int c, x, y;
 350         const float *cp;
 351         float deviation, squareddeviation, bestsquareddeviation;
 352         bestsquareddeviation = 0;
 353         for (y = 1;y < cpheight-1;y++)
 354         {
 355                 for (x = 0;x < cpwidth;x++)
 356                 {
 357                         squareddeviation = 0;
 358                         for (c = 0, cp = in + ((y * cpwidth) + x) * components;c < components;c++, cp++)
 359                         {
 360                                 deviation = cp[0] * 0.5f - cp[-cpwidth * components] * 0.25f - cp[cpwidth * components] * 0.25f;
 361                                 squareddeviation += deviation*deviation;
 362                         }
 363                         if (bestsquareddeviation < squareddeviation)
 364                                 bestsquareddeviation = squareddeviation;
 365                 }
 366         }
 367         return (float)sqrt(bestsquareddeviation);
 368 }
 369
 370 int QuadraticSplinePatchSubdivisionLevelForDeviation(float deviation, float level1tolerance, int levellimit)
 371 {
 372         int level;
 373         // count the automatic flatten step which reduces deviation by 50%
 374         deviation *= 0.5f;
 375         // count the levels to subdivide to come under the tolerance
 376         for (level = 0;level < levellimit && deviation > level1tolerance;level++)
 377                 deviation *= 0.25f;
 378         return level;
 379 }
 380
 381 int QuadraticSplinePatchSubdivisionLevelOnX(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
 382 {
 383         return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnX(cpwidth, cpheight, components, in), level1tolerance, levellimit);
 384 }
 385
 386 int QuadraticSplinePatchSubdivisionLevelOnY(int cpwidth, int cpheight, int components, const float *in, float level1tolerance, int levellimit)
 387 {
 388         return QuadraticSplinePatchSubdivisionLevelForDeviation(QuadraticSplinePatchLargestDeviationOnY(cpwidth, cpheight, components, in), level1tolerance, levellimit);
 389 }
 390
 391 /*
 392         // 1: flat (0th dimension)
 393         o = a
 394         // 2: linear (1st dimension)
 395         o = a * (1 - t) + b * t
 396         // 3: quadratic bspline (2nd dimension)
 397         o = a * (1 - t) * (1 - t) + 2 * b * (1 - t) * t + c * t * t
 398         // 4: cubic (bezier) bspline (3rd dimension)
 399         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
 400         // 5: quartic bspline (4th dimension)
 401         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
 402
 403         // n: arbitrary dimension bspline
 404 double factorial(int n)
 405 {
 406         int i;
 407         double f;
 408         f = 1;
 409         for (i = 1;i < n;i++)
 410                 f = f * i;
 411         return f;
 412 }
 413 double bsplinesample(int dimensions, double t, double *param)
 414 {
 415         double o = 0;
 416         for (i = 0;i < dimensions + 1;i++)
 417                 o += param[i] * factorial(dimensions)/(factorial(i)*factorial(dimensions-i)) * pow(t, i) * pow(1 - t, dimensions - i);
 418 }
 419 */
 420