diff --git a/solidify/Makefile b/solidify/Makefile index 1859f7f..5111edc 100644 --- a/solidify/Makefile +++ b/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 diff --git a/solidify/face.h b/solidify/face.h index 2c98e20..d524604 100644 --- a/solidify/face.h +++ b/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 */ diff --git a/solidify/matrix.c b/solidify/matrix.c new file mode 100644 index 0000000..63e90fd --- /dev/null +++ b/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 + +#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]; +} diff --git a/solidify/matrix.h b/solidify/matrix.h new file mode 100644 index 0000000..4b83404 --- /dev/null +++ b/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 */ diff --git a/solidify/overlap.c b/solidify/overlap.c index 0eea1c3..95b2a1a 100644 --- a/solidify/overlap.c +++ b/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; }