git.mutantstargoat.com Git - summerhack/blob - src/3dengfx/src/gfx/curves.cpp

   1 /*
   2 This file is part of the 3dengfx, realtime visualization system.
   3 Copyright (C) 2004, 2006 John Tsiombikas <nuclear@siggraph.org>
   4
   5 the 3dengfx library is free software; you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation; either version 2 of the License, or
   8 (at your option) any later version.
   9
  10 the 3dengfx library is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with the 3dengfx library; if not, write to the Free Software
  17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18 */
  19
  20 /* 3D Curves
  21  *
  22  * author: John Tsiombikas 2003
  23  * modified:
  24  *              John Tsiombikas 2004, 2006
  25  *              Mihalis Georgoulopoulos 2004
  26  */
  27
  28 #include <cstdio>
  29 #include <cstring>
  30 #include <cmath>
  31 #include <cctype>
  32 #include <cassert>
  33 #include "curves.hpp"
  34 #include "common/err_msg.h"
  35
  36 Curve::Curve() {
  37         arc_parametrize = false;
  38         ease_curve = 0;
  39         samples = 0;
  40         xform_cv = 0;
  41
  42         set_ease_sample_count(100);
  43 }
  44
  45 Curve::~Curve() {
  46         delete [] samples;
  47
  48 }
  49
  50 void Curve::set_arc_parametrization(bool state) {
  51         arc_parametrize = state;
  52 }
  53
  54 #define Param   0
  55 #define ArcLen  1
  56
  57 void Curve::sample_arc_lengths() {
  58         const int samplesPerSegment = 30;
  59         sample_count = get_segment_count() * samplesPerSegment;
  60
  61         arc_parametrize = false;        // to be able to interpolate with the original values
  62
  63         samples = new Vector2[sample_count];
  64         Vector3 prevpos;
  65         scalar_t step = 1.0f / (scalar_t)(sample_count-1);
  66         for(int i=0; i<sample_count; i++) {
  67                 scalar_t t = step * (scalar_t)i;
  68                 Vector3 pos = interpolate(t);
  69                 samples[i][Param] = t;
  70                 if(!i) {
  71                         samples[i][ArcLen] = 0.0f;
  72                 } else {
  73                         samples[i][ArcLen] = (pos - prevpos).length() + samples[i-1][ArcLen];
  74                 }
  75                 prevpos = pos;
  76         }
  77
  78         // normalize arc lenghts
  79         scalar_t maxlen = samples[sample_count-1][ArcLen];
  80         for(int i=0; i<sample_count; i++) {
  81                 samples[i][ArcLen] /= maxlen;
  82         }
  83
  84         arc_parametrize = true;
  85 }
  86
  87 static int binary_search(Vector2 *array, scalar_t key, int begin, int end) {
  88         int middle = begin + ((end - begin)>>1);
  89
  90         if(array[middle][ArcLen] == key) return middle;
  91         if(end == begin) return middle;
  92
  93         if(key < array[middle][ArcLen]) return binary_search(array, key, begin, middle);
  94         if(key > array[middle][ArcLen]) return binary_search(array, key, middle+1, end);
  95         return -1;      // just to make the compiler shut the fuck up
  96 }
  97
  98 scalar_t Curve::parametrize(scalar_t t) const {
  99         if(!samples) const_cast<Curve*>(this)->sample_arc_lengths();
 100
 101         int samplepos = binary_search(samples, t, 0, sample_count);
 102         scalar_t par = samples[samplepos][Param];
 103         scalar_t len = samples[samplepos][ArcLen];
 104
 105         // XXX: I can't remember the significance of this condition, I had xsmall_number here
 106         // previously and if t was 0.9999 it broke. I just changed the number blindly which fixed
 107         // the breakage but I should investigate further at some point.
 108         if((len - t) < 0.0005) return par;
 109
 110         if(len < t) {
 111                 if(!samplepos) return par;
 112                 scalar_t prevlen = samples[samplepos-1][ArcLen];
 113                 scalar_t prevpar = samples[samplepos-1][Param];
 114                 scalar_t p = (t - prevlen) / (len - prevlen);
 115                 return prevpar + (par - prevpar) * p;
 116         } else {
 117                 if(samplepos >= sample_count) return par;
 118                 scalar_t nextlen = samples[samplepos+1][ArcLen];
 119                 scalar_t nextpar = samples[samplepos+1][Param];
 120                 scalar_t p = (t - len) / (nextlen - len);
 121                 return par + (nextpar - par) * p;
 122         }
 123
 124         return par;     // not gonna happen
 125 }
 126
 127
 128 #define MIN(a, b) ((a) < (b) ? (a) : (b))
 129 #define MAX(a, b) ((a) > (b) ? (a) : (b))
 130
 131 scalar_t Curve::ease(scalar_t t) const {
 132         if(!ease_curve) return t;
 133
 134         const_cast<Curve*>(this)->ease_curve->set_arc_parametrization(true);
 135         scalar_t et = ease_curve->interpolate(t).y;
 136
 137         return MIN(MAX(et, 0.0f), 1.0f);
 138 }
 139
 140
 141 void Curve::add_control_point(const Vector3 &cp) {
 142         control_points.push_back(cp);
 143         delete [] samples;
 144         samples = 0;
 145 }
 146
 147 void Curve::remove_control_point(int index) {
 148         if(index < 0 || index >= control_points.size()) return;
 149
 150         ListNode<Vector3> *node = control_points.begin();
 151         for(int i=0; i<index; i++) {
 152                 node = node->next;
 153         }
 154
 155         control_points.erase(node);
 156 }
 157
 158 #define MIN(a, b)       ((a) < (b) ? (a) : (b))
 159 #define MAX(a, b)       ((a) > (b) ? (a) : (b))
 160
 161 Vector3 *Curve::get_control_point(int index) {
 162         index = MAX(0, MIN(index, control_points.size() - 1));
 163
 164         ListNode<Vector3> *node = control_points.begin();
 165         for(int i=0; i<index; i++) {
 166                 assert(node);
 167                 node = node->next;
 168         }
 169
 170         assert(node);
 171         return &node->data;
 172 }
 173
 174 int Curve::get_point_count() const {
 175         return control_points.size();
 176 }
 177
 178 void Curve::set_ease_curve(Curve *curve) {
 179         ease_curve = curve;
 180 }
 181
 182 void Curve::set_ease_sample_count(int count) {
 183         ease_sample_count = count;
 184         ease_step = (int)(1.0f / (scalar_t)ease_sample_count);
 185 }
 186
 187 Vector3 Curve::operator ()(scalar_t t) const {
 188         return xform_cv ? xform_cv(interpolate(t)) : interpolate(t);
 189 }
 190
 191 void Curve::set_xform_func(Vector3 (*func)(const Vector3&)) {
 192         xform_cv = func;
 193 }
 194
 195 ///////////////// B-Spline implementation ////////////////////
 196
 197 int BSpline::get_segment_count() const {
 198         return control_points.size() - 3;
 199 }
 200
 201 Vector3 BSpline::interpolate(scalar_t t) const {
 202         if(t > 1.0) t = 1.0;
 203         if(t < 0.0) t = 0.0;
 204
 205         if(control_points.size() < 4) return Vector3(0, 0, 0);
 206
 207         if(arc_parametrize) {
 208                 t = ease(parametrize(t));
 209         }
 210
 211         // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
 212         t = (scalar_t)(control_points.size() - 3) * t;
 213         int seg = (int)t;
 214         t -= (scalar_t)floor(t);
 215         if(seg >= get_segment_count()) {
 216                 seg = get_segment_count() - 1;
 217                 t = 1.0f;
 218         }
 219
 220         ListNode<Vector3> *iter = const_cast<BSpline*>(this)->control_points.begin();
 221         for(int i=0; i<seg; i++) iter = iter->next;
 222
 223         Vector3 Cp[4];
 224         for(int i=0; i<4; i++) {
 225         Cp[i] = iter->data;
 226                 iter = iter->next;
 227         }
 228
 229         Vector3 res;
 230         res.x = bspline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
 231         res.y = bspline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
 232         res.z = bspline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);
 233
 234         return res;
 235 }
 236
 237 //////////////// Catmull-Rom Spline implementation //////////////////
 238
 239 int CatmullRomSpline::get_segment_count() const {
 240         return control_points.size() - 1;
 241 }
 242
 243 Vector3 CatmullRomSpline::interpolate(scalar_t t) const {
 244         if(t > 1.0) t = 1.0;
 245         if(t < 0.0) t = 0.0;
 246
 247         if(control_points.size() < 2) return Vector3(0, 0, 0);
 248
 249         if(arc_parametrize) {
 250                 t = ease(parametrize(t));
 251         }
 252
 253         // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
 254         t = (scalar_t)(control_points.size() - 1) * t;
 255         int seg = (int)t;
 256         t -= (scalar_t)floor(t);
 257         if(seg >= get_segment_count()) {
 258                 seg = get_segment_count() - 1;
 259                 t = 1.0f;
 260         }
 261
 262         Vector3 cp[4];
 263         ListNode<Vector3> *iter = const_cast<CatmullRomSpline*>(this)->control_points.begin();
 264         for(int i=0; i<seg; i++) {
 265                 iter = iter->next;
 266         }
 267
 268         cp[1] = iter->data;
 269         cp[2] = iter->next->data;
 270
 271         if(!seg) {
 272                 cp[0] = cp[1];
 273         } else {
 274                 cp[0] = iter->prev->data;
 275         }
 276
 277         if(seg == control_points.size() - 2) {
 278                 cp[3] = cp[2];
 279         } else {
 280                 cp[3] = iter->next->next->data;
 281         }
 282
 283         Vector3 res;
 284         res.x = catmull_rom_spline(cp[0].x, cp[1].x, cp[2].x, cp[3].x, t);
 285         res.y = catmull_rom_spline(cp[0].y, cp[1].y, cp[2].y, cp[3].y, t);
 286         res.z = catmull_rom_spline(cp[0].z, cp[1].z, cp[2].z, cp[3].z, t);
 287
 288         return res;
 289 }
 290
 291
 292
 293 /* BezierSpline implementation - (MG) */
 294 int BezierSpline::get_segment_count() const
 295 {
 296         return control_points.size() / 4;
 297 }
 298
 299 Vector3 BezierSpline::interpolate(scalar_t t) const
 300 {
 301         if (!get_segment_count()) return Vector3(0, 0, 0);
 302
 303         if (arc_parametrize)
 304         {
 305                 t = ease(parametrize(t));
 306         }
 307
 308         t = (t * get_segment_count());
 309         int seg = (int) t;
 310         t -= (scalar_t) floor(t);
 311         if (seg >= get_segment_count())
 312         {
 313                 seg = get_segment_count() - 1;
 314                 t = 1.0f;
 315         }
 316
 317         seg *= 4;
 318         ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
 319         for (int i = 0; i < seg; i++) iter = iter->next;
 320
 321         Vector3 Cp[4];
 322         for (int i = 0; i < 4; i++)
 323         {
 324                 Cp[i] = iter->data;
 325                 iter = iter->next;
 326         }
 327
 328         // interpolate
 329         return bezier(Cp[seg], Cp[seg + 1], Cp[seg + 2], Cp[seg + 3], t);
 330 }
 331
 332 Vector3 BezierSpline::get_tangent(scalar_t t)
 333 {
 334         if (!get_segment_count()) return Vector3(0, 0, 0);
 335
 336         if (arc_parametrize)
 337         {
 338                 t = ease(parametrize(t));
 339         }
 340
 341         t = (t * get_segment_count());
 342         int seg = (int) t;
 343         t -= (scalar_t) floor(t);
 344         if (seg >= get_segment_count())
 345         {
 346                 seg = get_segment_count() - 1;
 347                 t = 1.0f;
 348         }
 349
 350         seg *= 4;
 351         ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
 352         for (int i = 0; i < seg; i++) iter = iter->next;
 353
 354         Vector3 Cp[4];
 355         for (int i = 0; i < 4; i++)
 356         {
 357                 Cp[i] = iter->data;
 358                 iter = iter->next;
 359         }
 360
 361         // interpolate
 362         return bezier_tangent(Cp[0], Cp[1], Cp[2], Cp[3], t);
 363 }
 364
 365 Vector3 BezierSpline::get_control_point(int i) const
 366 {
 367
 368         ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
 369         for (int j = 0; j < i; j++)
 370         {
 371                 if (!iter->next)
 372                 {
 373                         return Vector3(0, 0, 0);
 374                 }
 375                 iter = iter->next;
 376         }
 377
 378         return iter->data;
 379 }
 380
 381
 382 /* ------ polylines (JT) ------ */
 383 int PolyLine::get_segment_count() const {
 384         return control_points.size() - 1;
 385 }
 386
 387 Vector3 PolyLine::interpolate(scalar_t t) const {
 388         if(t > 1.0) t = 1.0;
 389         if(t < 0.0) t = 0.0;
 390
 391         if(control_points.size() < 2) return Vector3(0, 0, 0);
 392
 393         // TODO: check if this is reasonable for polylines.
 394         if(arc_parametrize) {
 395                 t = ease(parametrize(t));
 396         }
 397
 398         // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
 399         t = (scalar_t)(control_points.size() - 1) * t;
 400         int seg = (int)t;
 401         t -= (scalar_t)floor(t);
 402         if(seg >= get_segment_count()) {
 403                 seg = get_segment_count() - 1;
 404                 t = 1.0f;
 405         }
 406
 407         Vector3 cp[2];
 408         ListNode<Vector3> *iter = const_cast<PolyLine*>(this)->control_points.begin();
 409         for(int i=0; i<seg; i++) {
 410                 iter = iter->next;
 411         }
 412
 413         cp[0] = iter->data;
 414         cp[1] = iter->next->data;
 415
 416         return cp[0] + (cp[1] - cp[0]) * t;
 417 }
 418
 419
 420 bool save_curve(const char *fname, const Curve *curve) {
 421         FILE *fp = fopen(fname, "w");
 422         if(!fp) {
 423                 error("failed to save the curve %s", curve->name.c_str());
 424                 return false;
 425         }
 426
 427         fputs("curve_3dengfx\n", fp);
 428         fputs(curve->name.c_str(), fp);
 429         fputs("\n", fp);
 430
 431         if(dynamic_cast<const BSpline*>(curve)) {
 432                 fputs("bspline\n", fp);
 433         } else if(dynamic_cast<const CatmullRomSpline*>(curve)) {
 434                 fputs("catmullrom\n", fp);
 435         } else if(dynamic_cast<const BezierSpline*>(curve)) {
 436                 fputs("bezier\n", fp);
 437         } else if(dynamic_cast<const PolyLine*>(curve)) {
 438                 fputs("polyline\n", fp);
 439         } else {
 440                 error("unknown spline type, save failed %s", curve->name.c_str());
 441                 fclose(fp);
 442                 remove(fname);
 443                 return false;
 444         }
 445
 446         fprintf(fp, "%d\n", curve->control_points.size());
 447
 448         const ListNode<Vector3> *node = curve->control_points.begin();
 449         while(node) {
 450                 fprintf(fp, "%f %f %f\n", node->data.x, node->data.y, node->data.z);
 451                 node = node->next;
 452         }
 453
 454         fclose(fp);
 455         return true;
 456 }
 457
 458 Curve *load_curve(const char *fname) {
 459         FILE *fp = fopen(fname, "r");
 460         if(!fp) {
 461                 error("failed to open file %s", fname);
 462                 return 0;
 463         }
 464
 465         char buffer[256];
 466
 467         fgets(buffer, 256, fp);
 468         if(strcmp(buffer, "curve_3dengfx\n") != 0) {
 469                 error("load_curve failed, %s is not a curve file", fname);
 470                 fclose(fp);
 471                 return 0;
 472         }
 473
 474         Curve *curve;
 475
 476         fgets(buffer, 256, fp);
 477         std::string name = buffer;
 478
 479         fgets(buffer, 256, fp);
 480         if(!strcmp(buffer, "bspline\n")) {
 481                 curve = new BSpline;
 482         } else if(!strcmp(buffer, "catmullrom\n")) {
 483                 curve = new CatmullRomSpline;
 484         } else /*if(!strcmp(buffer, "bezier"))*/ {
 485                 error("unsupported curve type (%s) or not a curve file", buffer);
 486                 fclose(fp);
 487                 return 0;
 488         }
 489
 490         curve->name = name;
 491
 492         fgets(buffer, 256, fp);
 493         if(!isdigit(buffer[0])) {
 494                 error("load_curve failed, %s is not a valid curve file (count: %s)", fname, buffer);
 495                 delete curve;
 496                 fclose(fp);
 497                 return 0;
 498         }
 499         int count = atoi(buffer);
 500
 501         int failed = count;
 502         for(int i=0; i<count; i++, failed--) {
 503                 fgets(buffer, 256, fp);
 504                 if(!isdigit(buffer[0]) && buffer[0] != '.' && buffer[0] != '-') {
 505                         break;
 506                 }
 507                 float x = atof(buffer);
 508
 509                 char *ptr = strchr(buffer, ' ');
 510                 if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
 511                         break;
 512                 }
 513                 float y = atof(++ptr);
 514
 515                 ptr = strchr(ptr, ' ');
 516                 if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
 517                         break;
 518                 }
 519                 float z = atof(++ptr);
 520
 521                 curve->add_control_point(Vector3(x, y, z));
 522         }
 523
 524         fclose(fp);
 525
 526         if(failed) {
 527                 error("load_curve failed to read the data, %s is not a valid curve file", fname);
 528                 delete curve;
 529                 return 0;
 530         }
 531
 532         return curve;
 533 }