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