(* Copyright (c) 1992 The Geometry Center; University of Minnesota
   1300 South Second Street;  Minneapolis, MN  55454, USA;
   
This file is part of CirclePack. CirclePack is free software; you can
redistribute it and/or modify it only under the terms given in the
file COPYING, which you should have received along with this file.
This and other software may be obtained via anonymous ftp from
geom.umn.edu; email: software@geom.umn.edu. *)

(* Author: Oliver Goodman *)
(* $Id: DrawPacking.m,v 1.4 1992/07/17 16:09:52 oag Exp $ *)

(* Package to enable CirclePacking objects to be displayed. 
   The underlying triangulation is opened out and drawn in
   flat, hyperbolic or spherical geometry depending on 
   Euler characteristic. Disks are coloured to show which 
   are images of the same one. Hyperbolic
   geometry is shown using the Poincare disk model. 
   Spherical geometry may be shown in 3D or stereographically
   projected to the plane at the user's choice. 
*)

BeginPackage["DrawPacking`", {"CirclePack`", "Triangulation`", "Hyperbolic`",
   "Circle3D`"}]

DrawPacking::usage = "DrawPacking[circlepacking] draws a picture of
the packing in whichever is the appropriate geometry. Hyperbolic
geometry is shown using the Poincare disk model. Spherical geometry
is shown using stereographic projection to the plane."

DrawNerve::usage = "DrawNerve[circlepacking] draws the nerve of the
packing in whichever is the appropriate geometry, with the edges 
labeled."

DrawPacking3D::usage = "DrawPacking3D[circlepacking] draws a picture of
a spherical type circle packing on a sphere using 3D graphics."

Begin["`Private`"]

(* acSubgraph finds an acyclic subgraph of a graph which visits
   every vertex. The graph is given as a list of edges where each
   edge is a pair of vertices 
*)
acSubgraph[es_,first_:1]:=Module[{p,r,x,j,sg = {}},
   p[_]=x;
   r[v_]:= If[p[v]===x, v, p[v]=r[p[v]]];
   j[e:{a_,b_}]:= Block[{u=r[a],v=r[b]}, If[u=!=v, p[v]=u; AppendTo[sg,e]]];
   j[es[[first]]];
   Scan[j, es];
   sg
]

