Use rotation matrix

Reuse affine transformation code from variations.

4 files changed, 60 insertions, 26 deletions

diff --git a/math.h b/math.h
index defc7c6..96b9d5b 100644
--- a/math.h
+++ b/math.h
@@ -30,6 +30,55 @@
 
 #define clamp(a,min,max) (a > max ? max : (a < min ? min : a))
 
+/*	Apply affine coordinate transformation
+ */
+inline double2 apply_affine (const double2 in, const double2 matrix[3]) {
+	return matrix[0] * in[0] + matrix[1] * in[1] + matrix[2];
+}
+
+/*	Create affine rotation matrix, angle in degree
+ */
+inline void rotate (const double angle, double2 matrix[3]) {
+	double s, c;
+	sincos (angle * 2.0 * M_PI / 360.0, &s, &c);
+	matrix[0] = (double2) { c, s };
+	matrix[1] = (double2) { -s, c };
+	matrix[2] = (double2) { 0.0, 0.0 };
+}
+
+/*	Create affine translation matrix
+ */
+inline void translate (const double2 xy, double2 matrix[3]) {
+	matrix[0] = (double2) { 1.0, 0.0 };
+	matrix[1] = (double2) { 0.0, 1.0 };
+	matrix[2] = xy;
+}
+
+/*	Multiply two affine matrices a, b and store the result in c.
+ *
+ *	The last row of each matrix is assumed to be 0, 0, 1.
+ */
+inline void matrixmul (const double2 a[3], const double2 b[3], double2 c[3]) {
+	c[0] = a[0] * b[0][0] + a[1] * b[0][1];
+	c[1] = a[0] * b[1][0] + a[1] * b[1][1];
+	c[2] = a[0] * b[2][0] + a[1] * b[2][1] + a[2];
+}
+
+/* Create rotation around center. Note that matrix multiplication is
+ * right-associative, thus A*B*C == A*(B*C) */
+inline void rotate_center (const double2 center, const double angle, double2 out[3]) {
+	double2 rot[3], trans_a[3], trans_b[3], tmp[3];
+	translate (-1.0 * center, trans_a);
+	rotate (angle, rot);
+	translate (center, trans_b);
+	matrixmul (rot, trans_a, tmp);
+	matrixmul (trans_b, tmp, out);
+}
+
+inline double sum(const double2 in) {
+	return in[0] + in[1];
+}
+
 /*	Vector wrapping function, could be replaced by true vector functions later
  */
 
diff --git a/private.h b/private.h
index 1f27600..4cffb2d 100644
--- a/private.h
+++ b/private.h
@@ -48,7 +48,7 @@ typedef struct {
    unsigned short *xform_distrib;    /* Distribution of xforms based on weights */
    flam3_frame *spec; /* Frame contains timing information */
    double bounds[4]; /* Corner coords of viewable area */
-   double rot[2][2]; /* Rotation transformation */
+   double2 rot[3]; /* Rotation transformation */
    double size[2];
    int width, height; /* buffer width/height */
    double ws0, wb0s0, hs1, hb1s1; /* shortcuts for indexing */
diff --git a/rect.c b/rect.c
index d0b767e..0cea6ec 100644
--- a/rect.c
+++ b/rect.c
@@ -219,20 +219,17 @@ static void iter_thread(void *fth) {
 
       /* Put them in the bucket accumulator */
       for (j = 0; j < sub_batch_size; j++) {
-         double p0, p1, p00, p11;
-         const double4 p = fthp->iter_storage[j];
+         double4 p = fthp->iter_storage[j];
 
          if (fthp->cp.rotate != 0.0) {
-            p00 = p[0] - fthp->cp.rot_center[0];
-            p11 = p[1] - fthp->cp.rot_center[1];
-            p0 = p00 * ficp->rot[0][0] + p11 * ficp->rot[0][1] + fthp->cp.rot_center[0];
-            p1 = p00 * ficp->rot[1][0] + p11 * ficp->rot[1][1] + fthp->cp.rot_center[1];
-         } else {
-            p0 = p[0];
-            p1 = p[1];
+		 	const double2 p01 = (double2) { p[0], p[1] };
+		 	const double2 rotatedp = apply_affine (p01, ficp->rot);
+		 	p[0] = rotatedp[0];
+		 	p[1] = rotatedp[1];
          }
 
-         if (p0 >= ficp->bounds[0] && p1 >= ficp->bounds[1] && p0 <= ficp->bounds[2] && p1 <= ficp->bounds[3]) {
+         if (p[0] >= ficp->bounds[0] && p[1] >= ficp->bounds[1] &&
+		     p[0] <= ficp->bounds[2] && p[1] <= ficp->bounds[3]) {
 
             double logvis=1.0;
 
@@ -249,7 +246,7 @@ static void iter_thread(void *fth) {
 			   interpcolor *= logvis;
             }
 
-            ficp->buckets[(int)(ficp->ws0 * p0 - ficp->wb0s0) + ficp->width * (int)(ficp->hs1 * p1 - ficp->hb1s1)] += interpcolor;
+            ficp->buckets[(int)(ficp->ws0 * p[0] - ficp->wb0s0) + ficp->width * (int)(ficp->hs1 * p[1] - ficp->hb1s1)] += interpcolor;
 
          }
       }
@@ -514,10 +511,8 @@ int render_rectangle(flam3_frame *spec, void *out,
             fic.bounds[3] = corner1 + image_height / ppuy + t1 + shift;
             fic.size[0] = 1.0 / (fic.bounds[2] - fic.bounds[0]);
             fic.size[1] = 1.0 / (fic.bounds[3] - fic.bounds[1]);
-            fic.rot[0][0] = cos(cp.rotate * 2 * M_PI / 360.0);
-            fic.rot[0][1] = -sin(cp.rotate * 2 * M_PI / 360.0);
-            fic.rot[1][0] = -fic.rot[0][1];
-            fic.rot[1][1] = fic.rot[0][0];
+			rotate_center ((double2) { cp.rot_center[0], cp.rot_center[1] },
+						   cp.rotate, fic.rot);
             fic.ws0 = fic.width * fic.size[0];
             fic.wb0s0 = fic.ws0 * fic.bounds[0];
             fic.hs1 = fic.height * fic.size[1];
diff --git a/variations.c b/variations.c
index 44e4eb6..41445d5 100644
--- a/variations.c
+++ b/variations.c
@@ -2025,16 +2025,6 @@ int prepare_precalc_flags(flam3_genome *cp) {
    return(0);
 }
 
-/*	Apply affine coordinate transformation
- */
-static double2 apply_affine (const double2 in, const double2 matrix[3]) {
-	return matrix[0] * in[0] + matrix[1] * in[1] + matrix[2];
-}
-
-static double sum(const double2 in) {
-	return in[0] + in[1];
-}
-
 int apply_xform(flam3_genome *cp, int fn, const double4 p, double4 *q_ret, randctx *rc)
 {
    flam3_iter_helper f;