triangle_group_limit_set/main.c

#include <math.h>
#include <stdio.h>
#include <time.h>
#include <sys/time.h>
#include <cairo/cairo-pdf.h>
#include <X11/XKBlib.h>

#include "main.h"
#include "initcairo.h"
#include "triangle.h"
#include "linalg.h"

#define TOGGLE(a) do { (a) = !(a); } while(0)

DrawingContext *screen_context;

// setup everything except cairo and dim, which will be provided by the graphics system
void setupContext(DrawingContext *ctx, int argc, char *argv[])
{
	ctx->n_group_elements = NUM_GROUP_ELEMENTS;
	ctx->p[0] = atoi(argv[1]);
	ctx->p[1] = atoi(argv[2]);
	ctx->p[2] = atoi(argv[3]);
	ctx->k[0] = atoi(argv[4]);
	ctx->k[1] = atoi(argv[5]);
	ctx->k[2] = atoi(argv[6]);
	if(argc > 7)
		ctx->parameter = atof(argv[7]);
	else
		ctx->parameter = 1.0;
//	ctx->parameter = 2.77;
//	ctx->parameter = 0.1;
	ctx->show_boxes = 0;
	ctx->show_boxes2 = 0;
	ctx->show_attractors = 0;
	ctx->show_reflectors = 0;
	ctx->show_rotated_reflectors = 0;
	ctx->show_limit= 0;
	ctx->show_dual_limit= 0;
	ctx->show_text = 1;
	ctx->mode = 0;
	ctx->use_rotation_basis = 0;
	ctx->limit_with_lines = 1;
	ctx->use_repelling = 0;
	ctx->show_marking = 1;
	ctx->marking.x = -0.73679;
	ctx->marking.y = -0.01873;
	ctx->show_coxeter_orbit = 0;

	ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double));
	ctx->limit_curve_count = -1;

	ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t));
	generate_triangle_group(ctx->group, ctx->n_group_elements, ctx->p[0], ctx->p[1], ctx->p[2]);

	// the temporary stuff
	ctx->cartan = gsl_matrix_alloc(3, 3);
	ctx->cob = gsl_matrix_alloc(3, 3);
	ctx->ws = workspace_alloc(3);
}

void destroyContext(DrawingContext *ctx)
{
	free(ctx->limit_curve);
	free(ctx->group);

	gsl_matrix_free(ctx->cartan);
	gsl_matrix_free(ctx->cob);

	workspace_free(ctx->ws);
}

void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
{
	gsl_matrix *tmp = getTempMatrix(ctx->ws);
	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);

//	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");

	initializeTriangleGenerators(gen, ctx->cartan);
	gsl_matrix_set_identity(tmp);
	for(int i = 0; i < strlen(type); i++)
		multiply_right(tmp, gen[type[i]-'a'], ctx->ws);

	rotation_frame(tmp, result, ctx->ws);

	releaseTempMatrices(ctx->ws, 4);
}

void computeBoxTransform(DrawingContext *ctx, char *word1, char *word2, gsl_matrix *result)
{
	vector_t p[2][3],i[2];
	vector_t std[4] = {
		{-1, -1, 1},
		{-1, 1, 1},
		{1, 1, 1},
		{1, -1, 1}
	};

	gsl_vector **vertices = getTempVectors(ctx->ws, 4);
	gsl_vector **std_vertices = getTempVectors(ctx->ws, 4);
	gsl_matrix *tmp = getTempMatrix(ctx->ws);
	gsl_matrix *to_frame = getTempMatrix(ctx->ws);
	gsl_matrix *to_std_frame = getTempMatrix(ctx->ws);

	fixedPoints(ctx, word1, p[0]);
	fixedPoints(ctx, word2, p[1]);

	// intersect attracting line with neutral line of the other element
	for(int j = 0; j < 2; j++)
		i[j] = cross(cross(p[j%2][0],p[j%2][1]),cross(p[(j+1)%2][0],p[(j+1)%2][2]));

	// box consists of p[0][0], i[0], p[1][0], i[1]

	for(int i = 0; i < 4; i++)
		vectorToGsl(std[i], std_vertices[i]);

	vectorToGsl(p[0][0], vertices[0]);
	vectorToGsl(i[0],    vertices[1]);
	vectorToGsl(p[1][0], vertices[2]);
	vectorToGsl(i[1],    vertices[3]);

	projective_frame(std_vertices, to_std_frame, ctx->ws);
	projective_frame(vertices, to_frame, ctx->ws);
	invert(to_frame, tmp, ctx->ws);
	multiply(to_std_frame, tmp, result);

	/*
	LOOP(i) {
		LOOP(j) {
			printf("%.4f ", gsl_matrix_get(result, i, j));
		}
		printf("\n");
		}*/

	releaseTempVectors(ctx->ws, 8);
	releaseTempMatrices(ctx->ws, 3);
}

