#!/usr/local/bin/perl # # apollonius.pl : Draws a 3D Apollonius gasket # (c) 2013 jl_morel@bribes.org - http://http://bribes.org/perl use strict; use warnings; use OpenGL qw/ :all /; use Math::Trig; use Carp; #------ Base vectors my @i = ( 1, 0, 0 ); my @j = ( 0, 1, 0 ); my @k = ( 0, 0, 1 ); my @O = ( 0, 0, 0 ); #------ Circle closest point # ClosestPointC( @P, @C ); # Returns the closest point of the circle @C to the point @P sub ClosestPointC { my @P = @_[0..2]; # point P my @C = @_[3..5]; # circle center my @N = @_[6..8]; # normal vector my $r = $_[9]; # circle radius my @M = CrossProduct( @N, CrossProduct(SubVectors(@P, @C), @N) ); @M = ScaleVector(NormalizeVector(@M), $r); return AddVectors(@C, @M); } #------ Distance from a point to a circle # DistPC( @P, @C ); # Returns the distance between a point @P and a circle @C sub DistPC { my @P = @_[0..2]; # the point my @C = @_[3..9]; # the circle return Dist(@P, ClosestPointC(@P, @C)); } #------ Objective function (cost function) # The minimum of this positive function is zero. # At this minimum, the point P is equidistant of the three circles. sub ObjectiveFunction { my @P = @_[0..2]; # the point my @C1 = @_[3..9]; # the three circles my @C2 = @_[10..16]; my @C3 = @_[17..23]; my $d1 = DistPC(@P, @C1); my $d2 = DistPC(@P, @C2); my $d3 = DistPC(@P, @C3); return ($d1-$d2)**2+($d2-$d3)**2+($d3-$d1)**2; } #------ FindEquid # Being given three pairwise externally tangent circles on the sphere # of center O, and a good guess point, this function returns # a point P equidistant from the three circles. # We use the gradient descent to find the local minimum of the previous # objective function. sub FindEquid { my @P = @_[0..2]; # guess point my @C1 = @_[3..9]; # the three circles my @C2 = @_[10..16]; my @C3 = @_[17..23]; # $theta and $phi are the spherical coordinates of P @P = NormalizeVector(@P); my $phi = acos($P[2]); my $theta = acos($P[0]/sqrt($P[0]*$P[0]+$P[1]*$P[1])); $theta = -$theta if $P[1]<0; my $f = ObjectiveFunction(@P, @C1, @C2, @C3); my $ds = 1e-4; while ( $f > 1e-6 ) { # derivative of f with respect to theta and phi my $dft = (ObjectiveFunction(cos($theta+$ds)*sin($phi), sin($theta+$ds)*sin($phi), cos($phi), @C1, @C2, @C3)-$f)/$ds; my $dfp = (ObjectiveFunction(cos($theta)*sin($phi+$ds), sin($theta)*sin($phi+$ds), cos($phi+$ds) , @C1, @C2, @C3)-$f)/$ds; my $k = 0.1; my ($thetaN, $phiN, $fN); # New (and better!) values of ... $thetaN = $theta - $k*$dft; $phiN = $phi - $k*$dfp; $fN = ObjectiveFunction(cos($thetaN)*sin($phiN), sin($thetaN)*sin($phiN), cos($phiN) , @C1, @C2, @C3); $theta = $thetaN; $phi = $phiN; @P = ( cos($theta)*sin($phi), sin($theta)*sin($phi), cos($phi) ); $f = ObjectiveFunction(@P, @C1, @C2, @C3); } return @P; } #------ FindApolloniusCircle # Returns the Apollonius circle of three pairwise externally tangent circles # on a sphere of center O. sub FindApolloniusCircle { my @C1 = @_[0..6]; # the three circles my @C2 = @_[7..13]; my @C3 = @_[14..20]; # if one of the circles is too small then the Apollonius circle is too # small too and we do not draw small circles! return @C1 if ( $C1[-1] < 0.001 ); return @C2 if ( $C2[-1] < 0.001 ); return @C3 if ( $C3[-1] < 0.001 ); # points of tangency my @T12 = ClosestPointC( @C1[0..2], @C2); my @T13 = ClosestPointC( @C2[0..2], @C3); my @T23 = ClosestPointC( @C3[0..2], @C1); # we take as guess point the centroid of the points of tangency # (the centroid of the centers of the three circles is not a good guess point) my @P = GetBaryCenter( @T12, 1, @T13, 1, @T23, 1); @P = FindEquid(@P, @C1, @C2, @C3); # @P is equidistant of the circles # points of tangency of the Apollonius circle my @P1 = ClosestPointC( @P, @C1); my @P2 = ClosestPointC( @P, @C2); my @P3 = ClosestPointC( @P, @C3); return GetCircumCircle(@P1, @P2, @P3); } #------ Initializations #--- the four vertices of the tetrahedron my @S1 = (1, -1, -1); my @S2 = (1, 1, 1); my @S3 = (-1, -1, 1); my @S4 = (-1, 1, -1); #--- the three base circles my @C0 = GetInCircle(@S1, @S2, @S3); my @C1 = GetInCircle(@S1, @S2, @S4); my @C2 = GetInCircle(@S3, @S2, @S4); my @C3 = FindApolloniusCircle(@C0, @C1, @C2); #--- @gasket is the list of circles to draw my @gasket = ( [@C0], # circle C0 [@C1], # circle C1 [@C2], # circle C2 [@C3] # Apollonius Circle of C0, C1, C2 ); #--- @triads is a list of list of the indices of a triad of circles and # its Apollonius circle. my @triads = ([0, 1, 2, 3]); # the first triad at level 0 #--- %indexcolor is a hash whose keys are indexes where the level changes my %indexcolor; $indexcolor{4} = 1; #------ triad iterator sub nextstep { my @old_triads = @_; my $last_index = 1+$old_triads[-1]->[-1]; $indexcolor{$last_index} = 1; my @next_triads; foreach my $t ( @old_triads ) { push @next_triads, [$t->[0], $t->[1], $t->[3], $last_index++]; push @next_triads, [$t->[1], $t->[2], $t->[3], $last_index++]; push @next_triads, [$t->[2], $t->[0], $t->[3], $last_index++]; } return @next_triads; } #------ Initialization routine my $gasket_corner; # OpenGL display list for one corner of the gasket my $levelmax = 8; # Number of iterations sub init { glClearColor(0, 0, 0, 1 ); # Black background glShadeModel(GL_SMOOTH); # Smooth shading glEnable(GL_MULTISAMPLE); # Enable multisample antialiasing glEnable(GL_DEPTH_TEST); # Enable hidden surface removal # Generating the list of circles for (1..$levelmax) { @triads = nextstep @triads; foreach my $t ( @triads ) { $gasket[$t->[3]] = [ FindApolloniusCircle(@{$gasket[$t->[0]]}, @{$gasket[$t->[1]]}, @{$gasket[$t->[2]]}) ]; } } # Compilation of the display list $gasket_corner = glGenLists(1); glNewList( $gasket_corner, GL_COMPILE ); my $i = 0; # Current index of the circle to draw my $level = 0; # Current level foreach (@gasket) { $level++ if exists $indexcolor{$i++}; glColor3f( Rainbow($level/($levelmax+1)) ); if ( @{$_}[-1] > 0.001 ) { # One does not draw small circles my $r = @{$_}[-1] - 0.001; # Compensation line thickness DrawCircle(@{$_}[0..5], $r, 0.002); } } glEndList(); } #------ Draw the scene my $spin = 0; sub display { glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT ); glLoadIdentity(); gluLookAt( 2.0, 4.0, 10.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ); glPushMatrix(); my $scale = 3.8; glScalef( $scale, $scale, $scale ); glRotatef( $spin, 0.0, 1.0, 0.0 ); # draws the complete gasket using symmetries glCallList($gasket_corner); glRotatef( 180, @i ); glCallList($gasket_corner); glRotatef( 180, @j ); glCallList($gasket_corner); glRotatef( 180, @i ); glCallList($gasket_corner); glPopMatrix(); glutSwapBuffers(); } #------ GLUT Callback called when the window is resized sub reshape { my ( $w, $h ) = @_; glViewport( 0, 0, $w, $h ); glMatrixMode(GL_PROJECTION); glLoadIdentity(); # define the projection gluPerspective( 45.0, $h ? $w / $h : 0, 1.0, 20.0 ); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); } #------ Routine for rotating the scene my $WaitUntil = 0; sub spinDisplay { my $TimeNow = glutGet(GLUT_ELAPSED_TIME); if ( $TimeNow >= $WaitUntil ) { $spin += 1.0; $spin = $spin - 360.0 if ( $spin > 360.0 ); glutPostRedisplay(); $WaitUntil = $TimeNow + 1000 / 25; # 25 frames/s } } #------ GLUT callback for the mouse sub mouse { my ( $button, $state, $x, $y ) = @_; if ( $button == GLUT_LEFT_BUTTON ) { glutIdleFunc( \&spinDisplay ) if ( $state == GLUT_DOWN ); } elsif ( $button == GLUT_RIGHT_BUTTON ) { glutIdleFunc(undef) if ( $state == GLUT_DOWN ); } } #------ Main program glutInit(); glutInitDisplayMode( GLUT_DOUBLE # Double buffering | GLUT_RGB # RGB color mode | GLUT_DEPTH # Hidden surface removal | GLUT_MULTISAMPLE # Multisample antialiasing ); glutInitWindowSize( 300, 300 ); glutCreateWindow("Circles"); init(); glutDisplayFunc( \&display ); glutReshapeFunc( \&reshape ); glutMouseFunc( \&mouse ); glutIdleFunc( \&spinDisplay ); glutMainLoop(); # ====== Utility Functions #------ Add two vectors sub AddVectors { return ( $_[0] + $_[3], $_[1] + $_[4], $_[2] + $_[5] ); } #------ Substract two vector sub SubVectors { return ( $_[0] - $_[3], $_[1] - $_[4], $_[2] - $_[5] ); } #------ Returns the dot product of 2 vectors sub DotProduct { return $_[0] * $_[3] + $_[1] * $_[4] + $_[2] * $_[5]; } #------ Returns the cross product of 2 vectors sub CrossProduct { return ( $_[1] * $_[5] - $_[2] * $_[4], $_[3] * $_[2] - $_[0] * $_[5], $_[0] * $_[4] - $_[1] * $_[3] ); } #------ Returns a normalized vector (length = 1) sub NormalizeVector { return ( 0, 0, 0 ) if ( my $norm = GetVectorLength(@_) ) == 0; return ScaleVector( @_, 1 / $norm ); } #------ Returns the length of a vector sub GetVectorLength { return sqrt( $_[0] * $_[0] + $_[1] * $_[1] + $_[2] * $_[2] ); } #------ Returns the vector scaled by the last parameter sub ScaleVector { return ( $_[0] * $_[3], $_[1] * $_[3], $_[2] * $_[3] ); } #------ Returns the distance between two points sub Dist { return GetVectorLength(SubVectors(@_)); } #------ Returns the barycenter sub GetBaryCenter { croak "bad number of parameters in GetBaryCenter" if @_%4; my ($xG, $yG, $zG, $s) = (0, 0, 0, 0); while ( 0 < @_ ) { my ($x, $y, $z, $k) = splice @_, 0, 4; $xG += $k*$x; $yG += $k*$y; $zG += $k*$z; $s += $k; } return $xG/$s, $yG/$s, $zG/$s; } #------ Returns the unit normal vector of a triangle specified by the three # points P1, P2, and P3. sub FindUnitNormal { return NormalizeVector( CrossProduct( $_[0] - $_[3], $_[1] - $_[4], $_[2] - $_[5], # P1-P2 $_[3] - $_[6], $_[4] - $_[7], $_[5] - $_[8] # P2-P3 ) ); } #------ Draw a circle in space # We draw a circle in space as a torus with a small inner radius # Usage: Drawcircle( $cx, $cy, $cz, $nx, $ny, $nz, $R [, $r]); sub DrawCircle { croak "bad number of parameters in DrawCircle" if 7 > @_; my ( $cx, $cy, $cz, $nx, $ny, $nz, $R, $r) = @_; $r ||= 0.04; glPushMatrix(); my @axe = NormalizeVector CrossProduct(@k, $nx, $ny, $nz ); my $alpha = rad2deg acos(DotProduct($nx, $ny, $nz, @k )); glTranslatef( $cx, $cy, $cz ); glRotatef( $alpha, @axe ); glutSolidTorus( $r, $R, 10, 100 ); # The circle! glPopMatrix(); } #------ Returns the angle ABC = (BA, BC) sub GetAngle { return acos(DotProduct( NormalizeVector(SubVectors(@_[0..2], @_[3..5])), NormalizeVector(SubVectors(@_[6..8], @_[3..5])) )); } #------ Returns the circle circumscribed about a triangle sub GetCircumCircle { croak "bad number of parameters in GetCircumCircle" if 9 != @_; my @A = @_[0..2]; my @B = @_[3..5]; my @C = @_[6..8]; my @N = FindUnitNormal( @A, @B, @C ); # normale my $angleA = GetAngle( @B, @A, @C ); my $angleB = GetAngle( @A, @B, @C ); my $angleC = GetAngle( @A, @C, @B ); # the barycentric coordinates of the center my @G = GetBaryCenter( @A, sin(2*$angleA), @B, sin(2*$angleB), @C, sin(2*$angleC)); return @G, @N, Dist(@G, @A); } #------ Returns the circle inscribed in a triangle sub GetInCircle { croak "bad number of parameters in GetCircumCircle" if 9 != @_; my @A = @_[0..2]; my @B = @_[3..5]; my @C = @_[6..8]; my @N = FindUnitNormal( @A, @B, @C ); # normale my $angleA = GetAngle( @B, @A, @C ); my $angleB = GetAngle( @A, @B, @C ); my $angleC = GetAngle( @A, @C, @B ); # the barycentric coordinates of the center my @G = GetBaryCenter( @A, sin($angleA), @B, sin($angleB), @C, sin($angleC)); # find the radius my $a = Dist(@B, @C); my $b = Dist(@A, @C); my $c = Dist(@A, @B); my $s = ($a + $b + $c )/2; my $r = sqrt( ($s-$a)*($s-$b)*($s-$c)/$s ); return @G, @N, $r; } #------ Rainbow color map function # Usage: ($red, $green, $blue) = Rainbow( $x ); # $x must be between 0 and 1. # Returns the color of the rainbow (RGB list) associated with $x # from blue for $x = 0 to red for $x = 1. sub max { $_[0] < $_[1] ? $_[1] : $_[0] } # max auxiliary function sub Rainbow { my $dx = 0.8; my $s = ( 6 - 2 * $dx ) * $_[0] + $dx; return max( 0, ( 3 - abs( $s - 4 ) - abs( $s - 5 ) ) / 2 ), # Red max( 0, ( 4 - abs( $s - 2 ) - abs( $s - 4 ) ) / 2 ), # Green max( 0, ( 3 - abs( $s - 1 ) - abs( $s - 2 ) ) / 2 ); # Blue }