Programa del copo de Koch en OpenGL


El semestre pasado, que me toco dar clase de Graficación por computadora, en la FESA, un amigo mio (que no tomaba clase conmigo) me pidió que le ayudara en una tarea para su clase de Sistemas Dinámicos, como la verdad no tenia mucho tiempo, hice todo lo posible por zafarme, cuando ya no pude, le pregunte de que se trataba y me dijo que quería hacer un programa gráfico de un fractal, para ilustrar una exposición, que en concreto quería hacer el copo de Koch.

Al principio me gustó la idea, después de todo, traía fresco lo de la graficada y el programa parecía ser muy sencillo, le dije que solo preparara la recurrencia matemática y que después pasaría a ayudarle, después de una semana pase a verlo y como suele suceder no había hecho nada ni yo ni el, sin embargo el programa resulto ser lo bastante simple para que solo me llevara un par de horas terminarlo.

Pueden descargar el código aquí. Y para probar mis últimos dos plugines agregados al blog me decidí a escribir este post. los dos plugines de los que hablo son el WP-sintax que permite poner código de algún lenguaje de programación en tus post y el WP-LaTeX que agregar una fórmula escrita en LaTeX, (en forma de imagen) a tus post.

Copo de Koch

La primeras lineas solo incluyen las dos bibliotecas que vamos a usar, la primera es la glut estandar para levantar OpenGL y la otra es math.h, para hacer cálculos.

Veamos después de eso que hay:

6
7
8
typedef struct {
	float x, y, z;
}Punto;

Con estas lineas definimos un nuevo tipo de dato, lo bautizamos como Punto y tiene adentro tres flotantes, para representar coordenadas, aunque el programa es en esencia bidimensional, como el OpenGL trabaja en 3D, me gusta declarar todo como si estuviera dibujando en 3D.

Después se declaran algunas variables globales, y el prototipo de las funciones, las variables globales son:

  • xMax, yMax: Guardan el tope de las coordenadas del mundo, es decir de donde a donde corre la parte visible de la ventana en X y en Y, con respecto al origen (situada al centro de la ventana)
  • vertices: Es un arreglo de tres puntos, que guardara la única información que necesitamos para empezar a generar el copo de nieve, la posición de tres vértices.
  • numPtos: Este guarda el número de vértices proporcionados por el usuario, es un contador para poner la siguiente lógica: “Si ya tienes tres puntos, empieza a hacer el fractal, si tienes menos de tres puntos, continua esperando…”
  • fondo: Otro entero que guardará la profundidad máxima para hacer el fractal, donde un cero, pinta un triangulo, un uno la estrella de David y así en adelante. Es si lo quieren ver así el numero de veces que se manda a llamar recursivamente el fractal.
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
int main(int argc, char** argv) {
 
   glutInit(&argc, argv);
   glutInitDisplayMode (GLUT_SINGLE | GLUT_RGBA);
   glutInitWindowSize (640, 480);
   glutInitWindowPosition (100, 105);
   glutCreateWindow (argv[0]);
   inicializar();
 
   glutDisplayFunc(dibuja);
   glutReshapeFunc(redimensiona);
   glutKeyboardFunc(teclado);
   glutSpecialFunc(tecladoEspecial);
   glutMouseFunc(raton);
 
   glutMainLoop();
 
   return 0;
}

La función main, no mucho que decir, es muy estándar en el estilo glut, pone la maquina de estados en el estado adecuado, manda llamar a inicializar, registra cinco funciones de callback y por ultimo manda al loop principal de glut. Utilizo un solo buffer por que no voy a hacer animación, y ocupo el nombre del programa para el titulo de la ventana.

47
48
49
50
51
52
void inicializar () {
   glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
   glPointSize(5.0);
   numPtos = 0;
   fondo = 1;
}

La primer función en ser ejecutada, pone el color del fondo en negro, a los puntos les da tamaño 5, inicializa el contador de puntos y a profundidad.

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
void redimensiona (int ancho, int alto) {
   float aspect;
 
   glViewport (0, 0, ancho, alto);
   glMatrixMode (GL_PROJECTION);
   glLoadIdentity ();
 
   if (alto == 0){
   	  alto = 1;
   }
 
   aspect = (float) ancho / (float) alto;
 
   if (ancho <= alto) {
   	  glOrtho(-MUNDO, MUNDO, -MUNDO/aspect, MUNDO/aspect, -MUNDO, MUNDO);
   	  xMax = MUNDO;
   	  yMax = MUNDO / aspect;
 
   } else {
   	  glOrtho(-MUNDO * aspect, MUNDO * aspect, -MUNDO, MUNDO, -MUNDO, MUNDO);
   	  xMax = MUNDO * aspect;
   	  yMax = MUNDO;
   }
 
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();
}

La función redimensiona, es llamada al inicio del programa, (el estándar de glut así lo indica) y cada vez que hay un evento de resize a la ventana. es una función muy común en estilo glut, lo único raro, quizás es que guarda en variables globales el valor de la máxima coordenada Y y la máxima coordenada X. Por lo demás es muy normal, da una proyección ortogonal, i.e: pone el sistema de referencia.

