# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1450588892 -7200
# Node ID 84a647283237d7c9e13a2f515cc28569bb162bdf
# Parent  099fd7adb900d96ebd33aedc1e4c26152b7447d8
added all the extra functionality to the curve class
first pass at a project-to-curve function (needs more work)

diff -r 099fd7adb900 -r 84a647283237 src/app.cc
--- a/src/app.cc	Sat Dec 19 22:39:18 2015 +0200
+++ b/src/app.cc	Sun Dec 20 07:21:32 2015 +0200
@@ -30,6 +30,8 @@
 static float grid_size = 1.0;
 static SnapMode snap_mode;
 
+static bool show_bounds;
+
 static std::vector<Curve*> curves;
 static Curve *sel_curve;	// selected curve being edited
 static Curve *new_curve;	// new curve being entered
@@ -38,9 +40,7 @@
 
 static Label *weight_label;	// floating label for the cp weight
 
-#ifdef DRAW_MOUSE_POINTER
 static Vector2 mouse_pointer;
-#endif
 
 
 bool app_init(int argc, char **argv)
@@ -151,6 +151,20 @@
 	int numpt = curve->size();
 	int segm = numpt * 16;
 
+	if(show_bounds) {
+		Vector3 bmin, bmax;
+		curve->get_bbox(&bmin, &bmax);
+
+		glLineWidth(1.0);
+		glColor3f(0, 1, 0);
+		glBegin(GL_LINE_LOOP);
+		glVertex2f(bmin.x, bmin.y);
+		glVertex2f(bmax.x, bmin.y);
+		glVertex2f(bmax.x, bmax.y);
+		glVertex2f(bmin.x, bmax.y);
+		glEnd();
+	}
+
 	glLineWidth(curve == hover_curve ? 4.0 : 2.0);
 	if(curve == sel_curve) {
 		glColor3f(0.3, 0.4, 1.0);
@@ -162,7 +176,7 @@
 	glBegin(GL_LINE_STRIP);
 	for(int i=0; i<segm; i++) {
 		float t = (float)i / (float)(segm - 1);
-		Vector2 v = curve->interpolate(t);
+		Vector3 v = curve->interpolate(t);
 		glVertex2f(v.x, v.y);
 	}
 	glEnd();
@@ -183,10 +197,23 @@
 				glColor3f(0.2, 1.0, 0.2);
 			}
 		}
-		Vector2 pt = curve->get_point(i);
+		Vector2 pt = curve->get_point2(i);
 		glVertex2f(pt.x, pt.y);
 	}
 	glEnd();
+
+	// draw the projected mouse point on the selected curve
+	/*
+	if(curve == sel_curve) {
+		Vector3 pp = curve->proj_point(Vector3(mouse_pointer.x, mouse_pointer.y, 0.0));
+
+		glPointSize(5.0);
+		glBegin(GL_POINTS);
+		glColor3f(1, 0.8, 0.2);
+		glVertex2f(pp.x, pp.y);
+		glEnd();
+	}
+	*/
 	glPointSize(1.0);
 }
 
@@ -218,38 +245,29 @@
 			}
 			break;
 
-		case 'l':
-		case 'L':
+		case '1':
+		case '2':
+		case '3':
 			if(sel_curve) {
-				sel_curve->set_type(CURVE_LINEAR);
+				sel_curve->set_type((CurveType)((int)CURVE_LINEAR + key - '1'));
 				post_redisplay();
 			}
 			if(new_curve) {
-				new_curve->set_type(CURVE_LINEAR);
+				new_curve->set_type((CurveType)((int)CURVE_LINEAR + key - '1'));
 				post_redisplay();
 			}
 			break;
 
 		case 'b':
 		case 'B':
-			if(sel_curve) {
-				sel_curve->set_type(CURVE_BSPLINE);
-				post_redisplay();
-			}
-			if(new_curve) {
-				new_curve->set_type(CURVE_BSPLINE);
-				post_redisplay();
-			}
+			show_bounds = !show_bounds;
+			post_redisplay();
 			break;
 