void updateMatrices(DrawingContext *ctx)
{
	double angle[3];
	LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i];
	cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);

	gsl_matrix *tmp = getTempMatrix(ctx->ws);

	if(ctx->use_rotation_basis % 5 == 0) {
		gsl_matrix_set(tmp, 0, 0, 0.0);
		gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0);
		gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0);
		gsl_matrix_set(tmp, 1, 0, 1.0);
		gsl_matrix_set(tmp, 1, 1, -0.5);
		gsl_matrix_set(tmp, 1, 2, -0.5);
		gsl_matrix_set(tmp, 2, 0, 1.0);
		gsl_matrix_set(tmp, 2, 1, 1.0);
		gsl_matrix_set(tmp, 2, 2, 1.0);
		gsl_matrix_memcpy(ctx->cob, tmp);
	} else if(ctx->use_rotation_basis % 5 == 1) {
		gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason?
	} else if(ctx->use_rotation_basis % 5 == 2) {
		computeRotationMatrix(ctx, tmp, "ba");
		invert(tmp, ctx->cob, ctx->ws);
	} else if(ctx->use_rotation_basis % 5 == 3) {
		computeBoxTransform(ctx, "bca", "abc", ctx->cob);
//		computeBoxTransform(ctx, "cab", "bca", ctx->cob);
//		computeBoxTransform(ctx, "acb", "cba", ctx->cob);
	} else {
		cartanMatrix(tmp, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0);
		diagonalize_symmetric_form(tmp, ctx->cob, ctx->ws);
	}

	releaseTempMatrices(ctx->ws, 1);
}

void output_info(DrawingContext *ctx)
{
	vector_t p[4][3];
	point_t pt;

	fixedPoints(ctx, "abc", p[0]);
	fixedPoints(ctx, "bca", p[1]);
	fixedPoints(ctx, "cab", p[2]);

	pt = vectorToPoint(ctx, p[0][0]);
	printf("(abc)-+ = (%f %f)\n", pt.x, pt.y);
	pt = vectorToPoint(ctx, p[1][0]);
	printf("(bca)-+ = (%f %f)\n", pt.x, pt.y);
}

void print(DrawingContext *screen)
{
	DrawingContext file;
	DimensionsInfo dim;
	cairo_surface_t *surface;

	char filename[100];
	time_t t = time(NULL);
	strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t));

	memcpy(&file, screen, sizeof(file));

	dim.width = screen->dim->width;
	dim.height = screen->dim->width / sqrt(2.0);
	dim.matrix = screen->dim->matrix;
	dim.matrix.y0 += ((double)dim.height - (double)screen->dim->height) / 2.0; // recenter vertically
	updateDimensions(&dim);
	file.dim = &dim;

	surface = cairo_pdf_surface_create(filename, (double)dim.width, (double)dim.height);
    file.cairo = cairo_create(surface);

	draw(&file);

	cairo_destroy(file.cairo);
	cairo_surface_destroy(surface);

	printf("Wrote sceenshot to file: %s\n", filename);
}

