ben-scans/solidify/overlap.c

/*
 * overlap.c - Overlap two parallel faces
 *
 * 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 <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <gtk/gtk.h>

#include "util.h"
#include "face.h"
#include "solid.h"
#include "style.h"
#include "overlap.h"


#define	UNDEF_F	HUGE_VAL


static int has_osd;


static int sx(const struct solid *s)
{
	return (s->a->sx > s->b->sx ? s->a->sx : s->b->sx)+2*OVERLAP_BORDER;
}


static int sy(const struct solid *s)
{
	return (s->a->sy > s->b->sy ? s->a->sy : s->b->sy)+2*OVERLAP_BORDER;
}


static double r_center(const struct solid *s)
{
	return hypot(sx(s), sy(s))/OVERLAP_CENTER_DIV;
}


static double ramp(int z0, double w0, int z1, double w1)
{
	if (z0 != UNDEF && z1 != UNDEF)
		return z0*w0+z1*w1;
	if (z0 == UNDEF && z0 == UNDEF)
		return UNDEF_F;
	if (z0 == UNDEF && w0 < w1)
		return z1;
	if (z1 == UNDEF && w0 > w1)
		return z0;
	return UNDEF_F;
}


static double zmix(struct face *f, double x, double y)
{
	int xa, xb, ya, yb;
	double zx0, zx1;

	xa = floor(x);
	xb = xa+1;
	ya = floor(y);
	yb = ya+1;

	zx0 = ramp(
	    get_bounded(f->a, xa, ya), yb-y,
	    get_bounded(f->a, xa, yb), y-ya);
	zx1 = ramp(
	    get_bounded(f->a, xb, ya), yb-y,
	    get_bounded(f->a, xb, yb), y-ya);

	return ramp(zx0, xb-x, zx1, x-xa);
}


/*
 * Coordinate transformations, on the example of the x coordinate:
 *
 * - the x coordinate runs from 0 to sx(s)-1
 * - since we work relative to the screen center, this becomes x-sx(s)/2
 *   This is what we perform the coordinate transform on.
 * - 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,
    const struct matrix *ma, const struct matrix *mb)
{
	double za, zb, z;
	double xaf, xbf, yaf, ybf;

	matrix_map(x, y, ma, &xaf, &yaf);
	matrix_map(x, y, mb, &xbf, &ybf);

	za = zmix(s->a, xaf, yaf);
	zb = zmix(s->b, xbf, ybf);

	if (za == UNDEF_F && zb == UNDEF_F)
		return;

	if (za == UNDEF_F) {
		z = 128.0*(zb-s->b->a->min_z)/(s->b->a->max_z-s->b->a->min_z);
		if (z < 0)
			z = 0;
		if (z > 255)
			z = 255;
		p[0] = 255;
		p[1] = z;
		p[2] = z;
		return;
	}
	if (zb == UNDEF_F) {
		z = 128.0*(za-s->a->a->min_z)/(s->a->a->max_z-s->a->a->min_z);
		if (z < 0)
			z = 0;
		if (z > 255)
			z = 255;
		p[0] = z;
		p[1] = 255;
		p[2] = z;
		return;
	}

	z = za;
	za -= face_z0(s->a, xaf, yaf);
	zb -= face_z0(s->b, xbf, ybf);

	if (za+zb < -s->dist) {
		p[0] = 0;
		p[1] = 0;
		p[2] = 255;
		return;
	}

	z = 256.0*(z-s->a->a->min_z)/(s->a->a->max_z-s->a->a->min_z);
	if (z < 0)
		z = 0;
	if (z > 255)
		z = 255;
	p[0] = z;
	p[1] = z;
	p[2] = z;
}


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 } },
		.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, &ma, &mb);
			p += 3;
		}
	gdk_draw_rgb_image(widget->window,
	    widget->style->fg_gc[GTK_STATE_NORMAL],
	    0, 0, sx(s), sy(s), GDK_RGB_DITHER_MAX, rgbbuf, sx(s)*3);
	free(rgbbuf);
}


static void draw_image(GtkWidget *widget, struct solid *s, int osd)
{
	int cx = sx(s)/2;
	int cy = sy(s)/2;
	int p;

	draw_map(widget, s);
	has_osd = osd;
	if (!osd)
		return;
	draw_circle(widget->window, gc_osd, cx, cy, r_center(s));
	p = r_center(s)/sqrt(2);
	gdk_draw_line(widget->window, gc_osd, cx-p, cy-p, cx+p, cy+p);
	gdk_draw_line(widget->window, gc_osd, cx-p, cy+p, cx+p, cy-p);
}


/*
 * Rotate such that a point at distance "r" moves one unit. Rotate
 * counter-clockwise for r > 1, clockwise for r < 0.
 */

static void rotate(struct matrix *m, double r)
{
	struct matrix t;
	double s, c;

	s = 1/r;
	c = sqrt(1-s*s);
	t.a[0][0] = m->a[0][0]*c-m->a[1][0]*s;
	t.a[0][1] = m->a[0][1]*c-m->a[1][1]*s;
	t.a[1][0] = m->a[1][0]*c+m->a[0][0]*s;
	t.a[1][1] = m->a[1][1]*c+m->a[0][1]*s;
	t.b[0] = m->b[0]*c-m->b[1]*s;
	t.b[1] = m->b[0]*s+m->b[1]*c;
	*m = t;
}


static void do_shift(struct matrix *m, double dx, double dy)
{
	m->b[0] += dx;
	m->b[1] += dy;
}