-		case 'h':
-		case 'H':
+		case 'n':
+		case 'N':
 			if(sel_curve) {
-				sel_curve->set_type(CURVE_HERMITE);
-				post_redisplay();
-			}
-			if(new_curve) {
-				new_curve->set_type(CURVE_HERMITE);
+				sel_curve->normalize();
 				post_redisplay();
 			}
 			break;
@@ -263,14 +281,19 @@
 			printf("exported %d curves\n", (int)curves.size());
 			break;
 
-		case 'i':
-		case 'I':
+		case 'l':
+		case 'L':
 			{
 				std::list<Curve*> clist = load_curves("test.curves");
 				if(clist.empty()) {
 					fprintf(stderr, "failed to import curves\n");
 				}
 
+				for(size_t i=0; i<curves.size(); i++) {
+					delete curves[i];
+				}
+				curves.clear();
+
 				int num = 0;
 				std::list<Curve*>::iterator it = clist.begin();
 				while(it != clist.end()) {
@@ -279,6 +302,7 @@
 				}
 				printf("imported %d curves\n", num);
 			}
+			post_redisplay();
 			break;
 		}
 	}
@@ -402,10 +426,8 @@
 	if(!dx && !dy) return;
 
 	Vector2 uv = pixel_to_uv(x, y);
-#ifdef DRAW_MOUSE_POINTER
 	mouse_pointer = uv;
 	post_redisplay();
