Zoom! Demo Source
// (c) Bernie Freidin 1999-2000
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gl/glut.h>
#include "zoom.h"
#include "pack.h"
#define _min3(a,b,c) ((a)<(b)?((a)<(c)?(a):(c)):((b)<(c)?(b):(c)))
#define _max3(a,b,c) ((a)>(b)?((a)>(c)?(a):(c)):((b)>(c)?(b):(c)))
#define PI2 6.28318530717958647692528676655901 // Pi times 2
static packinfo_t *INFO;
static zoom_t *Z;
static int SHADOW_MODE;
static double VIEW_X;
static double VIEW_Y;
static double VIEW_ASPECT;
static pack_t *PACK[2] = {NULL, NULL};
static pack_t A = {0.0, +0.0, -1.0};
static pack_t B = {0.0, +1.0, +2.0};
static pack_t C = {0.0, -1.0, +2.0};
static pack_t *SEARCH_RESULT;
void PACK_Initialize(packinfo_t *packinfo)
{
INFO = packinfo;
INFO->obj_total = 3;
for(int i = 0; i < 3; i++)
{
A.child[i] = NULL, A.depth = -1;
B.child[i] = NULL, B.depth = 0;
C.child[i] = NULL, C.depth = 0;
}
}
void PACK_SetView(double view_x, double view_y)
{
VIEW_X = view_x;
VIEW_Y = view_y;
VIEW_ASPECT = view_x / view_y;
}
pack_t *PACK_Create(pack_t *a, pack_t *b, pack_t *c, int depth)
{
if(c->z > INFO->max_render_z*Z->z) return NULL;
if(c->z < 0.0) return NULL; //FIXME: really necessary?
double ax = a->x, ay = a->y, az = a->z;
double bx = b->x, by = b->y, bz = b->z;
double cx = c->x, cy = c->y, cz = c->z;
// ****************************************************
// Compute next circle (this is a 7th-order equation!!)
// ****************************************************
double dz = az + bz + cz + 2.0*sqrt(az*bz + bz*cz + cz*az);
double vx = ax*bz - az*bx;
double vy = ay*bz - az*by;
double mx = az*(ax*(1.0 - bx*bx + by*by) - 2.0*ay*(bx*by)) +
bz*(bx*(1.0 - ax*ax + ay*ay) - 2.0*by*(ax*ay)) ;
double my = az*(ay*(1.0 + bx*bx - by*by) - 2.0*ax*(bx*by)) +
bz*(by*(1.0 + ax*ax - ay*ay) - 2.0*bx*(ax*ay)) ;
double ex = mx*dz + 2.0*vx*(az-bz);
double ey = my*dz + 2.0*vy*(az-bz);
double fx = az*az*az*bx*(1.0 - bx*bx - by*by) +
bz*bz*bz*ax*(1.0 - ax*ax - ay*ay) ;
double fy = az*az*az*by*(1.0 - bx*bx - by*by) +
bz*bz*bz*ay*(1.0 - ax*ax - ay*ay) ;
double nx = az*bz*ex - dz*fx;
double ny = az*bz*ey - dz*fy;
double s = 4.0*az*bz*sqrt(az*bz + bz*dz + dz*az);
double nz = 2.0*az*bz*(az+bz)*(az+bz);
double dx = (nx - vy*s)/nz;
double dy = (ny + vx*s)/nz;
pack_t *d = (pack_t*)malloc(sizeof(pack_t));
d->x = dx;
d->y = dy;
d->z = dz;
d->parent[0] = a;
d->parent[1] = b;
d->parent[2] = c;
d->child[0] = PACK_Create(a, b, d, depth+1);
d->child[1] = PACK_Create(b, c, d, depth+1);
d->child[2] = PACK_Create(c, a, d, depth+1);
d->depth = depth;
d->error = (fabs(dx-floor(dx)) > 0.1 ||
fabs(dy-floor(dy)) > 0.1 ||
fabs(dz-floor(dz)) > 0.1 ) ? 1 : 0;
INFO->obj_created++;
INFO->obj_total++;
// *************************
// Compute bounding triangle
// *************************
double abx = (ax+bx)/(az+bz);
double aby = (ay+by)/(az+bz);
double bcx = (bx+cx)/(bz+cz);
double bcy = (by+cy)/(bz+cz);
double cax = (cx+ax)/(cz+az);
double cay = (cy+ay)/(cz+az);
// ********************
// Compute bounding box
// ********************
d->box[0] = _min3(abx, bcx, cax);
d->box[1] = _min3(aby, bcy, cay);
d->box[2] = _max3(abx, bcx, cax);
d->box[3] = _max3(aby, bcy, cay);
return d;
}
static void DrawText(char *s)
{
while(*s != '\0')
{
glutBitmapCharacter(GLUT_BITMAP_8_BY_13, (int)(*s++));
}
}
static void DrawCircleVector1(pack_t *pack, double r)
{
char s[10];
if(INFO->show_depth) sprintf(s, "%i", (int)pack->depth);
else sprintf(s, "%i", (int)pack->z);
double len = (double)strlen(s);
if(r*VIEW_Y*Z->z > 12.0*len)
{
double x = r*pack->x;
double y = r*pack->y;
glRasterPos2d(x - 8.0/(VIEW_Y*Z->z)*len,
y - 13.0/(VIEW_Y*Z->z));
DrawText(s);
}
}
static void DrawCircleVector3(pack_t *pack, double r)
{
char s1[10];
char s2[10];
char s3[10];
sprintf(s1, "%i", (int)pack->x);
sprintf(s2, "%i", (int)pack->y);
sprintf(s3, "%i", (int)pack->z);
double len1 = (double)strlen(s1);
double len2 = (double)strlen(s2);
double len3 = (double)strlen(s3);
double len = _max3(len1, len2, len3);
if(len1 + len2 + len3 > 4.0)
{
if(r*VIEW_Y*Z->z > 15.0*len)
{
double x = r*pack->x;
double y = r*pack->y;
double x_loc = x - 8.0/(VIEW_Y*Z->z)*len;
double y_loc = y - 39.0/(VIEW_Y*Z->z);
double dy = 26.0/(VIEW_Y*Z->z);
glRasterPos2d(x_loc + 16.0*(len-len3)/(VIEW_Y*Z->z), y_loc);
DrawText(s3); y_loc += dy;
glRasterPos2d(x_loc + 16.0*(len-len2)/(VIEW_Y*Z->z), y_loc);
DrawText(s2); y_loc += dy;
glRasterPos2d(x_loc + 16.0*(len-len1)/(VIEW_Y*Z->z), y_loc);
DrawText(s1);
}
}
else
{
char s4[40];
sprintf(s4, "%s,%s,%s", s1, s2, s3);
len = (double)strlen(s4);
if(r*VIEW_Y*Z->z > 12.0*len)
{
double x = r*pack->x;
double y = r*pack->y;
glRasterPos2d(x - 8.0/(VIEW_Y*Z->z)*len,
y - 13.0/(VIEW_Y*Z->z));
DrawText(s4);
}
}
}
void PACK_DrawCircle(pack_t *pack)
{
double r = 1.0/fabs(pack->z);
double x = r*pack->x;
double y = r*pack->y;
// **************
// Trivial reject
// **************
if((y - r + Z->y)*Z->z > +1.0) return;
if((y + r + Z->y)*Z->z < -1.0) return;
if((x - r + Z->x)*Z->z > +VIEW_ASPECT) return;
if((x + r + Z->x)*Z->z < -VIEW_ASPECT) return;
double q = 1.0 - pow(r*INFO->max_render_z*Z->z, INFO->fade_exponent);
if(!SHADOW_MODE && pack->z > 0.0)
{
if(pack->error) glColor3d(0.7, 0.3, 0.3);
else glColor3d(0.0, 1.0, 1.0);
if(INFO->show_vec_3)
{
DrawCircleVector3(pack, r);
}
else
{
DrawCircleVector1(pack, r);
}
}
if(pack->error)
{
if(SHADOW_MODE) glColor3d(0.3, 0.0, 0.0);
else glColor3d(1.0, 0.0, 0.0);
}
else
{
//NOTE: in a dark room: use q*0.1, q*0.3
if(SHADOW_MODE) glColor3d(q*0.0, q*0.3, q*0.6);
else glColor3d(q*0.2, q*0.2, q*1.0);
}
// *******************************************************
// Compute dTHETA based on LOD heuristic: radial deviation
// *******************************************************
double err = fabs(pack->z)/(INFO->lod_constant*Z->z);
double dtheta = acos(err>1.0?-1.0 : 1.0 - 4.0*err + 2.0*err*err);
// ***************************************
// Draw circle (this should be tight code)
// ***************************************
double xx = r, x2;
double yy = 0.0;
double c = cos(dtheta);
double s = sin(dtheta);
glBegin(GL_LINE_STRIP);
for(double theta = 0.0; theta <= PI2; theta += dtheta)
{
glVertex2d(x + xx, y + yy);
x2 = xx*c - yy*s;
yy = xx*s + yy*c;
xx = x2;
}
glVertex2d(x + r, y);
glEnd();
if(pack == SEARCH_RESULT)
{
// **********************************
// Draw lines to children and parents
// **********************************
glBegin(GL_LINES);
for(int i = 0; i < 3; i++)
{
if(pack->parent[i])
{
double pr = 1.0/pack->parent[i]->z;
double px = pr*pack->parent[i]->x;
double py = pr*pack->parent[i]->y;
glColor3d(1.0, 1.0, 1.0); glVertex2d(x, y);
glColor3d(0.0, 0.0, 0.0); glVertex2d(px, py);
}
if(pack->child[i])
{
double cr = 1.0/pack->child[i]->z;
double cx = cr*pack->child[i]->x;
double cy = cr*pack->child[i]->y;
glColor3d(0.0, 1.0, 0.0); glVertex2d(x, y);
glColor3d(0.0, 0.0, 0.0); glVertex2d(cx, cy);
}
}
glEnd();
}
INFO->obj_in_view++;
}
void PACK_Draw(pack_t *pack, pack_t *up)
{
if(pack)
{
if(pack->depth > 1)
{
// *************************************
// Trivial reject if circle is too small
// *************************************
if(pack->z > INFO->max_render_z*Z->z)
{
//don't count as node rejection
return;
}
// ************************************
// Trivial reject if box is out of view
// ************************************
double lft = (pack->box[0] + Z->x)*Z->z;
double top = (pack->box[1] + Z->y)*Z->z;
double rgt = (pack->box[2] + Z->x)*Z->z;
double bot = (pack->box[3] + Z->y)*Z->z;
if(lft > +VIEW_ASPECT || top > +1.0 ||
rgt < -VIEW_ASPECT || bot < -1.0 )
{
INFO->node_reject++;
return;
}
}
// ***********************
// Draw circle and recurse
// ***********************
PACK_DrawCircle(pack);
PACK_Draw(pack->child[0], pack);
PACK_Draw(pack->child[1], pack);
PACK_Draw(pack->child[2], pack);
}
else if(up)
{
if(!INFO->dynamic) return;
if(!INFO->dynamic_shadow && INFO->shadow_mode) return;
// ***********************************
// Dynamic generation of more TSP data
// ***********************************
for(int i = 0; i < 3; i++)
{
if(up->child[i]) continue; // already done it
pack_t *a = up->parent[i];
pack_t *b = up->parent[i<2?i+1:0];
up->child[i] = PACK_Create(up, a, b, up->depth+1);
}
}
}
void PACK_Draw(int shadow_mode, zoom_t *zoom)
{
SHADOW_MODE = shadow_mode, Z = zoom;
// **********************
// Root of circle packing
// **********************
INFO->obj_created = 0;
INFO->obj_in_view = 0;
INFO->node_reject = 0;
// ********************************
// Initialize packing, if necessary
// ********************************
if(!PACK[0]) PACK[0] = PACK_Create(&A, &B, &C, 1);
if(!PACK[1]) PACK[1] = PACK_Create(&A, &C, &B, 1);
// **************
// Render packing
// **************
PACK_DrawCircle(&A);
PACK_DrawCircle(&B);
PACK_DrawCircle(&C);
PACK_Draw(PACK[0], NULL);
PACK_Draw(PACK[1], NULL);
}
int PACK_Search(double x, double y, pack_t *pack)
{
if(!pack) return 0;
if(pack->depth > 0)
{
double lft = (pack->box[0] + Z->x)*Z->z;
double top = (pack->box[1] + Z->y)*Z->z;
double rgt = (pack->box[2] + Z->x)*Z->z;
double bot = (pack->box[3] + Z->y)*Z->z;
if(lft > x || top > y || rgt < x || bot < y )
{
return 0;
}
double r0 = Z->z/pack->z;
double x0 = pack->x*r0 + Z->x*Z->z;
double y0 = pack->y*r0 + Z->y*Z->z;
if((x-x0)*(x-x0) + (y-y0)*(y-y0) <= r0*r0)
{
SEARCH_RESULT = pack;
return 1;
}
}
if(PACK_Search(x, y, pack->child[0])) return 1;
if(PACK_Search(x, y, pack->child[1])) return 1;
if(PACK_Search(x, y, pack->child[2])) return 1;
return 0;
}
pack_t *PACK_Search(double x, double y, int show)
{
SEARCH_RESULT = NULL;
if(show)
{
double xd = +(2.0*(double)x/VIEW_X - 1.0)*VIEW_ASPECT;
double yd = -(2.0*(double)y/VIEW_Y - 1.0);
if(PACK_Search(xd, yd, PACK[0])) return SEARCH_RESULT;
if(PACK_Search(xd, yd, PACK[1])) return SEARCH_RESULT;
}
return NULL;
}
void PACK_ClearTSPMemory(pack_t *pack)
{
if(!pack) return;
for(int i = 0; i < 3; i++)
{
PACK_ClearTSPMemory(pack->child[i]);
}
free(pack);
INFO->obj_total--;
}
void PACK_ClearTSPMemory(void)
{
PACK_ClearTSPMemory(PACK[0]);
PACK_ClearTSPMemory(PACK[1]);
PACK[0] = NULL;
PACK[1] = NULL;
if(INFO->obj_total != 0)
{
printf("what???");
}
}
This page © Bernie Freidin, 2000.