(* treeWalk does a depth first walk around a tree (connected acyclic 
   graph). Its arguments are a list of edges from which to start the walk, 
   and a function which is assumed to return all edges emanating from
   a given edge.
*)
treeWalk[start_,c_]:= Scan[treeWalk[c[#],c]&, start]

(* UnionFind with transformations *)
unionFind[r_,es_,group[op_,id_,inv_]]:=Module[{p,root,j},
   p[_]=root;
   r[v_]:= Block[{u,U,w,W},
      If[p[v]===root, {v,id}, 
      {u,U}=p[v]; {w,W}=r[u]; p[v]={w,op[W,U]}]
   ];
   j[{Y_,x_,y_}]:= Block[{u,U,v,V},
      {u,U}=r[x]; {v,V}=r[y];
      If[u=!=v, p[u]={v,op[V,Y,inv[U]]}]
   ];
   Scan[j, es]
]

(* hT[u,#]& is the obvious Mobius transformation which takes u to 0 *)
hT[u_,v_]:= (v - u)/(1 - v Conjugate[u])
Cross[x_,y_] := RotateLeft[RotateLeft[y] x - RotateLeft[x] y]
Rotate3D[v_, w_, aa_]:= (1 - Cos[aa]) (v.w/w.w) w +
   Sin[aa]/Sqrt[w.w] Cross[v,w] + Cos[aa] v

(* relativePt calculates a point in terms of two points, a distance and 
   an angle. 
*)
relativePt[H2[u_],H2[v_],l_,a_]:= 
   hT[-u, Tanh[l/2] Exp[I (a + Arg[hT[u,v]])]] //N //H2
relativePt[E2[u_],E2[v_],l_,a_]:= 
   u + l Exp[I a](v - u)/Abs[v-u] //N //E2
relativePt[S2[u_],S2[v_],l_,a_]:= 
   Rotate3D[Rotate3D[u,Cross[v,u],l],u,a] //N //S2

(* DrawPacking will use the kids function to walk the tree.
   kids lists the (directed) edges pointing away from b, 
   and excluding {b,a}. As a side effect it calculates the 
   positions of all new endpoints.    
*)
kids[p_,an_][{a_,b_}]:= Map[(e={b,#};pos[#]=
   relativePt[pos[b],pos[a],lengths[p][[eN[p[[1]],e]]],an[e]-an[{b,a}]];e)&,
   DeleteCases[star[b],a]
]

getPoints[p_CirclePacking, ce_]:= 
Block[{a,b,c,ca,r,star,acg,S,pos,x,y},
   S = p[[1]];
   acg = acSubgraph[alledges[S],ce];
   (* list the angles between adjacent edges at each vertex *)
   ca = Table[{a,b,c} = corners[S][[i]];{angles[p][[i]], {a,c}, {a,b}}, 
      {i, nc[S]}];
   (* r will tell us the direction of any edge at a vertex *)
   unionFind[r,ca,group[Plus, 0, -#&]];
   star[_]={};
   Scan[({a,b}=#;AppendTo[star[a],b];AppendTo[star[b],a])&,acg];
   (* fix the position of the starting edge *)
   {a,b}=acg[[1]];
   l0 = lengths[p][[eN[S,{a,b}]]];
   {x,y} = Which[euler[S]<0, l0 = Tanh[l0/4]; {H2[-l0], H2[l0]},
      euler[S]==0, l0 = l0/2; {E2[-l0], E2[l0]},
      True, {S2[{1,0,0}], S2[Rotate3D[{1,0,0},{0,1,0},l0]]}
   ];
   pos[a]=x; pos[b]=y;
   (* treewalk around the graph to find the positions of all vertices *)
   treeWalk[{{a,b},{b,a}},kids[p,r[#][[2]]&]];
   (pos /@ allvertices[S]) /. 
      {H2[z_]:> PoincareBall[{Re[z],Im[z]}], E2[z_]:> {Re[z],Im[z]}}
]

(* Each disk gets a random pastel colour *)
randomColour[n_]:= randomColour[n] = Hue[Random[], .5, .5]

(* draw a spherical circle packing using 3D graphics *)
DrawPacking3D[p_CirclePacking, ce_:1]:=
Block[{pt, S = nerve[p]},
   pt=getPoints[p,ce];
   Show[Graphics3D[{
      Line /@ Map[pt[[#]]&, alledges[S], {2}],
      RGBColor[1,0,0],
      Circle[pt[[#]],radii[p][[vN[S, #]]]]& /@ allvertices[S]
      } /. ResolveSpherical3D,
      Boxed -> False
   ]]
] /; euler[nerve[p]]>0

(* draw a spherical circle packing stereographically projected *)
DrawPacking[p_CirclePacking, ce_:1]:=
Block[{pt, S = nerve[p]},
   pt=getPoints[p,ce];
   Show[Graphics[{
      Line /@ Map[pt[[#]]&, alledges[S], {2}],
      Thickness[.0005], RGBColor[1,0,0],
      Circle[pt[[#]],radii[p][[vN[S, #]]]]& /@ allvertices[S]
      } /. ResolveSpherical,
      AspectRatio->Automatic
   ]]
] /; euler[nerve[p]]>0 

(* draw flat and hyperbolic circle packings *)
DrawPacking[p_CirclePacking, ce_:1]:= 
Block[{pt, S = nerve[p]},
   pt=getPoints[p,ce];
   Show[Graphics[{
      If[euler[S]<0,UnitCircle,{}],
      {randomColour[vN[S, #]],Disk[pt[[#]], radii[p][[vN[S, #]]]]}& /@ 
         allvertices[S],
      Line /@ Map[pt[[#]]&, alledges[S], {2}],
      Thickness[.0005],
      Circle[pt[[#]],radii[p][[vN[S, #]]]]& /@ allvertices[S]},
      {AspectRatio->Automatic,
         If[euler[S]<0,PlotRange ->{{-1.,1.},{-1.,1.}},{}]}
   ]]
]

(* a lie but a useful ones for testing purposes *)
Mid[PoincareBall[a_], PoincareBall[b_]]:= PoincareBall[(a + b)/2]

Mid[S2[a_],S2[b_]]:= S2[unit[a+b]]
Mid[a_List, b_List]:= (a+b)/2

dolabel[e:{a_,b_}, S_, pt_]:= Text[eN[S,e], Mid[pt[[a]],pt[[b]]]]

(* draw the nerve of the packing with edge numbers labelled *)
DrawNerve[p_CirclePacking,ce_:1]:=
Block[{pt, S = nerve[p]},
   pt=getPoints[p,ce];
   Show[Graphics[{
      If[euler[S]<0,UnitCircle,{}],
      {RGBColor[0,1,0], Line /@ Map[pt[[#]]&, alledges[S], {2}]},
      dolabel[#,S,pt]& /@ alledges[S]
      } /. If[euler[S]>0, ResolveSpherical, {}],
      {AspectRatio->Automatic,
         If[euler[S]<0,PlotRange ->{{-1.,1.},{-1.,1.}},{}]}
   ]]
]

(* label c with text t inside a small disk *)
diskLabel[t_,c_]:= {{GrayLevel[1],Disk[c,Scaled[{.03,.03}]]},
   Circle[c,Scaled[{.03,.03}]],Text[t,c]}

(* normalize a vector *)
unit[v_List]:= v/Sqrt[N[v.v]]
(* find a set of orthonormal vectors from a given set using Gramm-Schmidt *)
Orthonormal[l_]:= Block[{add, n, prev = {}}, 
   add[v_]:= (AppendTo[prev, n = unit[v - Plus @@ ((#.v #)& /@ prev)]]; n);
   add /@ l]

(* stereographic projection from the sphere to the plane *)
sproj[{x1_,x2_,x3_}]:= If[x3==-1,spole,{x1,x2}/(1 + x3)]

(* more useful junk *)
angle[x_]:= Arg[x[[1]] + I x[[2]]]

(* return the image of a great circle arc under stereographic projn *)
parc[x1_,y1_]:= Block[{c, x = sproj[unit[x1]], y = sproj[unit[y1]], 
      n = unit[Cross[x1,y1]], r, ax, ay},
   If[x==spole, Return[Line[{y, 3 y}]]];
   If[y==spole, Return[Line[{x, 3 x}]]];
   Which[Abs[n[[-1]]]<.01, Return[Line[{x,y}]],
      c = Drop[n,-1]/n[[-1]]; 
      r = 1/Abs[n[[-1]]];
      n[[-1]] > 0, {ax, ay} = {angle[x-c],angle[y-c]},
      True, {ax, ay} = {angle[y-c],angle[x-c]}];
   If[ax>ay, ay += 2 N[Pi]]; 
   Circle[c, r, {ax,ay}]]

(* stuff for drawing 3D pictures in spherical geometry *)
ResolveSpherical3D = {
   Point[S2[x_]]:> Point[x],
   Line[{S2[a_],S2[b_]}]:> Block[{c = ArcCos[a.b], n, bb},
      n = Floor[20 c/Pi] + 1;
      bb = Orthonormal[{a, b}][[2]];
      Line[Table[a Cos[x] + bb Sin[x], {x, 0, c, c/n}]]],
   (ob:Circle|Disk)[S2[v_],r_]:>
      (ob /. {Circle->Circle3D,Disk->Disk3D})[Cos[r] v,Sin[r],v]}

(* stuff for drawing spherical geometry stereographically projected *)
ResolveSpherical = {
   Text[t_, S2[x_]]:> Text[t,sproj[x]],
   Line[{S2[a_], S2[b_]}]:> parc[a,b],
   (ob:Circle|Disk)[S2[{h__,v_}],r_]:> Block[
      {m = -(Cos[r] - v)/((1 - v Cos[r])^2 - {h}.{h} Sin[r]^2)},
      ob[m {h}, Abs[m Sin[r]]]],
   S2[v_]:> Point[sproj[v]]}

SetOptions[Convert, Model->PoincareBall]

End[]
EndPackage[]
