added 3dengfx into the repo, probably not the correct version for this
[summerhack] / src / 3dengfx / src / gfx / curves.cpp
diff --git a/src/3dengfx/src/gfx/curves.cpp b/src/3dengfx/src/gfx/curves.cpp
new file mode 100644 (file)
index 0000000..a71194f
--- /dev/null
@@ -0,0 +1,533 @@
+/*
+This file is part of the 3dengfx, realtime visualization system.
+Copyright (C) 2004, 2006 John Tsiombikas <nuclear@siggraph.org>
+
+the 3dengfx library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+the 3dengfx library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the 3dengfx library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+/* 3D Curves
+ * 
+ * author: John Tsiombikas 2003
+ * modified:
+ *             John Tsiombikas 2004, 2006
+ *             Mihalis Georgoulopoulos 2004
+ */
+
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <cctype>
+#include <cassert>
+#include "curves.hpp"
+#include "common/err_msg.h"
+
+Curve::Curve() {
+       arc_parametrize = false;
+       ease_curve = 0;
+       samples = 0;
+       xform_cv = 0;
+
+       set_ease_sample_count(100);
+}
+
+Curve::~Curve() {
+       delete [] samples;
+       
+}
+
+void Curve::set_arc_parametrization(bool state) {
+       arc_parametrize = state;
+}
+
+#define Param  0
+#define ArcLen 1
+
+void Curve::sample_arc_lengths() {
+       const int samplesPerSegment = 30;
+       sample_count = get_segment_count() * samplesPerSegment;
+
+       arc_parametrize = false;        // to be able to interpolate with the original values
+
+       samples = new Vector2[sample_count];
+       Vector3 prevpos;
+       scalar_t step = 1.0f / (scalar_t)(sample_count-1);
+       for(int i=0; i<sample_count; i++) {
+               scalar_t t = step * (scalar_t)i;
+               Vector3 pos = interpolate(t);
+               samples[i][Param] = t;
+               if(!i) {
+                       samples[i][ArcLen] = 0.0f;
+               } else {
+                       samples[i][ArcLen] = (pos - prevpos).length() + samples[i-1][ArcLen];
+               }
+               prevpos = pos;
+       }
+
+       // normalize arc lenghts
+       scalar_t maxlen = samples[sample_count-1][ArcLen];
+       for(int i=0; i<sample_count; i++) {
+               samples[i][ArcLen] /= maxlen;
+       }
+
+       arc_parametrize = true;
+}
+
+static int binary_search(Vector2 *array, scalar_t key, int begin, int end) {
+       int middle = begin + ((end - begin)>>1);
+
+       if(array[middle][ArcLen] == key) return middle;
+       if(end == begin) return middle;
+
+       if(key < array[middle][ArcLen]) return binary_search(array, key, begin, middle);
+       if(key > array[middle][ArcLen]) return binary_search(array, key, middle+1, end);
+       return -1;      // just to make the compiler shut the fuck up
+}
+
+scalar_t Curve::parametrize(scalar_t t) const {
+       if(!samples) const_cast<Curve*>(this)->sample_arc_lengths();
+
+       int samplepos = binary_search(samples, t, 0, sample_count);
+       scalar_t par = samples[samplepos][Param];
+       scalar_t len = samples[samplepos][ArcLen];
+
+       // XXX: I can't remember the significance of this condition, I had xsmall_number here
+       // previously and if t was 0.9999 it broke. I just changed the number blindly which fixed
+       // the breakage but I should investigate further at some point.
+       if((len - t) < 0.0005) return par;
+
+       if(len < t) {
+               if(!samplepos) return par;
+               scalar_t prevlen = samples[samplepos-1][ArcLen];
+               scalar_t prevpar = samples[samplepos-1][Param];
+               scalar_t p = (t - prevlen) / (len - prevlen);
+               return prevpar + (par - prevpar) * p;
+       } else {
+               if(samplepos >= sample_count) return par;
+               scalar_t nextlen = samples[samplepos+1][ArcLen];
+               scalar_t nextpar = samples[samplepos+1][Param];
+               scalar_t p = (t - len) / (nextlen - len);
+               return par + (nextpar - par) * p;
+       }
+
+       return par;     // not gonna happen
+}
+
+
+#define MIN(a, b) ((a) < (b) ? (a) : (b))
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+
+scalar_t Curve::ease(scalar_t t) const {
+       if(!ease_curve) return t;
+
+       const_cast<Curve*>(this)->ease_curve->set_arc_parametrization(true);
+       scalar_t et = ease_curve->interpolate(t).y;
+
+       return MIN(MAX(et, 0.0f), 1.0f);
+}
+
+
+void Curve::add_control_point(const Vector3 &cp) {
+       control_points.push_back(cp);
+       delete [] samples;
+       samples = 0;
+}
+
+void Curve::remove_control_point(int index) {
+       if(index < 0 || index >= control_points.size()) return;
+
+       ListNode<Vector3> *node = control_points.begin();
+       for(int i=0; i<index; i++) {
+               node = node->next;
+       }
+
+       control_points.erase(node);
+}
+
+#define MIN(a, b)      ((a) < (b) ? (a) : (b))
+#define MAX(a, b)      ((a) > (b) ? (a) : (b))
+
+Vector3 *Curve::get_control_point(int index) {
+       index = MAX(0, MIN(index, control_points.size() - 1));
+
+       ListNode<Vector3> *node = control_points.begin();
+       for(int i=0; i<index; i++) {
+               assert(node);
+               node = node->next;
+       }
+       
+       assert(node);
+       return &node->data;
+}
+
+int Curve::get_point_count() const {
+       return control_points.size();
+}
+
+void Curve::set_ease_curve(Curve *curve) {
+       ease_curve = curve;
+}
+
+void Curve::set_ease_sample_count(int count) {
+       ease_sample_count = count;
+       ease_step = (int)(1.0f / (scalar_t)ease_sample_count);
+}
+
+Vector3 Curve::operator ()(scalar_t t) const {
+       return xform_cv ? xform_cv(interpolate(t)) : interpolate(t);
+}
+
+void Curve::set_xform_func(Vector3 (*func)(const Vector3&)) {
+       xform_cv = func;
+}
+
+///////////////// B-Spline implementation ////////////////////
+
+int BSpline::get_segment_count() const {
+       return control_points.size() - 3;
+}
+
+Vector3 BSpline::interpolate(scalar_t t) const {
+       if(t > 1.0) t = 1.0;
+       if(t < 0.0) t = 0.0;
+
+       if(control_points.size() < 4) return Vector3(0, 0, 0);
+
+       if(arc_parametrize) {
+               t = ease(parametrize(t));
+       }
+
+       // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
+       t = (scalar_t)(control_points.size() - 3) * t;
+       int seg = (int)t;
+       t -= (scalar_t)floor(t);
+       if(seg >= get_segment_count()) {
+               seg = get_segment_count() - 1;
+               t = 1.0f;
+       }
+       
+       ListNode<Vector3> *iter = const_cast<BSpline*>(this)->control_points.begin();
+       for(int i=0; i<seg; i++) iter = iter->next;
+
+       Vector3 Cp[4];
+       for(int i=0; i<4; i++) {
+        Cp[i] = iter->data;
+               iter = iter->next;
+       }
+
+       Vector3 res;
+       res.x = bspline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
+       res.y = bspline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
+       res.z = bspline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);
+
+       return res;
+}
+
+//////////////// Catmull-Rom Spline implementation //////////////////
+
+int CatmullRomSpline::get_segment_count() const {
+       return control_points.size() - 1;
+}
+
+Vector3 CatmullRomSpline::interpolate(scalar_t t) const {
+       if(t > 1.0) t = 1.0;
+       if(t < 0.0) t = 0.0;
+
+       if(control_points.size() < 2) return Vector3(0, 0, 0);
+
+       if(arc_parametrize) {
+               t = ease(parametrize(t));
+       }
+
+       // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
+       t = (scalar_t)(control_points.size() - 1) * t;
+       int seg = (int)t;
+       t -= (scalar_t)floor(t);
+       if(seg >= get_segment_count()) {
+               seg = get_segment_count() - 1;
+               t = 1.0f;
+       }
+
+       Vector3 cp[4];
+       ListNode<Vector3> *iter = const_cast<CatmullRomSpline*>(this)->control_points.begin();
+       for(int i=0; i<seg; i++) {
+               iter = iter->next;
+       }
+
+       cp[1] = iter->data;
+       cp[2] = iter->next->data;
+       
+       if(!seg) {
+               cp[0] = cp[1];
+       } else {
+               cp[0] = iter->prev->data;
+       }
+       
+       if(seg == control_points.size() - 2) {
+               cp[3] = cp[2];
+       } else {
+               cp[3] = iter->next->next->data;
+       }
+
+       Vector3 res;
+       res.x = catmull_rom_spline(cp[0].x, cp[1].x, cp[2].x, cp[3].x, t);
+       res.y = catmull_rom_spline(cp[0].y, cp[1].y, cp[2].y, cp[3].y, t);
+       res.z = catmull_rom_spline(cp[0].z, cp[1].z, cp[2].z, cp[3].z, t);
+
+       return res;
+}
+
+
+
+/* BezierSpline implementation - (MG) */
+int BezierSpline::get_segment_count() const
+{
+       return control_points.size() / 4;
+}
+
+Vector3 BezierSpline::interpolate(scalar_t t) const
+{
+       if (!get_segment_count()) return Vector3(0, 0, 0);
+
+       if (arc_parametrize)
+       {
+               t = ease(parametrize(t));
+       }
+       t = (t * get_segment_count());
+       int seg = (int) t;
+       t -= (scalar_t) floor(t);
+       if (seg >= get_segment_count())
+       {
+               seg = get_segment_count() - 1;
+               t = 1.0f;
+       }
+
+       seg *= 4;
+       ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+       for (int i = 0; i < seg; i++) iter = iter->next;
+       
+       Vector3 Cp[4];
+       for (int i = 0; i < 4; i++)
+       {
+               Cp[i] = iter->data;
+               iter = iter->next;
+       }
+       
+       // interpolate
+       return bezier(Cp[seg], Cp[seg + 1], Cp[seg + 2], Cp[seg + 3], t);
+}
+
+Vector3 BezierSpline::get_tangent(scalar_t t)
+{      
+       if (!get_segment_count()) return Vector3(0, 0, 0);
+
+       if (arc_parametrize)
+       {
+               t = ease(parametrize(t));
+       }
+       t = (t * get_segment_count());
+       int seg = (int) t;
+       t -= (scalar_t) floor(t);
+       if (seg >= get_segment_count())
+       {
+               seg = get_segment_count() - 1;
+               t = 1.0f;
+       }
+
+       seg *= 4;
+       ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+       for (int i = 0; i < seg; i++) iter = iter->next;
+       
+       Vector3 Cp[4];
+       for (int i = 0; i < 4; i++)
+       {
+               Cp[i] = iter->data;
+               iter = iter->next;
+       }
+       
+       // interpolate
+       return bezier_tangent(Cp[0], Cp[1], Cp[2], Cp[3], t);
+}
+
+Vector3 BezierSpline::get_control_point(int i) const
+{      
+       
+       ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+       for (int j = 0; j < i; j++) 
+       {
+               if (!iter->next)
+               {
+                       return Vector3(0, 0, 0);
+               }
+               iter = iter->next;
+       }
+       
+       return iter->data;
+}
+
+
+/* ------ polylines (JT) ------ */
+int PolyLine::get_segment_count() const {
+       return control_points.size() - 1;
+}
+
+Vector3 PolyLine::interpolate(scalar_t t) const {
+       if(t > 1.0) t = 1.0;
+       if(t < 0.0) t = 0.0;
+
+       if(control_points.size() < 2) return Vector3(0, 0, 0);
+
+       // TODO: check if this is reasonable for polylines.
+       if(arc_parametrize) {
+               t = ease(parametrize(t));
+       }
+
+       // find the appropriate segment of the spline that t lies and calculate the piecewise parameter
+       t = (scalar_t)(control_points.size() - 1) * t;
+       int seg = (int)t;
+       t -= (scalar_t)floor(t);
+       if(seg >= get_segment_count()) {
+               seg = get_segment_count() - 1;
+               t = 1.0f;
+       }
+
+       Vector3 cp[2];
+       ListNode<Vector3> *iter = const_cast<PolyLine*>(this)->control_points.begin();
+       for(int i=0; i<seg; i++) {
+               iter = iter->next;
+       }
+
+       cp[0] = iter->data;
+       cp[1] = iter->next->data;
+
+       return cp[0] + (cp[1] - cp[0]) * t;
+}
+
+
+bool save_curve(const char *fname, const Curve *curve) {
+       FILE *fp = fopen(fname, "w");
+       if(!fp) {
+               error("failed to save the curve %s", curve->name.c_str());
+               return false;
+       }
+
+       fputs("curve_3dengfx\n", fp);
+       fputs(curve->name.c_str(), fp);
+       fputs("\n", fp);
+
+       if(dynamic_cast<const BSpline*>(curve)) {
+               fputs("bspline\n", fp);
+       } else if(dynamic_cast<const CatmullRomSpline*>(curve)) {
+               fputs("catmullrom\n", fp);
+       } else if(dynamic_cast<const BezierSpline*>(curve)) {
+               fputs("bezier\n", fp);
+       } else if(dynamic_cast<const PolyLine*>(curve)) {
+               fputs("polyline\n", fp);
+       } else {
+               error("unknown spline type, save failed %s", curve->name.c_str());
+               fclose(fp);
+               remove(fname);
+               return false;
+       }
+
+       fprintf(fp, "%d\n", curve->control_points.size());
+
+       const ListNode<Vector3> *node = curve->control_points.begin();
+       while(node) {
+               fprintf(fp, "%f %f %f\n", node->data.x, node->data.y, node->data.z);
+               node = node->next;
+       }
+
+       fclose(fp);
+       return true;
+}
+
+Curve *load_curve(const char *fname) {
+       FILE *fp = fopen(fname, "r");
+       if(!fp) {
+               error("failed to open file %s", fname);
+               return 0;
+       }
+
+       char buffer[256];
+
+       fgets(buffer, 256, fp);
+       if(strcmp(buffer, "curve_3dengfx\n") != 0) {
+               error("load_curve failed, %s is not a curve file", fname);
+               fclose(fp);
+               return 0;
+       }
+
+       Curve *curve;
+
+       fgets(buffer, 256, fp);
+       std::string name = buffer;
+
+       fgets(buffer, 256, fp);
+       if(!strcmp(buffer, "bspline\n")) {
+               curve = new BSpline;
+       } else if(!strcmp(buffer, "catmullrom\n")) {
+               curve = new CatmullRomSpline;
+       } else /*if(!strcmp(buffer, "bezier"))*/ {
+               error("unsupported curve type (%s) or not a curve file", buffer);
+               fclose(fp);
+               return 0;
+       }
+
+       curve->name = name;
+
+       fgets(buffer, 256, fp);
+       if(!isdigit(buffer[0])) {
+               error("load_curve failed, %s is not a valid curve file (count: %s)", fname, buffer);
+               delete curve;
+               fclose(fp);
+               return 0;
+       }
+       int count = atoi(buffer);
+
+       int failed = count;
+       for(int i=0; i<count; i++, failed--) {
+               fgets(buffer, 256, fp);
+               if(!isdigit(buffer[0]) && buffer[0] != '.' && buffer[0] != '-') {
+                       break;
+               }
+               float x = atof(buffer);
+
+               char *ptr = strchr(buffer, ' ');
+               if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
+                       break;
+               }
+               float y = atof(++ptr);
+               
+               ptr = strchr(ptr, ' ');
+               if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
+                       break;
+               }
+               float z = atof(++ptr);
+
+               curve->add_control_point(Vector3(x, y, z));
+       }
+
+       fclose(fp);
+
+       if(failed) {
+               error("load_curve failed to read the data, %s is not a valid curve file", fname);
+               delete curve;
+               return 0;
+       }
+
+       return curve;
+}