Program 13 - The Chaos Game

Program Code

(*  :Name:    ChaosGame  *)
(*  :Title:   *)
(*  :Author:  Daniel Mathews  *)
(*  :Context: ChaosGame`  *)
(*  :Package Version:  1.0  *)
(*  :Mathematica Version:  3.0  *)

Off[General::"spell1", General::"spell"];

BeginPackage["ChaosGame`"];

Unprotect[ChaosGame, MyPolygonList];

MyPolygonList = {Speck, Dipole, Triangle, Square, Pentagon, Hexagon, Heptagon, Octagon, Nonagon, Decagon};

ChaosGame::usage = "ChaosGame[start,n] generates a fractal image based upon the following map.
	Begin with a point \!\(P\_0\) start and a set of nodes. Choose a node at random and move to 
	the location halfway between the starting position and the chosen node to obtain the next 
	point \!\(P\_1\).  If this process is repeated and the successive iterations are plotted in 
	the plane, then the resulting figure will be a fractal.  By default, ChaosGame uses as nodes 
	the corners of an equilateral triangle, generating a representation of the Sierpinski triangle.  
	To change the vertex list, use the option Nodes->\!\({{x\_0, y\_0}, {x\_1, y\_1},  ... }\).  
	Additionally, ChaosGame uses a weighted average (with a default value of w=0.5) to calculate 
	the iterations.  This weight may be changed by using the Weight->w option.  ChaosGame also 
	recognizes standard graphics options, which it uses to either change the style of the points 
	(PlotStyle->stylename) or other aspects of the plot (all other options).";

Nodes::usage="The Nodes option of ChaosGame sets the nodes used in the iteration scheme.  The 
	standard usage is of the form Nodes->\!\({{x\_0, y\_0}, {x\_1, y\_1},  ... }\), a list of 
	points in the plane.  Nodes also recognizes the names of selected regular polygons 
	(Triangle through Decagon), so that Nodes->Hexagon will set the nodes to be the corners of 
	a regular hexagon which is inscribed in the unit circle.  The size and orientation of any 
	regular polygon can be changed using the notation <PolygonName>[radius,angle,center], which 
	can be used to change the radius of the circle, the angle of rotation (clockwise, in radians) 
	of the corners about the center, and the coordinates of the center of the polygon.  These 
	polygons may also be incorporated into lists with other points or shapes, allowing the 
	creation of complex lists such as Nodes->{Hexagon,{0,0},Square[0.5,Pi/8],Octagon[2]}.";

Weight::usage="The Weight->w option of ChaosGame sets the weight used to iterate the system.  
	The iteration function uses this value to compute a weighted average of the point \!\(P\_i\) 
	and the randomly-chosen node with weights (1-w) and w, respectively.  Thus a larger weight 
	will cause the points to be clustered more tightly around the nodes than a smaller weight.  
	Weight->list will cause the points to be weighted according to the values in list.  If 
	more weights than nodes are given, then the list of weights will be chopped, and if fewer 
	weights are given, then the list will be padded with values equal to the default weight(1/2).  
	Any weight specified in this manner must be in the range [0, 1], since other values will cause 
    the system to rapidly diverge from the set of nodes.";

Rotation::usage="The Rotation->\!\({a\_1, \ a\_2, \  ... }\) option of ChaosGame changes the 
	iteration function so that the point is rotated by a counterclockwise angle of \!\(a\_i\) 
	around node \!\(n\_i\) after the weighted average has been calculated.  If fewer angles 
	than nodes are given, then the remaining nodes will use no rotation.";

(*  ChaosGame Default Options  *)
Options[ChaosGame] = {
	Nodes -> {{0, 0}, {1, 0}, {0.5, Sqrt[3] / 2.0}},
	Weight -> 1/2,
	Rotation -> None,
	PlotStyle -> {}
};

Begin["`Private`"];

(*  ChaosGame messages  *)

ChaosGame::badw = "Weight for iteration \"`1`\" is not a real number (or list of real numbers) between 0 and 1.";
ChaosGame::badpt = "Initial point \"`1`\" is not a pair of real numbers.";
ChaosGame::badnl = "Node list `1` is not a list of pairs of real numbers.";
ChaosGame::badits = "Number of iterations must be a positive integer (instead of `1`) in position 2.";
ChaosGame::badopt = "Options expected (instead of `1`) beyond position 2 in argument list.";
ChaosGame::badang = "Rotation angle for iteration \"`1`\" is not a real number or list of real numbers.";

(*  Polygon primitives  *)

MyPolygonTest[s_] := Head[s] === Symbol && MemberQ[MyPolygonList, s];
MyPolygon[n_, r_:1, a_:0, c_:{0, 0}] := Sequence @@ Table[
	N[r #1[x + a + If[EvenQ[n], Pi / n, 0]]& /@ {Sin, Cos} + c],
	{x, 0, 2 Pi - 10^-6, 2 Pi / n}
];

MyShapes := {
	name_?MyPolygonTest[args___] -> Hold[MyPolygon[Flatten[Position[MyPolygonList, name]][[1]], args]],
	name_?MyPolygonTest -> Hold[MyPolygon[Flatten[Position[MyPolygonList, name]][[1]]]]
};

RealQ[r_] := NumericQ[r] && Im[r] == 0;
SetAttributes[RealQ,Listable];
RuleQ[r_] := Head[r] === Rule;

(*  The main program  *)

ChaosGame[start_,n_,args___] := Module[{points, len, nodelist, f, w, plotopts, pointopts, i, checklist},
	(*  Check initial point  *)
	If[!MatchQ[start, {_?RealQ, _?RealQ}], Message[ChaosGame::badpt, start]; Return[]];
	(* Check number of iterations  *)
	If[Head[n] =!= Integer || n < 1, Message[ChaosGame::badits, n]; Return[]];
	
	plotopts = {args};
	checklist[l_, test_, str_, mess_] := Module[{i, results = l},
		For[i = 1, i <= Length[results], Null,
			If[!test[results[[i]]],
				Message[mess, results[[i]], str];
				results = Delete[results, i],
				i++
			]
		];
		Return[results]
	];
	SetAttributes[checklist, HoldRest];
	plotopts = checklist[plotopts, RuleQ, ChaosGame::badopt, ""];
		
	(*  Extract settings from argument list  *)
	{w, nodelist, angles, pointopts} = {Weight, Nodes, Rotation, PlotStyle} /. {args} /. Options[ChaosGame];
	nodelist = nodelist /. MyShapes;
	nodelist = ReleaseHold[If[Head[nodelist] =!= List, {nodelist}, nodelist]];
	len=Length[nodelist];
	
	plotopts = DeleteCases[plotopts, (Weight -> _) | (Nodes -> _) | (Rotation -> _) | (PlotStyle -> _)];
	
	plotopts = checklist[plotopts, MemberQ[#[[1]]& /@ Options[Graphics], #1[[1]]]&, ChaosGame::badopt, "graphics "];
	
	(*  Verify node list  *)
	If[Head[nodelist] === Sequence, nodelist = List @@ nodelist];
	If[!MatchQ[nodelist,{{_?RealQ, _?RealQ}..}],
		Message[ChaosGame::badnl, nodelist]; Return[]];
	
	(*  Verify weight  *)
	If[!VectorQ[If[Head[w] === List, w, {w}], (RealQ[#1] && #1≤1 && #1≥0)&],
		Message[ChaosGame::badw, w]; Return[]];
	If[Head[w] === List && Length[w] < len,
		w = Join[w, Table[Weight /. Options[ChaosGame], {len - Length[w]}]]];
	
	(*  Verify rotation angles  *)
	If[!(angles === None || VectorQ[angles, RealQ] || RealQ[angles]),
		Message[ChaosGame::badang, angles]; Return[]];
	If[Head[angles] === List && Length[angles]<len,
		angles = Join[angles, Table[0, {len - Length[angles]}]]];
	
	(*  Create appropriate iteration function (for speed purposes) *)
	If[Head[angles] =!= List,
		If[Head[w] =!= List,
			f[x_] := (1 - w) x + w nodelist[[Random[Integer, {1, len}]]],
			f[x_] := (i = Random[Integer, {1, len}]; (1 - w[[i]]) x + w[[i]] nodelist[[i]])
		],
		angles = N[{{Cos[#1], Sin[#1]}, {-Sin[#1], Cos[#1]}}]& /@ angles;
		If[Head[w] =!= List,
			f[x_] := (i = Random[Integer, {1, len}]; ((1 - w)(x - nodelist[[i]])).angles[[i]] + nodelist[[i]]),
			f[x_] := (i = Random[Integer, {1, len}]; ((1 - w[[i]])(x - nodelist[[i]])).angles[[i]] + nodelist[[i]])
		]
	];
	
	(*  Make values numeric to speed computation  *)
	{nodelist, w} = N[{nodelist, w}];
	
	points=NestList[f,start,n];
	Return[Show[Graphics[Flatten[{pointopts, Point /@ points}, 1]], plotopts, Axes -> True, AspectRatio->1]]
];

End[];

SetAttributes[ChaosGame, ReadProtected];
Protect[chaosgame, ChaosGame, Nodes, Weight, Rotation, MyPolygonList, Evaluate[MyPolygonList]];

EndPackage[];

On[General::"spell1", General::"spell"];

Examples

All default values are used if options are omitted:  nodes at corners of equilateral triangle,weight 0.5, medium point size, and black dots.
This is the Sierpinski Triangle.

ChaosGame[{0, 0}, 1000];

Changing the point size and color.

ChaosGame[{0, 0}, 5000, 
      PlotStyle -> {PointSize[0.001], Hue[0.7]}];

Changing the nodes and weight.

ChaosGame[{0, 0}, 5000, PlotStyle -> PointSize[0.001], 
      
Nodes -> Pentagon, Weight -> 0.6];

Now a more complicated node list with varying weights.

ChaosGame[{0, 0}, 5000, PlotStyle -> PointSize[0.001], 
      
Nodes 
        -> {{0, 0}, Square[1, Pi/6, {0.1, 0.2}], Hexagon[2]}, 
      Weight -> {0.7, 0.81, 0.85, 0.9, 0.85, 0.91, 0.8, 0.84, 0.87, 
          0.89, 0.8}];

A red fractal with hexagonal symmetry.

ChaosGame[{0, 0}, 10000, 
      PlotStyle -> {PointSize[0.001], Hue[0]}, 
      
Nodes -> {{0, 0}, Hexagon}, 
      Weight -> 2/3];

Other plot options can be used as well.

ChaosGame[{0, 0}, 1000, Axes -> False];

Zooming in on the bottom left corner to see the fractal structure in more detail.

ChaosGame[{0, 0}, 10000, PlotRange -> {{0, 0.7}, {0, 0.5}}];

ChaosGame returns a graphics object which can be used later in a  Mathematica session.

p1 = ChaosGame[{0, 0}, 1000, Nodes -> Triangle];

p2 = Graphics[Circle[{0, 0}, 1]];

Show[p2, p1,  AspectRatio -> 1]

ChaosGame can also be used to generate the Cantor Set.

p1 = 
      ChaosGame[{0, 1}, 2000, PlotStyle -> PointSize[0.001], 
        Nodes -> {{0, 1}, {1, 1}}, 
        
Weight -> 2/3, AspectRatio -> 1/5, 
        ImageSize -> {500, 100}, PlotRange -> {{0, 1}, {0, 2}}, 
        Ticks -> {{0, 1, 2, 3, 6, 7, 8, 9}/9, None}, 
        Axes -> {True, False}];

Show[p1, PlotRange -> {{0, 1/3}, {0, 2}}, 
      Ticks -> {{0, 1, 2, 3, 6, 7, 8, 9}/27, None}];

Show[p1, PlotRange -> {{2/27, 1/9}, {0, 2}}, 
      Ticks -> {{2/27, 1/9}, None}];