82
83
84
85
86
87
88
89
90
91
92
93
void teclado (unsigned char key, int cx, int cy) {
 
	switch (key) {
		case 27:
		   exit(0);
		break;
 
		default: break;
	}
 
	glutPostRedisplay();
}

La función teclado es muy simple, si detecta un evento de la tecla “Esc”, termina el programa con codigo 0.

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
void tecladoEspecial (int key, int cy, int cx) {
 
	switch (key) {
		case GLUT_KEY_PAGE_UP:
		   fondo++;
		break;
 
		case GLUT_KEY_PAGE_DOWN:
		   if (fondo > 1) fondo--;
		break;
 
		default: break;
	}
 
	glutPostRedisplay();
}

Está función sirve para cachar eventos de las teclas PgUp y PgDn, y en respuesta aumenta el valor de la variable global fondo, es decir hace que el fractal itere mas o menos veces, después de poner el nuevo valor en fondo, avisa a la ventana que necesita volver a dibujarse.

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
void raton (int boton, int estado, int cx, int cy) {
	float x, y;
 
	/* Calculamos las coordenadas de mundo (relativas) apartir de 
	 * las de ventana (absolutas) */
 
	x = ((GLfloat)cx * 2.0 * xMax) / glutGet(GLUT_WINDOW_WIDTH) - xMax;
	y = yMax - ((GLfloat)cy * 2.0 * yMax) / glutGet(GLUT_WINDOW_HEIGHT);
 
	if (boton == GLUT_LEFT_BUTTON && estado == GLUT_DOWN) {
 
		if (numPtos < 3) {
		   vertices[numPtos].x = x;
		   vertices[numPtos].y = y;
		   vertices[numPtos].z = 0.0;
		   ++numPtos;
	    }
	    else {
	       numPtos = 0;
	       vertices[numPtos].x = x;
		   vertices[numPtos].y = y;
		   vertices[numPtos].z = 0.0;
		   numPtos++;
	    }
	}
 
	glutPostRedisplay();
}

Este código cacha las acciones del ratón, las lineas 118 y 119, calculan la posicion en x y en y en donde el usuario realizo el click, recordemos que el glut, nos pasa como parámetros dos enteros que representan el pixel en el que se hizo el click, pero para la lógica de nuestro programa lo queremos en coordenadas de mundo.

En las lineas siguientes, se pregunta si el mouse dio click izquierdo de ser así, preguntamos cuantos puntos llevamos, si son menos de tres, seguimos esperando si no que este sea el primero, de esta manera cada cuarto punto se resetea, el estado del programa, así que si después de pintar un fractal queremos pintar otro simplemente damos un nuevo click y ese es el primer punto de un nuevo fractal.

142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
void dibuja (void) {
	int i;
 
	glClear(GL_COLOR_BUFFER_BIT);
 
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
 
 
	if (numPtos != 3) {
	   glColor3f(1.0, 1.0, 1.0);
	   glBegin(GL_POINTS);
	      for (i = 0; i <numPtos; i++) {
	   	     glVertex3f(vertices[i].x, vertices[i].y, vertices[i].z);
	      }
	   glEnd();
	}
 
	if (numPtos == 3) {
	   glColor3f(1.0, 1.0, 0.0);
	   glBegin(GL_LINE_STRIP);
	      koch(vertices[0], vertices[1], 0);
	      koch(vertices[1], vertices[2], 0);
	      koch(vertices[2], vertices[0], 0);
	   glEnd();
	}	
 
	glFlush();
}

Esta función se encarga de dibujar el fractal, primero limpia la pantalla, luego pone el estado de dibujar y pregunta si tenemos algún numero de puntos que no sea tres (por la lógica de la función anterior esto solo puede pasar cuando son menos que tres) de ser así solo pone los puntos que llevamos. En caso contrario se prepara a dibujar el fractal, para esto hace uso de la función auxiliar koch, que da los vértices de una curva de koch, así que para hacer el copo de koch, pinta tres curvas de koch: Una del vértice 0, al vértice 1, otra del vértice 1 al vértice 2 y por ultimo una que cierra el copo del vértice 2 al vértice 0.

Dibujando la curva de Koch

Ahora viene la parte mas interesante: ¿Como se dibuja la curva de Koch?
La función encargada de hacerlo es la función koch, pide como argumento dos puntos y un entero, que es el nivel de profundidad en la recursividad del algoritmo, es decir: la primera vez que se llame tendrá valor 0, la segunda vez que se llame tendrá valor 1 y así sucesivamente. Aquí esta el código de la función:

