(* :Copyright: 1989-1999, Wolfram Research, Inc. *)

(* :Mathematica Version: Mathematica 2.0 *)

(* :Package Version: 1.1 *)

(* :Title: Shapes of Common 3D Solids *)

(* :Author: Roman Maeder *)

(* :History:
	V1.0 by Roman Maeder
	V1.1 John M. Novak -- added operations for Polyhedra demo, August 96
*)

(* :Keywords:
	Shapes, cylinder, cone, torus, sphere, MoebiusStrip,
	helix, DoubleHelix, wireframe
*)

(* :Requirements: None. *)

(* :Warnings: None. *)

(* :Sources:
	Roman E. Maeder: Programming in Mathematica, 2nd Ed.,
	Addison-Wesley, 1991.
*)

(* :Summary:
This package provides lists of polygons for some common
three-dimensional shapes.  Functions for translating,
rotating, and affine transforming Graphics3D objects and
converting a Graphics3D object to its wire frame 
representation are also included.
*)


BeginPackage["Graphics`Shapes`", "Geometry`Rotations`"]


Cylinder::usage = "Cylinder[(r:1, h:1, (n:20r))] is a list of n polygons
	approximating an open cylinder centered around the z-axis with
	radius r and half height h."

Cone::usage = "Cone[(r:1, h:1, (n:20r))] is a list of n polygons
	approximating a cone centered around the z-axis with
	radius r and extending from -h to h."

Torus::usage = "Torus[(r1:1, r2:0.5, (n:20r1, m:20r2))] is a list of n*m
	polygons approximating a torus centered around the z-axis with
	radii r1 and r2."

Sphere::usage = "Sphere[(r:1, (n:20r, m:15r))] is a list of n*(m-2)+2
	polygons approximating a sphere with radius r."

MoebiusStrip::usage = "MoebiusStrip[(r1:1, r2:0.5, (n:20r1))] is a list
	of 2n polygons approximating a moebius strip centered around
	the z-axis with radii r1 and r2."

Helix::usage = "Helix[(r:1, h:0.5, (m:2, n:20r))] is a list of n*m polygons
	approximating a helix with half height h and m turns."

DoubleHelix::usage = "DoubleHelix[(r:1, h:0.5, (m:2, n:20r))] is a list
	of n*m polygons approximating a double helix with half height h
	and m turns."

RotateShape::usage = "RotateShape[graphics3D, phi, theta, psi] rotates
	the three-dimensional graphics object by the specified Euler angles."

TranslateShape::usage = "TranslateShape[graphics3D, {x, y, z}] translates
	the three-dimensional graphics object by the specified vector."

AffineShape::usage = "AffineShape[graphics3D, {x, y, z}] multiplies
	all coordinates of the three-dimensional graphics object by
	the respective scale factors x, y, and z."

WireFrame::usage = "WireFrame[graphics3D] replaces all
	polygons in the three-dimensional graphics object by outlines."

OutlinePolygons::usage =
"OutlinePolygons[graphics3D, inner, outer] outlines the polygons in the
graphic by beam-like representations. The width of each beam is formed
by creating a hole in each polygon by scaling each point in the polygon
inward by ratio inner. The height of each beam oriented with respect to
the origin is given by scaling the polygon inward by
ratio outer. This works best with planar convex polygons."

ShrinkPolygons::usage =
"ShrinkPolygons[graphics3D, ratio] shrinks the polygons in the graphic
by the given ratio. Each polygon remains centered in the same location,
with its distance of each point from the center multiplied by ratio.
This works best with planar convex polygons."

PerforatePolygons::usage =
"PerforatePolygons[graphics3D, ratio] creates a hole in the middle of each
polygon with radius ratio relative to the distance of each point from
the center of the polygon. Note that each holed polygon is formed by
a number of polygons equivalent to the number of sides of the original
polygon. This works best with planar convex polygons."

Unprotect[Cylinder, Cone, Torus, Sphere, MoebiusStrip, Helix, DoubleHelix,
	RotateShape, TranslateShape, AffineShape, WireFrame]

Begin["`Private`"]