-#endif
 
 	/* when entering a new curve, have the last (extra) point following
 	 * the mouse until it's entered by a click (see on_click).
diff -r 099fd7adb900 -r 84a647283237 src/curve.cc
--- a/src/curve.cc	Sat Dec 19 22:39:18 2015 +0200
+++ b/src/curve.cc	Sun Dec 20 07:21:32 2015 +0200
@@ -1,10 +1,39 @@
 #include <float.h>
+#include <assert.h>
 #include <algorithm>
 #include "curve.h"
 
 Curve::Curve(CurveType type)
 {
 	this->type = type;
+	bbvalid = true;
+}
+
+Curve::Curve(const Vector4 *cp, int numcp, CurveType type)
+	: Curve(type)
+{
+	this->cp.resize(numcp);
+	for(int i=0; i<numcp; i++) {
+		this->cp[i] = cp[i];
+	}
+}
+
+Curve::Curve(const Vector3 *cp, int numcp, CurveType type)
+	: Curve(type)
+{
+	this->cp.resize(numcp);
+	for(int i=0; i<numcp; i++) {
+		this->cp[i] = Vector4(cp[i].x, cp[i].y, cp[i].z, 1.0f);
+	}
+}
+
+Curve::Curve(const Vector2 *cp, int numcp, CurveType type)
+	: Curve(type)
+{
+	this->cp.resize(numcp);
+	for(int i=0; i<numcp; i++) {
+		this->cp[i] = Vector4(cp[i].x, cp[i].y, 0.0f, 1.0f);
+	}
 }
 
 void Curve::set_type(CurveType type)
@@ -17,9 +46,20 @@
 	return type;
 }
 
+void Curve::add_point(const Vector4 &p)
+{
+	cp.push_back(p);
+	inval_bounds();
+}
+
+void Curve::add_point(const Vector3 &p, float weight)
+{
+	add_point(Vector4(p.x, p.y, p.z, weight));
+}
+
 void Curve::add_point(const Vector2 &p, float weight)
 {
-	cp.push_back(Vector3(p.x, p.y, weight));
+	add_point(Vector4(p.x, p.y, 0.0f, weight));
 }
 
 bool Curve::remove_point(int idx)
@@ -28,23 +68,14 @@
 		return false;
 	}
 	cp.erase(cp.begin() + idx);
+	inval_bounds();
 	return true;
 }
 
-int Curve::nearest_point(const Vector2 &p)
+void Curve::clear()
 {
-	int res = -1;
-	float bestsq = FLT_MAX;
-
-	for(size_t i=0; i<cp.size(); i++) {
-		float d = (get_point(i) - p).length_sq();
-		if(d < bestsq) {
-			bestsq = d;
-			res = i;
-		}
-	}
-
-	return res;
+	cp.clear();
+	inval_bounds();
 }
 
 bool Curve::empty() const
@@ -57,29 +88,45 @@
 	return (int)cp.size();
 }
 
-Vector3 &Curve::operator [](int idx)
+Vector4 &Curve::operator [](int idx)
+{
+	inval_bounds();
+	return cp[idx];
+}
+
+const Vector4 &Curve::operator [](int idx) const
 {
 	return cp[idx];
 }
 
-const Vector3 &Curve::operator [](int idx) const
+const Vector4 &Curve::get_point(int idx) const
 {
 	return cp[idx];
 }
 
-const Vector3 &Curve::get_homo_point(int idx) const
+Vector3 Curve::get_point3(int idx) const
 {
-	return cp[idx];
+	return Vector3(cp[idx].x, cp[idx].y, cp[idx].z);
 }
 
-Vector2 Curve::get_point(int idx) const
+Vector2 Curve::get_point2(int idx) const
 {
 	return Vector2(cp[idx].x, cp[idx].y);
 }
 
 float Curve::get_weight(int idx) const
 {
-	return cp[idx].z;
+	return cp[idx].w;
+}
+
+bool Curve::set_point(int idx, const Vector3 &p, float weight)
+{
+	if(idx < 0 || idx >= (int)cp.size()) {
+		return false;
+	}
+	cp[idx] = Vector4(p.x, p.y, p.z, weight);
+	inval_bounds();
+	return true;
 }
 
 bool Curve::set_point(int idx, const Vector2 &p, float weight)
@@ -87,7 +134,8 @@
 	if(idx < 0 || idx >= (int)cp.size()) {
 		return false;
 	}
-	cp[idx] = Vector3(p.x, p.y, weight);
+	cp[idx] = Vector4(p.x, p.y, 0.0, weight);
+	inval_bounds();
 	return true;
 }
 
@@ -96,7 +144,17 @@
 	if(idx < 0 || idx >= (int)cp.size()) {
 		return false;
 	}
-	cp[idx].z = weight;
+	cp[idx].w = weight;
+	return true;
+}
+
+bool Curve::move_point(int idx, const Vector3 &p)
+{
+	if(idx < 0 || idx >= (int)cp.size()) {
+		return false;
+	}
+	cp[idx] = Vector4(p.x, p.y, p.z, cp[idx].w);
+	inval_bounds();
 	return true;
 }
 
@@ -105,22 +163,207 @@
 	if(idx < 0 || idx >= (int)cp.size()) {
 		return false;
 	}
-	cp[idx] = Vector3(p.x, p.y, cp[idx].z);
+	cp[idx] = Vector4(p.x, p.y, 0.0f, cp[idx].w);
+	inval_bounds();
 	return true;
 }
 
-Vector2 Curve::interpolate(float t, CurveType type) const
+
+int Curve::nearest_point(const Vector3 &p) const
 {
-	if(cp.empty()) {
-		return Vector2(0, 0);
+	int res = -1;
+	float bestsq = FLT_MAX;
+
+	for(size_t i=0; i<cp.size(); i++) {
+		float d = (get_point3(i) - p).length_sq();
+		if(d < bestsq) {
+			bestsq = d;
+			res = i;
+		}
+	}
+	return res;
+}
+
+int Curve::nearest_point(const Vector2 &p) const
+{
+	int res = -1;
+	float bestsq = FLT_MAX;
+
+	for(size_t i=0; i<cp.size(); i++) {
+		float d = (get_point2(i) - p).length_sq();
+		if(d < bestsq) {
+			bestsq = d;
+			res = i;
+		}
+	}
+	return res;
+}
+
+
+void Curve::inval_bounds() const
+{
+	bbvalid = false;
+}
+
+void Curve::calc_bounds() const
+{
+	calc_bbox(&bbmin, &bbmax);
+	bbvalid = true;
+}
+
+void Curve::get_bbox(Vector3 *bbmin, Vector3 *bbmax) const
+{
+	if(!bbvalid) {
+		calc_bounds();
+	}
+	*bbmin = this->bbmin;
+	*bbmax = this->bbmax;
+}
+
+void Curve::calc_bbox(Vector3 *bbmin, Vector3 *bbmax) const
+{
+	if(empty()) {
+		*bbmin = *bbmax = Vector3(0, 0, 0);
+		return;
+	}
+
+	Vector3 bmin = cp[0];
+	Vector3 bmax = bmin;
+	for(size_t i=1; i<cp.size(); i++) {
+		const Vector4 &v = cp[i];
+		for(int j=0; j<3; j++) {
+			if(v[j] < bmin[j]) bmin[j] = v[j];
+			if(v[j] > bmax[j]) bmax[j] = v[j];
+		}
+	}
+	*bbmin = bmin;
+	*bbmax = bmax;
+}
+
+void Curve::normalize()
+{
+	if(!bbvalid) {
+		calc_bounds();
+	}
+
+	Vector3 bsize = bbmax - bbmin;
+	Vector3 boffs = (bbmin + bbmax) * 0.5;
+
+	Vector3 bscale;
+	bscale.x = bsize.x == 0.0f ? 1.0f : 1.0f / bsize.x;
+	bscale.y = bsize.y == 0.0f ? 1.0f : 1.0f / bsize.y;
+	bscale.z = bsize.z == 0.0f ? 1.0f : 1.0f / bsize.z;
+
+	for(size_t i=0; i<cp.size(); i++) {
+		cp[i].x = (cp[i].x - boffs.x) * bscale.x;
+		cp[i].y = (cp[i].y - boffs.y) * bscale.y;
+		cp[i].z = (cp[i].z - boffs.z) * bscale.z;
+	}
+	inval_bounds();
+}
+
+/* Projection to the curve is not correct, but should be good enough for
+ * most purposes.
+ *
+ * First we find the nearest segment (pair of control points), and then
+ * we subdivide between them to find the nearest interpolated point in that
+ * segment.
+ * The incorrect assumption here is that the nearest segment as defined by
+ * the distance of its control points to p, will contain the nearest point
+ * on the curve. This only holds for polylines, and *possibly* bsplines, but
+ * certainly not hermite splines.
+ */
+Vector3 Curve::proj_point(const Vector3 &p) const
+{
+	// TODO fix: select nearest segment based on point-line distance, not distance from the CPs
+	int num_cp = size();
+	if(num_cp <= 0) return p;
+	if(num_cp == 1) return cp[0];
+
+	int idx0 = nearest_point(p);
+	int next_idx = idx0 + 1;
+	int prev_idx = idx0 - 1;
+
+	float next_distsq = next_idx >= num_cp ? FLT_MAX : (get_point3(next_idx) - p).length_sq();
+	float prev_distsq = prev_idx < 0 ? FLT_MAX : (get_point3(prev_idx) - p).length_sq();
+	int idx1 = next_distsq < prev_distsq ? next_idx : prev_idx;
+	assert(idx1 >= 0 && idx1 < num_cp - 1);
+	if(idx0 > idx1) std::swap(idx0, idx1);
+
+	float t0 = 0.0f, t1 = 1.0f;
+	Vector3 pp0 = interpolate_segment(idx0, idx1, 0.0f);
+	Vector3 pp1 = interpolate_segment(idx0, idx1, 1.0f);
+	float dist0 = (pp0 - p).length_sq();
+	float dist1 = (pp1 - p).length_sq();
+	Vector3 pp;
+
+	for(int i=0; i<32; i++) {	// max iterations
+		float t = (t0 + t1) / 2.0;
+		pp = interpolate_segment(idx0, idx1, t);
+		float dist = (pp - p).length_sq();
+
+		// mid point more distant than both control points, nearest cp is closest
+		if(dist > dist0 && dist > dist1) {
+			pp = dist0 < dist1 ? pp0 : pp1;
+			break;
+		}
+
+		if(dist0 < dist1) {
+			t1 = t;
+			dist1 = dist;
+			pp1 = pp;
+		} else {
+			t0 = t;
+			dist0 = dist;
+			pp0 = pp;
+		}
+
+		if(fabs(dist0 - dist1) < 1e-4) {
+			break;
+		}
+	}
+	return pp;
+}
+
+Vector3 Curve::interpolate_segment(int a, int b, float t) const
+{
+	int num_cp = size();
+
+	if(t < 0.0) t = 0.0;
+	if(t > 1.0) t = 1.0;
+
+	Vector4 res;
+	if(type == CURVE_LINEAR || num_cp == 2) {
+		res = lerp(cp[a], cp[b], t);
+	} else {
+		int prev = a <= 0 ? a : a - 1;
+		int next = b >= num_cp - 1 ? b : b + 1;
+
+		if(type == CURVE_HERMITE) {
+			res = catmull_rom_spline(cp[prev], cp[a], cp[b], cp[next], t);
+		} else {
+			res = bspline(cp[prev], cp[a], cp[b], cp[next], t);
+			if(res.w != 0.0f) {
+				res.x /= res.w;
+				res.y /= res.w;
+			}
+		}
+	}
+
+	return Vector3(res.x, res.y, res.z);
+}
+
+Vector3 Curve::interpolate(float t) const
+{
+	if(empty()) {
+		return Vector3(0, 0, 0);
 	}
 
 	int num_cp = (int)cp.size();
 	if(num_cp == 1) {
-		return Vector2(cp[0].x, cp[0].y);
+		return Vector3(cp[0].x, cp[0].y, cp[0].z);
 	}
 
-	Vector3 res;
 	int idx0 = std::min((int)floor(t * (num_cp - 1)), num_cp - 2);
 	int idx1 = idx0 + 1;
 
@@ -129,35 +372,17 @@
 	float t1 = (float)idx1 * dt;
 
 	t = (t - t0) / (t1 - t0);
-	if(t < 0.0) t = 0.0;
-	if(t > 1.0) t = 1.0;
 
-	if(type == CURVE_LINEAR || num_cp <= 2) {
-		res = lerp(cp[idx0], cp[idx1], t);
-	} else {
-		int idx_prev = idx0 <= 0 ? idx0 : idx0 - 1;
-		int idx_next = idx1 >= num_cp - 1 ? idx1 : idx1 + 1;
+	return interpolate_segment(idx0, idx1, t);
+}
 
-		if(type == CURVE_HERMITE) {
-			res = catmull_rom_spline(cp[idx_prev], cp[idx0], cp[idx1], cp[idx_next], t);
-		} else {
-			res = bspline(cp[idx_prev], cp[idx0], cp[idx1], cp[idx_next], t);
-			if(res.z != 0.0f) {
-				res.x /= res.z;
-				res.y /= res.z;
-			}
-		}
-	}
-
+Vector2 Curve::interpolate2(float t) const
+{
+	Vector3 res = interpolate(t);
 	return Vector2(res.x, res.y);
 }
 