173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
void koch (Punto p0, Punto p1, int profundidad) {
	Punto a, b, p2;
 
	if (profundidad < fondo) {
		a.x = 2.0 / 3.0 * p0.x + 1.0 / 3.0 * p1.x;
		a.y = 2.0 / 3.0 * p0.y + 1.0 / 3.0 * p1.y;
		a.z = 0.0;
 
		b.x = 1.0 / 3.0 * p0.x + 2.0 / 3.0 * p1.x;
		b.y = 1.0 / 3.0 * p0.y + 2.0 / 3.0 * p1.y;
		b.z = 0.0;
 
		p2.x = 0.5 * p0.x + 0.5 * p1.x + a.y - b.y;
		p2.y = 0.5 * p0.y + 0.5 * p1.y + b.x - a.x;
		p2.z = 0.0;
 
		glVertex3f(p0.x, p0.y, p0.z);
		koch(p0, a, profundidad + 1);
		glVertex3f(a.x, a.y, a.z);
		koch(a, p2, profundidad + 1);
		glVertex3f(p2.x, p2.y, p2.z);
		koch(p2, b, profundidad + 1);
		glVertex3f(b.x, b.y, b.z);
		koch(b, p1, profundidad + 1);
		glVertex3f(p1.x, p1.y, p1.z);
 
	} 
 
	return;
}

En la linea 176, se pregunta cuantas veces se ha ejecutado esta función en la recursión, si es mas que lo que vale la variable fondo, regresa sin hacer ninguna operación, poniendo fin a la recursión. Sin embargo si aun es menor, se ejecutan los cálculos necesarios (se hace otra iteración) y dentro del procedimiento se llama recursivamente la función koch aumentando el nivel de profundidad.

La lógica de la función koch se explica a continuación:

La idea principal es que como entrada tenemos dos puntos (p0 y p1), los cuales representan un segmento y queremos a partir de ellos calcular tres puntos mas (a, b y p2). Como se ve en la figura:

Segmento de recta

Segmento de recta

Para calcular cada uno de estos puntos, se hace uso de la ecuación paramétrica de un segmento de recta. Si tenemos dos puntos p y q entonces la ecuación paramétrica en t, del semento de recta entre p y q es: s(t) = (1-t)\vec{p} + t\vec{q} con 0 \leq t \leq 1 con esta ecuación podemos calcular fácilmente tanto a, como b. Sin embargo para calcular p2 tendremos que hacer algunas cosas mas.

Puntos a calcular en el nuevo segmento de koch

Puntos a calcular en el nuevo segmento de koch

Primero calculamos el punto medio del segmento, a ese vector hay que sumarle la elevación del nuevo punto, para calcularla necesitamos sumarle un nuevo vector que sea perpendicular al segmento, para esto utilizamos el hecho de que el producto punto de dos vectores perpendiculares entre si, siempre es cero. Así que construimos el vector que va de a a b y de ahí lo volteamos es decir:

Dado un vector:

\left( \begin{array}{c} x  \\ y \end{array} \right)

Es perpendicular a el vector v dado por:

\vec{v} = \left( \begin{array}{c} -y \\ x \end{array} \right)

Por lo que podemos encontrar el punto p2 de la siguiente manera:

\vec{p_2} = \dfrac{1}{2}\vec{p_0} + \dfrac{1}{2}\vec{p_1} + \left( \begin{array}{c} a_{y}-b_{y} \\ b_{x}-a_{x} \end{array} \right)

Recordemos que el vector que va de a a b se puede ver como b-a.

Descargas

¿Por que solo pongo el ejecutable para Windows, si yo siempre uso GNU/Linux?, por que los usuarios del lado luminoso de la fuerza les conviene mas compilar el fuente, que al cabo es solo uno, recuerden que deberán tener instalado OpenGL y glut.

Pantallazos

, , , , , ,

  1. #1 by Anny López on 3 Mayo 2009 - 21:47 pm

    Hola Nemediano, mi nombre es Anny y estoy haciendo un proyecto de investigación acerca de fractales y obviamente de Koch, entre otros, sabes que me demore más que tu haciendo el bendito código, jijijiji :)
    Sabes quisiera hacer un simulador, que dibuje koch, Sierpinski, y Cantor. básicamente los tres.
    Me gustaria por favor me ayudarás :)
    agregame correo@hotmail.com, asi podriamos hablar mejor. Besos y desde ya! muchas gracias :)

  2. #2 by nemediano on 4 Mayo 2009 - 14:55 pm

    Hola Anny, encantado de poder ayudar siempre y cuando se entienda que ayudar, no es lo mismo que hacer.
    los tres fractales que quieres hacer son relativamente fáciles de hacer, dependiendo de a que detalle los quieras.
    edite tu post de arriba, para que no este mi nombre y para que tampoco este tu correo,
    saludos

  3. #3 by mavipasi on 3 Marzo 2010 - 2:27 am

    hola intento hacer funcionar tu codigo pero me sale un error:
    Error 1 error C2381: ‘exit’ : nueva definición; __declspec(noreturn) es diferente c:\archivos de programa\microsoft visual studio 9.0\vc\include\stdlib.h 371 PROYECCIONES
    me puedes ayudar a arreglarlo gracias.

  4. #4 by nemediano on 3 Marzo 2010 - 12:33 pm

    Hola! Se refiere a que usas un compilador que no sigue el estandar (el de visual estudio me imagino). Dado que el ansi c, dice claramente que exit() esta en stdlib. Y se diefie como yo la uso void exit (int error_code)
    Puedes intentar compilarlo comentando la llamada a exit en la linea 86 del código y luego comentar también la inclusion de la libreria stdlib.

(No será publicado)