MakeShape[vl_List, c1_Integer, c2_Integer] :=
    Block[{l = vl,
    	   l1 = RotateLeft /@ vl,
    	   mesh},
	mesh = {l, l1, RotateLeft[l1], RotateLeft[l]};
	If[c1 == 1, mesh = Map[Drop[#, -1]&, mesh, {1}] ];
	If[c2 == 1, mesh = Map[Drop[#, -1]&, mesh, {2}] ];
	(*Graphics3D[*) Polygon /@ Transpose[ Map[Flatten[#, 1]&, mesh] ] (*]*)
    ]					/; TensorRank[vl] >= 2
(* c1 = 0 closes the surface in the first dimension, c1 = 1 leaves it open
   and analogous for the second dimension *)

Cylinder[r_?NumericQ, h_?NumericQ, n_Integer] :=
    MakeShape[
	Block[{rcphi, rsphi},
	    Table[rcphi = N[r Cos[2Pi i/n]]; rsphi = N[r Sin[2Pi i/n]];
		  {{rcphi, rsphi, h}, {rcphi, rsphi, -h}}, {i,n}] ],
	0, 1]						/; n>2
Cylinder[r_:1, h_:1] := Cylinder[r, h, Round[20r]]

Cone[r_?NumericQ, h_?NumericQ, n_Integer] :=
	(*Graphics3D[*)N[Table[Polygon[{{r Cos[2Pi i/n], r Sin[2Pi i/n], -h},
			 {r Cos[2Pi (i+1)/n], r Sin[2Pi (i+1)/n], -h},
			 {0, 0, h}}], {i, 0, n-1}]](*]*)/; n > 2
Cone[r_:1, h_:1] := Cone[r, h, Round[20r]]

Torus[r1_?NumericQ, r2_?NumericQ, n_Integer, m_Integer] :=
    MakeShape[
	Block[{cphi, sphi, s},
	    Table[cphi = N[Cos[2Pi i/n]]; sphi = N[Sin[2Pi i/n]];
		  Table[s = N[r1 + r2 Cos[2Pi j/m]];
			{cphi s, sphi s, N[r2 Sin[2Pi j/m]]},
		       {j, m}],
	         {i, n}]],
	0, 0]						/; n>2 && m>2
Torus[r1_:1.0, r2_:0.5] := Torus[r1, r2, Round[20r1], Round[20r2]]

Sphere[r_?NumericQ, n_Integer, m_Integer] :=
    (*Graphics3D[Join[#[[1]], #[[2]]]]& [*)
    MakeShape[
	Block[{cphi, sphi, s},
	  Table[cphi = N[Cos[2Pi i/n]]; sphi = N[Sin[2Pi i/n]];
		Table[s = N[r Cos[-Pi/2+ Pi j/m]];
		      {cphi s, sphi s, N[r Sin[-Pi/2 + Pi j/m]]},
		     {j, 1, m-1}],
	       {i, 0, n-1}]],
	0, 1] ~Join~	(* pole patches *)
    Block[{s = N[r Cos[-Pi/2 + Pi/m]], z = N[r Sin[-Pi/2 + Pi/m]]},
	(*Graphics3D[*){Polygon[Table[{N[s Cos[2Pi i/n]],
			N[s Sin[2Pi i/n]],
			z},
		      {i, 0, n-1}]],
	 	    z=-z;
		    Polygon[Table[{N[s Cos[2Pi i/n]],
	 		N[s Sin[2Pi i/n]],
	 		z},
	 	      {i, n-1, 0, -1}]]
	}(*]*)
    ](*]*)				/; n>2 && m>2
Sphere[r_:1.0] := Sphere[r, Round[20r], Round[15r]]

MoebiusStrip::notes = "We go around it twice, so that shading comes out right."
MoebiusStrip[r1_?NumericQ, r2_?NumericQ, n_Integer] :=
    MakeShape[
	Block[{cphi, sphi, h, rs},
	  Table[rs = N[r2 Cos[Pi i/n]]; h = N[r2 Sin[Pi i/n]];
	        cphi = N[Cos[2Pi i/n]]; sphi = N[Sin[2Pi i/n]];
		{{(r1 + rs) cphi, (r1 + rs) sphi,  h},
		 {(r1 - rs) cphi, (r1 - rs) sphi, -h}},
	       {i,0,2n-1}]],
	0, 1]							/; n>2
MoebiusStrip[r1_:1.0, r2_:0.5] := MoebiusStrip[r1, r2, Round[20r1]]

Helix[r_?NumericQ, h_?NumericQ, m_?NumericQ, n_?NumericQ] :=
    MakeShape[
    	Block[{in, pin, hh = N[2h/m]},
	  Table[in = N[i/n]; pin = N[2Pi in];
	        {{0, 0, hh in}, {r Cos[pin], r Sin[pin], hh in}},
	       {i, -n m / 2, n m / 2}]],
	1, 1]
Helix[r_:1.0, h_:1.0, m_:2] := Helix[r, h, m, Round[20r]]

DoubleHelix[r_?NumericQ, h_?NumericQ, m_?NumericQ, n_?NumericQ] :=
    MakeShape[
    	Block[{rc, rs, ih, hh = N[2h/m]},
	  Table[rc = N[r Cos[2Pi i/n]]; rs = N[r Sin[2Pi i/n]]; ih = hh i/n;
		{{-rc, -rs, ih}, {0, 0, ih}, {rc, rs, ih}},
	       {i, -n m / 2, n m / 2}]],
	1, 1]
DoubleHelix[r_:1.0, h_:1.0, m_:2] := DoubleHelix[r, h, m, Round[20r]]

RotateShape[ shape_, phi_, theta_, psi_ ] :=
    Block[{rotmat = RotationMatrix3D[N[phi], N[theta], N[psi]]},
	shape /. { poly:Polygon[_] :> Map[(rotmat . #)&, poly, {2}],
		   line:Line[_]    :> Map[(rotmat . #)&, line, {2}],
		   point:Point[_]  :> Map[(rotmat . #)&, point,{1}] }
    ]

TranslateShape[shape_, vec_List] :=
    Block[{tvec = N[vec]},
	shape /. { poly:Polygon[_] :> Map[(tvec + #)&, poly, {2}],
		   line:Line[_]    :> Map[(tvec + #)&, line, {2}],
		   point:Point[_]  :> Map[(tvec + #)&, point,{1}] }
    ] /; Length[vec] == 3


AffineShape[ shape_, vec_List ] :=
    Block[{tvec = N[vec]},
	shape /. { poly:Polygon[_] :> Map[(tvec * #)&, poly, {2}],
		   line:Line[_]    :> Map[(tvec * #)&, line, {2}],
		   point:Point[_]  :> Map[(tvec * #)&, point,{1}] }
    ] /; Length[vec] == 3


WireFrame[shape_Graphics3D] := shape /. {Polygon[x_] :> Line[ Append[x, First[x]] ],
		Cuboid[x__] :> cuboidtolines[x]}

WireFrame[shape_] := WireFrame[Graphics3D[shape]]

(* convert a cuboid to a wireframe.  Note that this isn't the fastest
	method; it would be faster to simply enumerate the lines involved,
	as there are only 12. But this is more elegant, makes for a better
	programming example :-)
*)
cuboidtolines[v1_] := cuboidtolines[v1, v1 + {1,1,1}]

cuboidtolines[v1_,v2_] :=
	With[{segs = Outer[List, ##]& @@ Transpose[{v1,v2}]},
		Map[Line,
			Flatten[
				{segs, Transpose[segs,{3,1,2}],
					Transpose[segs,{1,3,2}]},
			2]
		]
	]

(* find the midpoint of a polygon by averaging the points *)
averagepoints[{a_, b___, c_}] := averagepoints[{a,b}]/;a == c

averagepoints[points_] :=
    Map[(Plus @@ #)&, Transpose[points]]/Length[points]

(* scale the points relative to their center by some ratio *)
scalepoints[points_, center_, ratio_] :=
    Map[ratio (#  - center)+ center &, points]

(* shrink polygons by a given ratio; if ratio is >1, polygon expands, if <0,
   mirrors and scales *)
ShrinkPolygons[shape_, ratio_:0.6] :=
    shape/.Polygon[p_] :> Polygon[scalepoints[p, averagepoints[p], ratio]]

(* make holes in polygon of size ratio with respect to center of polygon *)
PerforatePolygons[shape_, ratio_:0.5] :=
    shape/.Polygon[p_] :> perforateaux[p, ratio, True]

perforateaux[{first_, middle___, last_}, r_, e_] :=
    perforateaux[{first, middle, last, first}, r, e]/; first != last

(* if edgeflag is True, polygons are returned with outlining lines,
   if False just the points of the polygons are returned *)
perforateaux[points_, ratio_, edgeflag_] :=
    Block[{newpoints = scalepoints[points, averagepoints[points], ratio],
           polys},
        polys = Map[Flatten[#,1]&,
                    Transpose[{Partition[points,2,1], 
                               Map[Reverse, Partition[newpoints,2,1]]}]
                ];
        If[TrueQ[edgeflag],
            {EdgeForm[], Map[Polygon, polys], Line[points], Line[newpoints]},
            polys
        ]
	]

(* OutlinePolygons -- create 'beams' outlining each polygon. This
   is written with provision for orienting the beams with respect to
   a point other than the origin, but it is currently incompletely
   implemented. One choice might be the average of all points in the
   graphic, which could be found by
   averagepoints[Flatten[Cases[shape, Polygon[p_] :> p, Infinity], 1]]
   Based on an idea/code by Michael Trott, WRI, modified by JMN.
*)
OutlinePolygons[shape_, inner_:0.8, outer_:0.9] :=
    shape/.Polygon[p_] :> outlineaux[p, inner, outer, {0,0,0}]

outlineaux[{first_, middle___, last_}, i_, o_, c_] :=
    outlineaux[{first, middle, last, first}, i, o, c]/; first != last

outlineaux[points_, inner_, outer_, center_] :=
    Block[{outerpoints, outerfunction, innerpoints, radialpoints},
        innerpoints = Apply[{##,#1}&,
                perforateaux[points, inner, False],
            {1}];
        If[center =!= {0,0,0},
            outerfunction = (((outer (# - center)) + center)&),
            outerfunction = ((outer #)&)
        ];
        outerpoints = Map[outerfunction, innerpoints, {2}];
        radialpoints = Map[Flatten[#,1]&,
                    Transpose[{Flatten[Map[Partition[#,2,1]&, innerpoints], 1],
                               Flatten[Map[Reverse,
                                   Map[Partition[#,2,1]&, outerpoints],
                                   {2}
                                ], 1]}]];
        Map[Polygon, Join[innerpoints, outerpoints, radialpoints]]
    ]

End[]   (* Graphics`Shapes`Private` *)

Protect[Cylinder, Cone, Torus, Sphere, MoebiusStrip, Helix, DoubleHelix,
	RotateShape, TranslateShape, AffineShape, WireFrame]

EndPackage[]   (* Graphics`Shapes` *)

(*:Limitations:
*)

(*:Tests:
*)

(*:Examples:

Show[ Graphics3D[ Cylinder[ 0.5,0.5]]]

Show[ Graphics3D[ Cone[]]]

Show[ Graphics3D[ Torus[2,0.7,15,14]]]

Show[ Graphics3D[Sphere[]]]

Show[ Graphics3D[ MoebiusStrip[ 2,1,80]]]

Show[ Graphics3D[ Helix[] ] ]

Show[ Graphics3D[ DoubleHelix[] ] ]

Show[ RotateShape[ Graphics3D[ MoebiusStrip[] ], Pi/4, Pi/3, Pi/2 ] ]

Show[ TranslateShape[ RotateShape[ Graphics3D[ MoebiusStrip[] ], Pi/4, Pi/3,
   Pi/2 ], {2,3,4}] ]

Show[ AffineShape[ Graphics3D[ Cone[] ],{1,2,3} ] ]

Show[ WireFrame[ Graphics3D[ Cone[] ] ] ]

*)
