Instead of performing the tranformations stepwise for each point, pre-calculate

them once.

- solidify/matrix.h, solidify/matrix.c: 2D matrix operations
- solidify/Makefile (OBJS): added matrix.o
- solidify/face.h (struct matrix): moved to solidify/matrix.h
- solidify/overlap.c (point, merge_matrix, draw_map): precalculate a single
  vA+b transformation instead of stepwise transforming the vector inside
  the loop
- solidify/overlap.c (do_shift): now that the math is right (or at least
  better), reverse the shift

solidify/Makefile

@@ -12,8 +12,8 @@
 
 SHELL = /bin/bash
 
-OBJS = array.o face.o histo.o level.o overlap.o solid.o solidify.o style.o \
-       util.o
+OBJS = array.o face.o histo.o level.o matrix.o overlap.o solid.o solidify.o \
+       style.o util.o
 
 CFLAGS_WARN = -Wall -Wshadow -Wmissing-prototypes \
 	      -Wmissing-declarations -Wno-format-zero-length

solidify/face.h

@@ -14,21 +14,9 @@
 #define FACE_H
 
 #include "array.h"
+#include "matrix.h"
 
 
-/*
- * 2D transformation:
- *
- * x' = x*a[0][0]+y*a[0][1]+b[0]
- * y' = x*a[1][0]+y*a[1][1]+b[1]
- */
-
-
-struct matrix {
-	double a[2][2];
-	double b[2];
-};
-
 struct face {
 	struct array *a;
 	int sx, sy;	/* size */
@@ -47,4 +35,4 @@ static inline double face_z0(const struct face *f, int x, int y)
 
 struct face *read_face(const char *name);
 
-#endif /* FACE_H */
+#endif /* !FACE_H */

solidify/matrix.c

@@ -0,0 +1,66 @@
+/*
+ * matrix.c - 2D matrix operations
+ *
+ * Written 2010 by Werner Almesberger
+ * Copyright 2010 by Werner Almesberger
+ *
+ * This program 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.
+ */
+
+
+#include <assert.h>
+
+#include "matrix.h"
+
+
+void matrix_identity(struct matrix *m)
+{
+	m->a[0][0] = m->a[1][1] = 1;
+	m->a[0][1] = m->a[1][0] = 0;
+	m->b[0] = m->b[1] = 0;
+}
+
+
+void matrix_invert(const double m[2][2], double res[2][2])
+{
+	double det = m[0][0]*m[1][1]-m[0][1]*m[1][0];
+
+	assert(res != (void *) m);
+	res[0][0] = m[1][1]/det;
+	res[0][1] = -m[0][1]/det;
+	res[1][0] = -m[1][0]/det;
+	res[1][1] = m[0][0]/det;
+}
+
+
+void matrix_mult(double a[2][2], double b[2][2], double res[2][2])
+{
+	assert(res != a);
+	assert(res != b);
+	res[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0];
+	res[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1];
+	res[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0];
+	res[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1];
+}
+
+
+void matrix_multv(const double v[2], double m[2][2], double res[2])
+{
+	double tmp;
+
+	tmp = v[0]*m[0][0]+v[1]*m[0][1];
+	res[1] = v[0]*m[1][0]+v[1]*m[1][1];
+	res[0] = tmp;
+}
+
+
+void matrix_copy(double from[2][2], double to[2][2])
+{
+	to[0][0] = from[0][0];
+	to[0][1] = from[0][1];
+	to[1][0] = from[1][0];
+	to[1][1] = from[1][1];
+}

solidify/matrix.h

@@ -0,0 +1,45 @@
+/*
+ * matrix.h - 2D matrix operations
+ *
+ * Written 2010 by Werner Almesberger
+ * Copyright 2010 by Werner Almesberger
+ *
+ * This program 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.
+ */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+/*
+ * 2D transformation:
+ *
+ * x' = x*a[0][0]+y*a[0][1]+b[0]
+ * y' = x*a[1][0]+y*a[1][1]+b[1]
+ */
+
+
+struct matrix {
+	double a[2][2];
+	double b[2];
+};
+
+
+void matrix_identity(struct matrix *m);
+void matrix_invert(const double m[2][2], double res[2][2]);
+void matrix_mult(double a[2][2], double b[2][2], double res[2][2]);
+void matrix_multv(const double v[2], double m[2][2], double res[2]);
+void matrix_copy(double from[2][2], double to[2][2]);
+
+
+static inline void matrix_map(int x, int y, const struct matrix *m,
+    double *res_x, double *res_y)
+{
+	*res_x = x*m->a[0][0]+y*m->a[0][1]+m->b[0];
+	*res_y = x*m->a[1][0]+y*m->a[1][1]+m->b[1];
+
+}
+
+#endif /* !MATRIX_H */

solidify/overlap.c

@@ -92,22 +92,16 @@ static double zmix(struct face *f, double x, double y)
  * - our model runs from min_x to max_x. Its center is at cx.
  */
 
-static void point(const struct solid *s, int x, int y, guchar *p)
+static void point(const struct solid *s, int x, int y, guchar *p,
+    const struct matrix *ma, const struct matrix *mb)
 {
 	double za, zb, z;
-	int xa, xb, ya, yb;
 	double xaf, xbf, yaf, ybf;
 
-	xa = x-sx(s)/2;
-	ya = y-sy(s)/2;
-	xaf = xa*s->a->m.a[0][0]+ya*s->a->m.a[0][1]+s->a->m.b[0]+s->a->cx;
-	yaf = xa*s->a->m.a[1][0]+ya*s->a->m.a[1][1]+s->a->m.b[1]+s->a->cy;
-	za = zmix(s->a, xaf, yaf);
+	matrix_map(x, y, ma, &xaf, &yaf);
+	matrix_map(x, y, mb, &xbf, &ybf);
 
-	xb = x-sx(s)/2;
-	yb = (sy(s)-1)/2-y;
-	xbf = xb*s->b->m.a[0][0]+yb*s->b->m.a[0][1]+s->b->m.b[0]+s->b->cx;
-	ybf = xb*s->b->m.a[1][0]+yb*s->b->m.a[1][1]+s->b->m.b[1]+s->b->cy;
+	za = zmix(s->a, xaf, yaf);
 	zb = zmix(s->b, xbf, ybf);
 
 	if (za == UNDEF_F && zb == UNDEF_F)
@@ -158,19 +152,115 @@ static void point(const struct solid *s, int x, int y, guchar *p)
 }
 
 
+static void merge_matrix(struct matrix *m, const struct solid *s,
+    const struct face *f)
+{
+	double tm[2][2], tm2[2][2];
+	double tv[2];
+	double f_x, f_y;
+
+	/*
+	 * Finally, we convert to model matrix coordinates.
+	 *
+	 * v' = v+c
+	 */
+
+	m->b[0] += f->cx;
+	m->b[1] += f->cy;
+
+	/*
+	 * Apply shrinkage caused by rotation out of z0. We use that
+	 * cos a = sqrt(1-sin^2 a)
+	 */
+
+	f_x = 1.0/sqrt(1-f->fx*f->fx);
+	f_y = 1.0/sqrt(1-f->fy*f->fy);
+
+	m->a[0][0] *= f_x;
+	m->a[0][1] *= f_x;
+	m->b[0] *= f_x;
+	m->a[1][0] *= f_y;
+	m->a[1][1] *= f_y;
+	m->b[1] *= f_y;
+
+	/*
+	 * The transformation matrix f->m describes a transformation of
+	 * (centered) model coordinates. We therefore have to reverse it:
+	 *
+	 * v = v'A+b
+	 * v-b = v'A
+	 * (v-b)A^-1 = v'
+	 * vA^-1-bA^-1 = v'
+	 */
+
+	matrix_invert(f->m.a, tm);
+	matrix_multv(f->m.b, tm, tv);
+	tv[0] = -tv[0];
+	tv[1] = -tv[1];
+
+	/*
+	 * Merge with the transformations we have so far:
+	 *
+	 * v' = vA1+b1		the transformation we have so far
+	 * v'' = v'A2+b2	the transformation we apply
+	 *
+	 * v'' = (vA1+b1)A2+b2
+	 * v'' = vA1A2+b1A2+b2
+	 */
+
+	/*
+	 * So far, the theory. To make it really work, we have to calculate
+	 * v'' = vA1A2+b1+b2
+	 * duh ?!?
+	 */
+
+	matrix_mult(m->a, tm, tm2);	/* A1A2 */
+	matrix_copy(tm2, m->a);
+//	matrix_multv(m->b, tm, m->b);	/* b1A2 */
+	m->b[0] += tv[0];		/* b2 */
+	m->b[1] += tv[1];
+
+	/*
+	 * Our input is a screen coordinate, its origin is in a corner so we
+	 * first have to make it center-based:
+	 *
+	 * v' = (v-s/2)A+b
+	 * v' = vA+(b-s/2*A)
+	 */
+
+	tv[0] = sx(s)/2;
+	tv[1] = sy(s)/2;
+	matrix_multv(tv, m->a, tv);
+	m->b[0] -= tv[0];
+	m->b[1] -= tv[1];
+}
+
+
 static void draw_map(GtkWidget *widget, struct solid *s)
 {
 	guchar *rgbbuf, *p;
 	int x, y;
+	struct matrix ma = {
+		.a = { { 1, 0 }, { 0, 1 } },
+		.b = { 0, 0 },
+	};
+	struct matrix mb = {
+		.a = { { 1, 0 }, { 0, 1 } },	/* @@@ why not a[1][1] = -1 ? */
+		.b = { 0, 0 },
+	};
 
 	rgbbuf = p = calloc(sx(s)*sy(s), 3);
 	if (!rgbbuf) {
 		perror("calloc");
 		exit(1);
 	}
+
+	merge_matrix(&ma, s, s->a);
+	merge_matrix(&mb, s, s->b);
+
 	for (y = sy(s)-1; y >= 0; y--)
 		for (x = 0; x != sx(s) ; x++) {
-			point(s, x, y, p);
+			point(s, x, y, p, &ma, &mb);
 			p += 3;
 		}
 	gdk_draw_rgb_image(widget->window,
@@ -221,8 +311,8 @@ static void rotate(struct matrix *m, double r)
 
 static void do_shift(struct matrix *m, int dx, int dy)
 {
-	m->b[0] -= dx;
-	m->b[1] += dy;
+	m->b[0] += dx;
+	m->b[1] -= dy;
 }