ptrude/: many major math fixes, especially in stretch_path

ptrude/Makefile

@@ -44,6 +44,9 @@ $(MAIN):		$(OBJS)
 try:		$(MAIN)
 		./$(MAIN) -d try 1 0.1 | tee out
 
+t2:		$(MAIN)
+		./$(MAIN) tp ts 25 0.1 >out.stl
+
 clean:
 		rm -f $(OBJS) $(OBJS:.o=.d)
 

ptrude/extrude.c

@@ -14,14 +14,14 @@
  * Known bugs:
  *
  * - negative x coordinates in shape produce overlapping faces
- * - the rounding is only accurate for x = 0. For x > 0, we need to add
- *   more points. Tricky.
  */
 
 
 #include <stdlib.h>
 #include <math.h>
+#include <assert.h>
 
+#include "ptrude.h"
 #include "path.h"
 #include "extrude.h"
 
@@ -34,6 +34,26 @@ static void cvt_3d(struct point *p, const struct vertex *v, double z)
 }
 
 
+static void size(const struct path *p)
+{
+	const struct vertex *v, *next;
+	double s = 0, sl = 0;
+	int n = 0;
+
+	if (!debug)
+		return;
+	for (v = p->vertices; v; v = next) {
+		next = v->next;
+		if (next)
+			s += hypot(v->x-next->x, v->y-next->y);
+		sl += v->len;
+		n++;
+	}
+	fprintf(stderr, "%d virt %g real %g\n", n, sl, s);
+	
+}
+
+
 static void mesh(const struct path *pa, double za,
     const struct path *pb, double zb,
     void (*face)(void *data, struct point a, struct point b, struct point c),
@@ -47,6 +67,7 @@ static void mesh(const struct path *pa, double za,
 
 	cvt_3d(&a3, a, za);
 	cvt_3d(&b3, b, zb);
+
 	while (a->next || b->next) {
 		da = fabs(sa+a->len-sb);
 		db = fabs(sb+b->len-sa);
@@ -83,14 +104,18 @@ void extrude(const struct path *path, const struct path *shape,
 	if (!v || !v->next)
 		return;
 
-	tmp = stretch_path(path, v->x);
+	tmp = stretch_path(path, v->x, r);
+	size(tmp);
 	prev = round_path(tmp, r, d);
 	free_path(tmp);
+	size(prev);
 
 	while (v->next) {
 		next = v->next;
-		tmp = stretch_path(path, next->x);
+		tmp = stretch_path(path, next->x, r);
+		size(tmp);
 		curr = round_path(tmp, r, d);
+		size(curr);
 		free_path(tmp);
 
 		mesh(prev, v->y, curr, next->y, face, data);

ptrude/path.c

@@ -14,6 +14,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <math.h>
+#include <assert.h>
 
 #include "ptrude.h"
 #include "path.h"
@@ -127,14 +128,16 @@ static void adjust_length(struct vertex *from, struct vertex *to, double len)
 	struct vertex *v;
 	double sum, f;
 
+	if (from == to)
+		return;
 	sum = 0;
-	for (v = from; v != to; v = v->next) {
+	for (v = from->next; v != to; v = v->next) {
 		v->len = hypot(v->x-v->next->x, v->y-v->next->y);
 		sum += v->len;
 	}
 
 	f = len/sum;
-	for (v = from; v != to; v = v->next)
+	for (v = from->next; v != to; v = v->next)
 		v->len *= f;
 }
 
@@ -178,6 +181,7 @@ static struct vertex *corner(struct path *path, struct vertex *a,
 	int n;		/* number of connecting segments (0 if none) */
 
 	double f;	/* scale factor; various uses */
+	double fa, fb;	/* scale factors for first and last vertex */
 	double ang;	/* current angle, for iteration */
 	double x, y;	/* current position; for iteration */
 	int i;		/* segment; for iteration */
@@ -226,15 +230,21 @@ static struct vertex *corner(struct path *path, struct vertex *a,
 	 * available length in half, one for the inbound arc, the other for the
 	 * outbound arc.
 	 */
+
+	/*
+	 * @@@ Our error checking is a bit overzealous and doesn't provide
+	 * enough information to debug any problems. Turn errors into warnings
+	 * for now.
+	 */
 	if (aa/2 < s) {
 		fprintf(stderr, "first vector is too short (%g/2 < %g)\n",
 		    aa, s);
-		exit(1);
+//		exit(1);
 	}
 	if (bb/2 < s) {
 		fprintf(stderr, "second vector is too short (%g/2 < %g)\n",
 		    bb, s);
-		exit(1);
+//		exit(1);
 	}
 
 	/*
@@ -279,11 +289,11 @@ static struct vertex *corner(struct path *path, struct vertex *a,
 	 * Step 4: emit the starting point of the arc
 	 */
 
-	f = s/aa;
+	fa = s/aa;
-	x = b->x-f*ax;
+	x = b->x-fa*ax;
-	y = b->y-f*ay;
+	y = b->y-fa*ay;
 	v0 = add_vertex(path, x, y, b->r, b->d, b->tag);
-	a->len = a->len-s;
+	v0->len = a->len*(1-fa);
 
 	/*
 	 * Step 5: determine if we need intermediate points. If yes, how many,
@@ -337,15 +347,15 @@ static struct vertex *corner(struct path *path, struct vertex *a,
 	 * Step 6: emit the finishing point of the arc
 	 */
 
-	f = s/bb;
+	fb = s/bb;
-	v1 = add_vertex(path, b->x+f*bx, b->y+f*by, 0, 0, NULL);
+	v1 = add_vertex(path, b->x+fb*bx, b->y+fb*by, 0, 0, NULL);
-	v1->len = b->len-s;
+	v1->len = b->len*(1-fb);
 
 	/*
 	 * Step 7: adjust the nominal length of the segments
 	 */
 
-	adjust_length(v0, v1, 2*s);
+	adjust_length(v0, v1, a->len*fa+b->len*fb);
 
 
 	return v1;
@@ -386,18 +396,19 @@ struct path *round_path(const struct path *path, double r, double d)
 
 
 static void move_vertex(struct path *path, const struct vertex *v,
-    double nx, double ny, double d)
+    double nx, double ny, double dist, double r)
 {
 	struct vertex *new;
 
 	new = clone_vertex(v);
-	new->x += nx*d;
+	new->x += nx*dist;
-	new->y += ny*d;
+	new->y += ny*dist;
+	new->r = r;
 	append_vertex(path, new);
 }
 
 
-struct path *stretch_path(const struct path *path, double d)
+struct path *stretch_path(const struct path *path, double dist, double r)
 {
 	struct path *new;		/* new path */
 	const struct vertex *v;		/* current vertex (for iteration) */
@@ -413,9 +424,13 @@ struct path *stretch_path(const struct path *path, double d)
 	nx = b->y-a->y;
 	ny = a->x-b->x;
 	f = hypot(nx, ny);
-	move_vertex(new, a, nx/f, ny/f, d);
+	if (a->r)
+		r = a->r;
+	move_vertex(new, a, nx/f, ny/f, dist, r);
 
 	for (v = path->vertices->next; v->next; v = v->next) {
+		double tmp;
+
 		b = v;
 		c = v->next;
 
@@ -425,12 +440,40 @@ struct path *stretch_path(const struct path *path, double d)
 		nx = tx/f;
 		ny = ty/f;
 
+		tmp = f;
+
 		tx = c->y-b->y;
 		ty = b->x-c->x;
 		f = hypot(tx, ty);
 		nx += tx/f;
 		ny += ty/f;
-		move_vertex(new, b, nx/2, ny/2, d);
+		if (b->r)
+			r = b->r;
+
+		f = hypot(nx, ny);
+		nx /= f;
+		ny /= f;
+
+		/*
+		 * We have this far:
+		 * nx, ny = normal on corner, normalized
+		 * tmp = |a|, length of vector "a" (A -> B)
+		 * dist = the distance by which we stretch
+		 *
+		 * As shown in stretch.fig, we the length we need is
+		 * d' = d/cos(90-t)
+		 *
+		 * With
+		 * http://en.wikipedia.org/wiki/Trigonometric_identities#Symmetry
+		 * cos(90-t) = sin t = (n x a)/(|n|*|a|)
+		 *
+		 * Thus
+		 * d' = d/sin(t) - d*(|n|*|a|)/(n x a)
+		 *    = d/sin(t) - d*|a|/(n x a)
+		 */
+		tmp = dist*tmp/(nx*(b->y-a->y)-ny*(b->x-a->x));
+
+		move_vertex(new, b, nx, ny, tmp, r+dist);
 
 		a = v;
 	}
@@ -438,7 +481,9 @@ struct path *stretch_path(const struct path *path, double d)
 	nx = v->y-a->y;
 	ny = a->x-v->x;
 	f = hypot(nx, ny);
-	move_vertex(new, v, nx/f, ny/f, d);
+	if (v->r)
+		r = v->r;
+	move_vertex(new, v, nx/f, ny/f, dist, r);
 
 	return new;
 }
@@ -485,6 +530,7 @@ struct path *load_path(FILE *file)
 		tag = NULL;
 	}
 
+	path_set_length(path);
 	return path;
 }
 

ptrude/path.h

@@ -42,7 +42,7 @@ void free_path(struct path *path);
 
 double path_set_length(struct path *path);
 struct path *round_path(const struct path *path, double r, double d);
-struct path *stretch_path(const struct path *path, double d);
+struct path *stretch_path(const struct path *path, double dist, double r);
 
 struct path *load_path(FILE *file);
 void save_path(FILE *file, const struct path *path);

ptrude/stretch.fig

@@ -0,0 +1,62 @@
+#FIG 3.2  Produced by xfig version 3.2.5b
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+5 1 0 1 0 7 50 -1 -1 0.000 0 1 0 0 4275.476 4948.766 4725 4770 4507 4524 4275 4465
+5 1 0 1 0 7 50 -1 -1 0.000 0 1 0 0 3184.134 5404.540 4275 5400 4160 4917 3912 4592
+5 1 0 1 0 7 50 -1 -1 0.000 0 1 0 0 3189.289 5402.777 8865 5400 8394 3139 7314 1504
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 2.00 90.00 120.00
+	 4275 7875 4275 4950
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 2.00 90.00 120.00
+	 4275 4950 2250 2925
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 1.00 90.00 120.00
+	 5647 4927 6592 3982
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 2.00 90.00 120.00
+	 4275 4950 6570 4005
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 7 1 0 2
+	0 0 1.00 90.00 120.00
+	 4275 4950 5625 4950
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 7 1 0 2
+	0 0 1.00 90.00 120.00
+	 4275 4950 5220 4005
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 7 1 0 2
+	0 0 1.00 90.00 120.00
+	 5220 4005 6570 4005
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 1.00 90.00 120.00
+	 4275 4950 4275 2475
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 1 2
+	0 0 1.00 90.00 120.00
+	0 0 1.00 90.00 120.00
+	 4320 5850 8820 5850
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 0 0 3
+	 8865 6705 8865 3060 7245 1440
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 1 2
+	0 0 1.00 90.00 120.00
+	0 0 1.00 90.00 120.00
+	 4365 4995 8775 3150
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 3195 5310 3195 5490
+2 1 0 2 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 3105 5400 3285 5400
+2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 3195 5400 4275 5400
+4 1 0 50 -1 18 18 0.0000 4 150 165 4095 6570 a\001
+4 1 0 50 -1 18 18 0.0000 4 210 180 3465 3825 b\001
+4 1 0 50 -1 18 18 0.0000 4 150 330 4995 5220 na\001
+4 1 0 50 -1 18 18 0.0000 4 150 165 5490 4365 n\001
+4 1 0 50 -1 18 18 0.0000 4 210 345 4770 4230 nb\001
+4 1 0 50 -1 18 18 0.0000 4 195 90 4410 4770 t\001
+4 1 0 50 -1 18 18 0.0000 4 270 1455 7470 3330 d/cos(90-t)\001
+4 1 0 50 -1 18 18 0.0000 4 210 180 6525 5805 d\001
+4 0 0 50 -1 18 18 0.0000 4 270 4845 4950 7875 cos(90-t) = sin t = n x a (normalized)\001
+4 1 0 50 -1 18 18 0.0000 4 150 105 3735 5310 r\001

ptrude/tp

@@ -7,4 +7,4 @@
 0 300
 -100 200
 -100 100
-0 100
+50 100