static void shift(struct matrix *m, int dx, int dy, double dist)
{
	/*
	 * Wheeling "up" in each quadrant shifts in the respective direction,
	 * wheeling "down" in the opposite direction.
	 *
	 * No rule without exception: we treat the "down" quadrant like the
	 * "up" quadrant, because it would be extremely counter-intuitive to
	 * wheel "up" to move "down".
	 */

	if (dx > 0 && dy < dx && dy > -dx)
		do_shift(m, dist, 0);
	if (dx < 0 && dy < -dx && dy > dx)
		do_shift(m, -dist, 0);
	if (dy > 0 && dx < dy && dx > -dy)
		do_shift(m, 0, dist);
	if (dy < 0 && dx < -dy && dx > dy)
		do_shift(m, 0, dist);	/* exception ! */
}


static int osd_proximity(const struct solid *s, int dx, int dy)
{
	double r = hypot(dx, dy);
	double rc = r_center(s);

	if (fabs(r-rc) < OSD_PROXIMITY)
		return 1;
	if (r > rc)
		return 0;
	if (abs(abs(dx)-abs(dy)) < OSD_PROXIMITY)
		return 1;
	return 0;
}


static gboolean scroll_event(GtkWidget *widget, GdkEventScroll *event,
    gpointer data)
{
	GtkWidget *darea = gtk_bin_get_child(GTK_BIN(widget));
	struct solid *s = data;
	int dx = event->x-sx(s)/2;
	int dy = event->y-sy(s)/2;
	double r = hypot(dx, dy);
	double rc = r_center(s);
	double rs, rot, dist;
	int center = r < rc;
	int osd = osd_proximity(s, dx, dy);

	if (r < 1)
		return TRUE;

	/*
	 * rot goes exponentially from SLOWEST_ROT*rs to FASTEST_ROT for
	 * r = rc to rs, with rs being half the canvas diagonal.
	 *
	 * The values are picked such that we achieve sufficient precision at
	 * a reasonably large distance from the circle (for accidently entering
	 * the circle would change the mode) but can also spin quickly, e.g.,
	 * when a 180 degrees rotation is needed.
	 *
	 * First, normalize to 0 ... 1
	 * Then, we start at exp(0) and end at
	 * exp(ln(SLOWEST_ROT*rs/FASTEST_ROT)))
	 */
	rs = hypot(sx(s), sy(s))/2;
	rot = (r-rc)/(rs-rc);
	rot = SLOWEST_ROT*rs*exp(-rot*log(SLOWEST_ROT*rs/FASTEST_ROT));

	/*
	 * dist stays at 1 from 0...rc/DIST_STEPS, then linearly goes up to
	 * DIST_STEPS from rc/DIST_STEPS...rc
	 */
	dist = r/rc*DIST_STEPS;
	if (dist < 0)
		dist = 1;

	switch (event->direction) {
	case GDK_SCROLL_UP:
		if (center)
			shift(&s->a->m, dx, dy, dist);
		else
			rotate(&s->a->m, dx > 0 ? rot : -rot);
		draw_image(darea, s, osd);
		break;
	case GDK_SCROLL_DOWN:
		if (center)
			shift(&s->a->m, dx, dy, -dist);
		else
			rotate(&s->a->m, dx > 0 ? -rot : rot);
		draw_image(darea, s, osd);
		break;
	default:
		/* ignore */;
	}
	return TRUE;
}


static gboolean expose_event(GtkWidget *widget, GdkEventExpose *event,
    gpointer user_data)
{
	draw_image(widget, user_data, has_osd);
	return TRUE;
}


static gboolean motion_notify_event(GtkWidget *widget, GdkEventMotion *event,
    gpointer data)
{
	struct solid *s = data;
	int dx = event->x-sx(s)/2;
	int dy = event->y-sy(s)/2;
	int osd = osd_proximity(s, dx, dy);

	if (osd != has_osd)
		draw_image(widget, s, osd);
	return FALSE;
}


void overlap(GtkWidget *canvas, struct solid *s)
{
	GtkWidget *evbox, *darea;

	evbox = gtk_event_box_new();
	darea = gtk_drawing_area_new();

	gtk_widget_set_events(darea,
	    GDK_EXPOSE | GDK_KEY_PRESS_MASK |
	    GDK_BUTTON_PRESS_MASK | GDK_BUTTON_RELEASE_MASK |
	    GDK_SCROLL |
	    GDK_POINTER_MOTION_MASK);

	gtk_widget_set_size_request(darea, sx(s), sy(s));
	gtk_container_add(GTK_CONTAINER(canvas), evbox);
	gtk_container_add(GTK_CONTAINER(evbox), darea);

	draw_image(darea, s, 0);

	g_signal_connect(G_OBJECT(evbox), "scroll-event",
	    G_CALLBACK(scroll_event), s);
	g_signal_connect(G_OBJECT(darea), "expose-event",
	    G_CALLBACK(expose_event), s);
	g_signal_connect(G_OBJECT(darea), "motion-notify-event",
	    G_CALLBACK(motion_notify_event), s);

if (0) {
int i;
long t0 = time(NULL);
gtk_widget_show_all(canvas);
for (i = 0; i != 1000; i++) {
	rotate(&s->a->m, 100);
	draw_image(darea, s, 0);
	while (gtk_events_pending())
		gtk_main_iteration();
}
fprintf(stderr, "%lu\n", time(NULL)-t0);
}

}