-Vector2 Curve::interpolate(float t) const
-{
-	return interpolate(t, type);
-}
-
-Vector2 Curve::operator ()(float t) const
+Vector3 Curve::operator ()(float t) const
 {
 	return interpolate(t);
 }
diff -r 099fd7adb900 -r 84a647283237 src/curve.h
--- a/src/curve.h	Sat Dec 19 22:39:18 2015 +0200
+++ b/src/curve.h	Sun Dec 20 07:21:32 2015 +0200
@@ -12,37 +12,72 @@
 
 class Curve {
 private:
-	std::vector<Vector3> cp;
+	std::vector<Vector4> cp;
 	CurveType type;
 
+	// bounding box
+	mutable Vector3 bbmin, bbmax;
+	mutable bool bbvalid;
+
+	void calc_bounds() const;
+	void inval_bounds() const;
+
 public:
 	Curve(CurveType type = CURVE_HERMITE);
+	Curve(const Vector4 *cp, int numcp, CurveType type = CURVE_HERMITE); // homogenous
+	Curve(const Vector3 *cp, int numcp, CurveType type = CURVE_HERMITE); // 3D points, w=1
+	Curve(const Vector2 *cp, int numcp, CurveType type = CURVE_HERMITE); // 2D points, z=0, w=1
 
 	void set_type(CurveType type);
 	CurveType get_type() const;
 
+	void add_point(const Vector4 &p);
+	void add_point(const Vector3 &p, float weight = 1.0f);
 	void add_point(const Vector2 &p, float weight = 1.0f);
 	bool remove_point(int idx);
 
-	int nearest_point(const Vector2 &p);
+	void clear();		// remove all control points
+	bool empty() const;	// true if 0 control points
+	int size() const;	// returns number of control points
+	// access operators for control points
+	Vector4 &operator [](int idx);
+	const Vector4 &operator [](int idx) const;
+	const Vector4 &get_point(int idx) const;
 
-	bool empty() const;
-	int size() const;
-	Vector3 &operator [](int idx);
-	const Vector3 &operator [](int idx) const;
-
-	const Vector3 &get_homo_point(int idx) const;	// homogeneous point
-	Vector2 get_point(int idx) const;
+	Vector3 get_point3(int idx) const;
+	Vector2 get_point2(int idx) const;
 	float get_weight(int idx) const;
 
+	bool set_point(int idx, const Vector3 &p, float weight = 1.0f);
 	bool set_point(int idx, const Vector2 &p, float weight = 1.0f);
 	bool set_weight(int idx, float weight);
 	// move point without changing its weight
+	bool move_point(int idx, const Vector3 &p);
 	bool move_point(int idx, const Vector2 &p);
+ 
+	int nearest_point(const Vector3 &p) const;
+	// nearest control point on the 2D plane z=0
+	int nearest_point(const Vector2 &p) const;
 
-	Vector2 interpolate(float t, CurveType type) const;
-	Vector2 interpolate(float t) const;
-	Vector2 operator ()(float t) const;
+	/* get_bbox returns the axis-aligned bounding box of the curve's
+	 * control points.
+	 * NOTE: hermite curves can go outside of the bounding box of their control points
+	 * NOTE: lazy calculation of bounds is performed, use calc_bbox in multithreaded programs
+	 */
+	void get_bbox(Vector3 *bbmin, Vector3 *bbmax) const;
+	void calc_bbox(Vector3 *bbmin, Vector3 *bbmax) const;
+	// normalize the curve's bounds to coincide with the unit cube
+	void normalize();
+
+	// project a point to the curve (nearest point on the curve)
+	Vector3 proj_point(const Vector3 &p) const;
+	// equivalent to (proj_point(p) - p).length();
+	float distance(const Vector3 &p) const;
+
+	Vector3 interpolate_segment(int a, int b, float t) const;
+	Vector3 interpolate(float t) const;
+	Vector2 interpolate2(float t) const;
+	Vector3 operator ()(float t) const;
 };
 
 #endif	// CURVE_H_