int processEvent(GraphicsInfo *info, XEvent *ev)
{
	int state;
	unsigned long key;
	char filename[100];

//	fprintf(stderr, "Event: %d\n", ev->type);

	switch(ev->type) {
	case ButtonPress:
		state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);

		if(ev->xbutton.button == 1 && state & ShiftMask) {
			screen_context->marking.x = (double)ev->xbutton.x;
			screen_context->marking.y = (double)ev->xbutton.y;
			printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
			cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y);
			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
			return STATUS_REDRAW;
		}
		break;

	case KeyPress:
		state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
		key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
		printf("Key pressed: %ld\n", key);

		switch(key) {
		case XK_Down:
			screen_context->parameter /= exp(0.002);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Up:
			screen_context->parameter *= exp(0.002);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Left:
			screen_context->parameter /= exp(0.00002);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Right:
			screen_context->parameter *= exp(0.00002);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Page_Down:
			screen_context->parameter /= exp(0.02);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Page_Up:
			screen_context->parameter *= exp(0.02);
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case ' ':
			screen_context->parameter = 5.57959706;
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case XK_Return:
//			screen_context->parameter = 2.76375163;
			screen_context->parameter = 5.29063366;
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case 'm':
			printf("matrix.xx = %f;\n", info->dim->matrix.xx);
			printf("matrix.xy = %f;\n", info->dim->matrix.xy);
			printf("matrix.x0 = %f;\n", info->dim->matrix.x0);
			printf("matrix.yx = %f;\n", info->dim->matrix.yx);
			printf("matrix.yy = %f;\n", info->dim->matrix.yy);
			printf("matrix.y0 = %f;\n", info->dim->matrix.y0);
			break;
		case 'i':
			output_info(screen_context);
			break;
		case 'b':
			TOGGLE(screen_context->show_boxes);
			break;
		case 'B':
			TOGGLE(screen_context->show_boxes2);
			break;
		case 'a':
			TOGGLE(screen_context->show_attractors);
			break;
		case 'r':
			TOGGLE(screen_context->show_reflectors);
			break;
		case 'x':
			TOGGLE(screen_context->show_rotated_reflectors);
			break;
		case 'L':
			TOGGLE(screen_context->limit_with_lines);
			break;
		case 'l':
			TOGGLE(screen_context->show_limit);
			break;
		case 'd':
			TOGGLE(screen_context->show_dual_limit);
			break;
		case 'R':
			screen_context->use_rotation_basis++;
			updateMatrices(screen_context);
			computeLimitCurve(screen_context);
			break;
		case 'p':
			print(screen_context);
			break;
		case 'M':
			/*
			screen_context->limit_with_lines = 0;
			double parameter_start = screen_context->parameter;
			for(int i = 0; i <= 1300; i++) {
				if(i < 400)
					screen_context->parameter = exp(log(parameter_start)+0.002*i);
				else if(i < 500)
					screen_context->parameter = exp(log(parameter_start)+0.002*400);
				else
					screen_context->parameter = exp(log(parameter_start)+0.002*(900-i));
				updateMatrices(screen_context);
				computeLimitCurve(screen_context);
				draw(screen_context);
				sprintf(filename, "movie3/test%03d.png", i);
				cairo_surface_write_to_png(info->buffer_surface, filename);
				printf("Finished drawing %s\n", filename);
			}
			*/
			screen_context->limit_with_lines = 0;
			double parameter_start = screen_context->parameter;
			for(int i = 0; i <= 1300; i++) {
				if(i < 400)
					screen_context->parameter = exp(0.003*i);
				else if(i < 500)
					screen_context->parameter = exp(0.003*400);
				else
					screen_context->parameter = exp(0.003*(900-i));
				updateMatrices(screen_context);
				computeLimitCurve(screen_context);
				draw(screen_context);
				sprintf(filename, "movie5/test%03d.png", i);
				cairo_surface_write_to_png(info->buffer_surface, filename);
				printf("Finished drawing %s\n", filename);
			}

		case 'f':
			TOGGLE(screen_context->use_repelling);
			computeLimitCurve(screen_context);
			break;
		case 't':
			TOGGLE(screen_context->show_text);
			break;
		case 'c':
			TOGGLE(screen_context->show_coxeter_orbit);
			break;
		case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0':
			screen_context->mode = key - '0';
			break;
		}

		return STATUS_REDRAW;
	}

	return STATUS_NOTHING;
}

int main(int argc, char *argv[])
{
	GraphicsInfo *info;

	screen_context = malloc(sizeof(DrawingContext));
	setupContext(screen_context, argc, argv);
	updateMatrices(screen_context);
	computeLimitCurve(screen_context);

	info = initCairo(0, KeyPressMask, 200, 200, "Triangle group");
	if(!info)
		return 1;

/*
	info->dim->matrix.xx = 274.573171;
	info->dim->matrix.xy = 0.000000;
	info->dim->matrix.x0 = 583.073462;
	info->dim->matrix.yx = 0.000000;
	info->dim->matrix.yy = 274.573171;
	info->dim->matrix.y0 = 777.225293;
*/

	info->dim->matrix.xx = 274.573171;
	info->dim->matrix.xy = 0.000000;
	info->dim->matrix.x0 = 910.073462;
	info->dim->matrix.yx = 0.000000;
	info->dim->matrix.yy = 274.573171;
	info->dim->matrix.y0 = 509.225293;

	updateDimensions(info->dim);

	screen_context->dim = info->dim;
	screen_context->cairo = info->buffer_context;

	startTimer(info);

	while(1) {
		int result = checkEvents(info, processEvent, NULL);
		if(result == STATUS_QUIT)
			return 0;
		else if(result == STATUS_REDRAW) {
			struct timeval current_time;
			double start_time, intermediate_time, end_time;
			gettimeofday(&current_time, 0);
			start_time = current_time.tv_sec + current_time.tv_usec*1e-6;

			draw(screen_context);

			gettimeofday(&current_time, 0);
			intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6;

			cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0);
			cairo_paint(info->front_context);

			gettimeofday(&current_time, 0);
			end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
			printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
		}
		waitUpdateTimer(info);
	}

	free(screen_context);
	destroyCairo(info);
	destroyContext(screen_context);

	return 0;
}