diff -r 099fd7adb900 -r 84a647283237 src/curvefile.cc
--- a/src/curvefile.cc	Sat Dec 19 22:39:18 2015 +0200
+++ b/src/curvefile.cc	Sun Dec 20 07:21:32 2015 +0200
@@ -46,8 +46,8 @@
 	fprintf(fp, "    type %s\n", curve_type_str(curve->get_type()));
 	fprintf(fp, "    cpcount %d\n", curve->size());
 	for(int i=0; i<curve->size(); i++) {
-		Vector3 cp = curve->get_homo_point(i);
-		fprintf(fp, "    cp %g %g %g\n", cp.x, cp.y, cp.z);
+		Vector4 cp = curve->get_point(i);
+		fprintf(fp, "    cp %g %g %g %g\n", cp.x, cp.y, cp.z, cp.w);
 	}
 	fprintf(fp, "}\n");
 	return true;
@@ -83,8 +83,6 @@
 		}
 		s.push_back(c);
 	}
-
-	printf("TOKEN: \"%s\"\n", s.c_str());
 	return s;
 }
 
@@ -159,13 +157,13 @@
 			if(tok != "cp") {
 				goto err;
 			}
-			Vector3 cp;
-			for(int i=0; i<3; i++) {
+			Vector4 cp;
+			for(int i=0; i<4; i++) {
 				if(!expect_float(fp, &cp[i])) {
 					goto err;
 				}
 			}
-			curve->add_point(Vector2(cp.x, cp.y), cp.z);
+			curve->add_point(Vector3(cp.x, cp.y, cp.z), cp.w);
 		}
 	}