This file contains self-contained code for generating each piece of art in the article.

Figure : Alphabet Example

font="Times";
ToGraphic[n_String,font_,opts___]:=ToGraphic[n,font]=(grays=ImageData[
ColorConvert[Rasterize[Graphics[Text[Style[ToString[n],50,FontFamily->font],{0,0}]],RasterSize->200,ImageSize->200],"GrayScale"]];L=ListContourPlot[grays,Contours->{0.5},ContourShading->False,AspectRatio->Automatic,ContourStyle->Black,PlotRange->All];cc=Cases[Normal[L],Line[z_]:>Polygon[z],∞]/.Polygon[z_]:>Polygon[z[[Range[1,Length[z],5]]]]/.{x_Real,y_}:>{x,-y}/.{Polygon[z_],Polygon[u_]}->{Polygon[z],Gray,Polygon[u]}/.{Polygon[z_],Polygon[u_],Polygon[v_]}->{Polygon[z],Gray,Polygon[u],Polygon[v]};pts=Cases[cc,{x_,y_Real},∞];{{xmin,xmax},{ymin,ymax}}={{Min[First/@pts],Max[First/@pts]},{Min[Last/@pts],Max[Last/@pts]}};Graphics[cc/.{x_,y_}:>(1/35) {x-100,y+96.5},opts,PlotRange->All])

Clear[polyPts, linePts];
f[x_, y_]:= {((((x^2)+(y^2)))^(1/4)) Cos[(1/2) Arg[x+I y]], ((((x^2)+(y^2)))^(1/4)) Sin[(1/2) Arg[x+I y]]};
g[x_, y_]:= {-((((x^2)+(y^2)))^(1/4)) Cos[(1/2) Arg[x+I y]], -((((x^2)+(y^2)))^(1/4)) Sin[(1/2) Arg[x+I y]]};

polyPts[tt_]:=Join @@ Table[
Table[(1-t) tt[[i]]+t tt[[i+1/. Length[tt]+1->1]],{t,0,.99,.01}],{i,1,Length[tt]}]
linePts[tt_]:=Line[Join @@ Table[
Table[(1-t) tt[[i]]+t tt[[i+1]],{t,0,.99,.01}],{i,1,Length[tt]-1}]]
st={.1,.1};
{nresx,nresy}={6,4};
res=.1;st={.1,.1};
cols =Lighter/@Lighter/@Table[ ColorData["Rainbow"][x], {x,0, 1, 1/(nresx nresy-1)}];
colsX=White&/@ Range[100];
letters=Characters["ABCDEFGHIJKLMNOPQRSTUVWX"];

SortByPolarFromCentroid[pts_]:=(cen=Mean[pts];
perm=SortBy[Range[Length[pts]],
 N@Arg[Complex@@((pts[[#]]-cen))]&];
pts[[perm]])

basicGrid=Graphics[Table[{EdgeForm[Black],FaceForm[cols[[(ij[[2]]) nresx + ij[[1]]+1]]],
 Polygon[polyPts[#+st+{ij[[1]],-ij[[2]]+nresy-1}&/@ {{-nresx/2,-nresy/2},{-nresx/2+1,-nresy/2},{-nresx/2+1,-nresy/2+1},{-nresx/2,-nresy/2+1}}]]},{ij, Tuples[{Range[0,nresx-1],Range[0,nresy-1]}]}]];

whiteData={EdgeForm[White], FaceForm[White],Polygon@Table[.09{Cos[t], Sin[t]}, {t,0, 2π, 2π/100}]};
blackData={EdgeForm[Black], FaceForm[Black],Polygon@Table[.012{Cos[t], Sin[t]}, {t,0, 2π, 2π/100}]};

transformedGrid = Graphics[Table[{EdgeForm[Black],FaceForm[cols[[(ij[[2]]) nresx + ij[[1]]+1]]],
 pp=polyPts[#+st+{ij[[1]],-ij[[2]]+nresy-1}&/@ {{-nresx/2,-nresy/2},{-nresx/2+1,-nresy/2},{-nresx/2+1,-nresy/2+1},{-nresx/2,-nresy/2+1}}];
Which[ij=={2,2},
pp1=Flatten[pp/. {x_, y_?NumericQ}:>{ f[x,y], g[x,y]},1];
pp1=SortByPolarFromCentroid[pp1];Polygon[pp1],

!MemberQ[{{0,2},{1,2}}, ij],{Polygon[pp/. {x_, y_?NumericQ}:> f[x,y]],Polygon[pp/. {x_, y_?NumericQ}:> g[x,y]]},

MemberQ[{{0,2},{1,2}}, ij], 
 
ppA=Take[ppT=polyPts[{{-ij[[1]]-2+st[[1]], -.001},{-ij[[1]]-2, -1}+st,{-ij[[1]]-1, -1}+st,{-ij[[1]]-1+st[[1]], -.001} }],300];

ppB=Take[polyPts[{{-ij[[1]]-1+st[[1]], .001},{-ij[[1]]-1, 0}+st,{-ij[[1]]-2, 0}+st,{-ij[[1]]-2, 0}+st , {-ij[[1]]-1+st[[1]], .001}}],300];
ppA1=ppA/.{x_, y_?NumericQ}:> f[x,y];
ppB1=ppB/.{x_, y_?NumericQ}:> g[x,y];
pp1=Join[ppA1,ppB1];{Polygon[pp1],Polygon[pp1/. {x_,y_}:>{-x,-y}]}],
PointSize[.01],EdgeForm[],FaceForm[White],
whiteData;Point[{0,0}];Line[{{{-2,0}, {2,0}},Reverse/@ {{-2,0}, {2,0}}}];{}},{ij, Tuples[{Range[0,nresx-1],Range[0,nresy-1]}]}],Frame->!True, FrameTicks->Automatic]

transformedWhiteData={whiteData/. {x_,y_?NumericQ}:>f@@{x,y},whiteData/. {x_,y_?NumericQ}:>g@@{x,y}};
transformedBlackData={blackData/. {x_,y_?NumericQ}:>f@@{x,y},blackData/. {x_,y_?NumericQ}:>g@@{x,y}};
letterData=Table[ToGraphic[letters[[ij[[2]] nresx+ij[[1]]+1]],font][[1]]/. Gray->cols[[ij[[2]] nresx+ij[[1]]+1]]/. {x_Real,y_}:>((2 {x,y})/4)+{.6+ij[[1]]-(nresx/2),1.6-ij[[2]]+(nresy/2)-2},{ij,Tuples[{Range[0,nresx-1],Range[0,nresy-1]}]}];

transformedGridAndLetters=Show[transformedGrid,Graphics[{
letterData/. {x_, y_?NumericQ}:> f@@{x,y},letterData/. {x_, y_?NumericQ}:> g@@{x,y}}], Frame->True, PlotRange->1.9, FrameTicks-> {{{},{}}, {{0, 1.732,1.8},{}}}, GridLines->{{},{}}]

GraphicsRow[{Show[basicGrid,Graphics[{letterData,PointSize[.01],FaceForm[GrayLevel@1],Disk[{0,0}, .14 ] ;whiteData,Point[{0,0}];blackData}], PlotRange-> {{-3.1, 3.3}, {-2.1, 2.3}},Frame->True,ImageSize->200 1.3,Axes->False,PlotRegion->{{0,1}, {0,.95}},
FrameTicks->{{Range[-1,2],{}},{Range[-2,3],{}}}], Show[Graphics[Rectangle[.95{-2,-2}, .95 {2,2}]],transformedGridAndLetters, 
Graphics[{transformedWhiteData, Point[{0,0}];transformedBlackData}],Frame->True,PlotRange-> 2.0,ImageSize->1.2 150, Axes->False,
PlotRegion->{{0,1}, {0,.95}},FrameTicks-> {{Join[Range[-1,2], {}],{}},{Join[Range[-2,2], {}],{}}}]}, Spacings->-30]

Figure : Statue of Liberty via Quadratic Camera

scale = 1;

pic=;
dims=ImageDimensions[pic]
pic=ImageResize[pic, scale dims];
dims=ImageDimensions[pic]

complexSize=2.4 ;

 im=Image[Rasterize[ImagePad[pic,{{0,0},{0,0}}],RasterSize->(800;Automatic), ImageResolution->300]];
im=ImageResize[im , {3558,2668}];
ImageDimensions[im]

Quiet[LaunchKernels[]];

c=0;  p=0.;
limitColor = {0.15,0.95,0.95};
{xmin,xmax}={-complexSize,complexSize};
{ymin,ymax}={-complexSize,complexSize};
{m,n}=ImageDimensions[im]
imageData = ImageData[im];

color1 = With[{c=c,imageData=imageData,m=m,n=n,p=p,limitColor=limitColor,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax},Compile[{{z,_Complex}},Module[{z1},z1=(z^2)+c;If[!(xmin≤Re[z1]≤xmax&&ymin≤Im[z1]≤ymax),{0,0,0},p limitColor+(1-p) imageData[[Round[Max[1,((n (Im[z1]-ymax))/(ymin-ymax))]],Round[Max[1,((m (Re[z1]-xmin))/(xmax-xmin))]]]]]],RuntimeAttributes->{Listable},Parallelization->True]];
colors=color1[Table[x+I y,{y,ymax,ymin,((ymin-ymax)/n)},{x,xmin,xmax,((xmax-xmin)/m)}]];

im = ImageTake[Image[colors], {250,2400}, {400, 3100}];
Show[im, ImageSize->400]

Code for Figure : Julia Inverse Iteration Algorithm

c = -0.12256+0.74486 I;c=.32+.043 I;c=.15+.43 I;

cols={Red,Green,Lighter@Blue,Yellow,Purple,Cyan,Orange, Pink, Brown, Lighter@Blue, Red, Lighter@Green, Lighter@Orange, Lighter@Cyan};

Clear[invImage, invImage2, reducedInvImage];
invImage[z_?NumericQ] := {1,-1}Sqrt[z-c];
invImage[zs:{__?NumericQ}] := Flatten[invImage /@ zs];
cPts = NestList[invImage,{N[1+I]},15];
res=200;
invImage2:=(N[{1,-1} Floor[res Sqrt[#-c]]/res]&/@#//Flatten)&;
reducedInvImage[points_]:=Module[{newPoints},newPoints=Complement[invImage2[points],pointsSoFar];
pointsSoFar=Union[newPoints,pointsSoFar];
newPoints];
pointsSoFar=Nest[invImage2,N[{1+I}],15];
FixedPoint[reducedInvImage,pointsSoFar];
cPts = Join[cPts,{pointsSoFar}];
pts = Take[Map[{Re[#],Im[#]}&,cPts,{2}],10];
Graphics[{Table[
{{EdgeForm[{Which[i≤3, Black,i<50,cols[[i]]//Darker,True, White],
Which[i==1,Thickness[.005], i≤3, Thickness[.005], 
i≤7, Thickness[.004], True, Thickness[.008]]}],FaceForm[{cols[[i]]}],
Disk[#, Which[i==1, .06,i==2, .04, i≥2,.1/i^1.2]]&/@ pts[[i]]}},{i,1,Length[pts]}],
Text[Style["1",FontSize->18, White],pts[[1,1]]];
Text[Style["2",FontSize->14, White],pts[[2,1]]];
Text[Style["2",FontSize->14, White],pts[[2,2]]];{}},Frame->True, FrameTicks->{{{0, {1,I},{-1,"-I"}},{}},{{0,1,-1},{}}}, PlotRange->{{-1.4,1.4}, {-1.25, 1.25}},ImageSize->280]

Code for Figure : Filled Julia Set

The JuliaSet Package (from S. Wagon, Mathematica in Action, 3rd Edition, SpringerChapter 11)

(* :Name: JuliaSet` *)
(* :Author: Mark McClure, April 1998 *)  (* some additions by Stan Wagon, August, 2007 *)
(* :Copyright: Copyright Mark McClure, 1998 *)
(* :Package Version: 1.1 *)
(* :Mathematica Version: 2.2 *)
(* :Summary:
	Generates images of Julia sets using several
	variants of the inverse iteration algorithm. 
*)

BeginPackage["JuliaSet`"];


InitialIterations::usage = "InitialIterations is an option to JuliaModified that sets the number of times the inverse iteration method is run.";


ComplexTicks::usage = "ComplexTicks[s] turns the complex number (or list) s into specifications for complex tick marks.";

TraceQ::usage = "TraceQ is an option to JuliaModified that prints the iteration number.";
TraceStep::usage = "";
JuliaSet::usage = "JuliaSet is a package containing several
implementations of inverse iteration algorithms for the
generation of Julia sets.";

Julia::usage = "Julia[formula,var] generates the Julia set of 
the rational function defined by formula, using the modified 
inverse iteration algorithm.";
 
JuliaSimple::usage = "JuliaSimple[c,depth] generates 
2^depth points of the Julia set for z^2 + c using a simple 
deterministic algorithm.";
 	
JuliaModified::usage = "JuliaModified[c] generates the Julia set 
for z^2 + c using a modified determinstic algorithm.";

JuliaStochastic::usage = "JuliaStochastic[c,numpoints] generates 
numpoints points of the Julia set for z^2 + c using a stochastic
algorithm.";

FilledJuliaSet::usage = "FilledJuliaSet[c, {x0,x1},{y0,y1}, opts] generates a filled Julia set for z^2 + c using an escape time algorithm.";

Resolution::usage = "Resolution is a positive integer valued 
Option of Julia and JuliaModified which controls the level of 
detail of the generated set.";

Bound::usage = "Bound is a positive real valued Option of 
Julia which specifies the maximum absolute vale of points
plotted in the Julia set of a rational function.";

ColorQ::usage = "ColorQ is an option to FilledJuliaSet that chooses between a colored or a gray color scheme.";

Iterations::usage = "Iterations is an option to FilledJuliaSet that sets the number of iterations to try.";


OneBias::usage = "OneBias is an Option of JuliaStochastic which 
determines the probability of choosing the positive inverse of 
z^2 + c.  OneBias must be a real number stricly between 0 and 1.";

JuliaSetOnlyQ::usage = "JuliaSetOnlyQ is an option to FilledJuliaSet that, when True, causes only the black Julia set to appear.";
 	
Begin["`Private`"];

Attributes[ComplexTicks] = Listable;

ComplexTicks[s_?NumericQ] := 
  {s, s I} /. Thread[{-1. I, 0. I, 1. I} ->  {-I, 0, I}];


(* fR is like FilterOptions *) 

fR[sym_, opts___] := Sequence @@ FilterRules[{opts}, First/@Options[sym]];

(* Set the Options *)

(* For JuliaSimple *)
JuliaSimpleOptions = Join[{AspectRatio -> Automatic, Axes -> False,
		PlotStyle -> {PointSize[0.0005], Black}}, Options[ListPlot]];
keywords = Complement[Union[First /@ JuliaSimpleOptions], 
	{DisplayFunction, DefaultFont}];
vals = keywords /.  JuliaSimpleOptions;
special = {DisplayFunction :> $DisplayFunction, DefaultFont :> $DefaultFont};
Options[JuliaSimple] = Union[Apply[Rule,Transpose[{keywords,vals}],{1}],
	special];

(* For JuliaModified *)
Options[JuliaStochastic] = Join[{OneBias -> .5}, Options[JuliaSimple]];
Options[JuliaModified] = Join[{Resolution -> 200, TraceQ->False,TraceStep->20, InitialIterations->10}, Options[JuliaSimple]];
Options[Julia]    = Join[{Bound -> 4}, Options[JuliaModified]];


(* The functions *)

(* Error Messages *)
JuliaSet::argu = "`1` called with `2`; `3` expected.";
JuliaSet::badNumber = "Number expected at position 1 in `1`.";
JuliaSet::badInt = "Positive Integer expected at position 2 in `1`.";
JuliaSet::badOpt = "Options expected as optional arguments in `1`.";
JuliaSet::optx = "Unknown Option `1`.";
Julia::badFormula = "Rational function of degree at least two 
expected at position 1 in `1`.";
Julia::noInverses = "Can't find the the inverse functions.";
OneBias::outOfRange = "OneBias must be stricly between 0 and 1. Using 0.5.";
Resolution::outOfRange = "Resolution must be a positive integer. Using 200.";
Bound::outOfRange = "Bound must be a positive number.  Using 4.";


JuliaSimple[c_, depth_, opts___] := Module[
	{invImage, points, valid = First /@ Options[JuliaSimple]},
	
	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		Flatten[{opts}]
	];
	
	invImage := 
		( N[{1,-1} Sqrt[#-c]]& /@ # // Flatten ) &;
	points = {Re[#], Im[#]}& /@ Nest[invImage, {1}, depth];
	ListPlot[points, 
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaSimple]]
	]
] /;  
  ( NumberQ[N[c]] || Message[JuliaSet::badNumber, JuliaSimple, c] ) &&
	( (IntegerQ[depth] && depth > 0 ) ||
	 Message[JuliaSet::badInt, JuliaSimple, depth] ) &&
	( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaSimple] )

(* Special cases for JuliaSimple *)
JuliaSimple[]   := Message[JuliaSet::argu, "JuliaSimple", 
	"0 arguments", "2 arguments are"]
JuliaSimple[c_] := Message[JuliaSet::argu, "JuliaSimple", 
	"1 argument",  "2 arguments are"]


JuliaModified[c_, opts___] := Module[
	{cc,invImage,reducedInvImage,pointsSoFar,res, 
	valid = First /@ Options[JuliaModified]},

	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}];
	{res, trQ, ninit,step} = {Resolution, TraceQ, InitialIterations, TraceStep} /. {opts} /. Options[JuliaModified];
  cc=ninit;
        If[!(IntegerQ[res] && res>0),
		{Message[Resolution::outOfRange];
		res = Resolution /. Options[JuliaModified]}];
	
	invImage := 
		( N[{1,-1} Floor[N[res Sqrt[#-c]]]/res]&  /@ # // Flatten )&;

	reducedInvImage[points_] := Module[{newPoints},
		newPoints = Complement[invImage[points], pointsSoFar];
		pointsSoFar = Join[newPoints, pointsSoFar];
		If[trQ && Divisible[cc,step], Print[StringForm["`` levels, `` points, `` new points",cc,Length[pointsSoFar],Length[newPoints]]]];
		newPoints];  (* Mark used Union here *) 
    pointsSoFar = Nest[invImage, {1}, ninit];
FixedPoint[(cc++;reducedInvImage[#])&, pointsSoFar];

If[trQ,Print[StringForm["`` levels and `` points",cc,Length[pointsSoFar]]]];
	ListPlot[{Re[#], Im[#]}& /@ pointsSoFar, 
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaModified]]
	]] /;  
	(NumberQ[c] || Message[JuliaSet::badNumber, JuliaModified, c]) &&
    ( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaModified] )

(* Special case for JuliaModified *)
JuliaModified[] := Message[JuliaSet::argu, "JuliaModified", 
	"0 arguments", "1 argument is"]
	
	

Julia[formula_, var_Symbol, opts___] := Module[
	{z0,inverses,inverseRules,funcs,image,
	reducedImage,pointsSoFar, res,bound,
	valid = First /@ Options[Julia]},

	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}
	];

	res = Resolution /. {opts} /. Options[Julia];
	If[!(IntegerQ[res] && res>0),
		{Message[Resolution::outOfRange];
		res = Resolution /. Options[Julia]}];
	bound = Bound /. {opts} /. Options[Julia];
	If[!(bound>0),
		{Message[Bound::outOfRange];
		bound = Bound /. Options[Julia]}];
	
	inverseRules = NSolve[#==formula,var];
	If[And @@ ListQ /@ inverseRules, If[Length[inverseRules] > 1,
	
		inverses = var /. inverseRules;
		funcs = Function[anInverse,N[Floor[N[anInverse res]]/res] &] 
			/@ inverses;
		image = Flatten[Through[funcs[#]]& /@ #,1]&;
		If[PolynomialQ[formula,var],
			reducedImage[points_] := Module[{newPoints},
				newPoints = Complement[image[points], pointsSoFar];
				pointsSoFar = (* Union *) Join [newPoints, pointsSoFar];
				newPoints],
			reducedImage[points_] := Module[{newPoints},
				newPoints = Complement[image[points], pointsSoFar];
				newPoints = Select[newPoints,(N[Abs[#]]<=bound)&];
				pointsSoFar = Union[newPoints, pointsSoFar];
				newPoints]
		];
		z0 = N[Nest[funcs[[1]][funcs[[2]][#]]&, N[Pi],10]];
		pointsSoFar = Select[Nest[image,{z0},3],(N[Abs[#]]<=bound)&];
		FixedPoint[reducedImage, pointsSoFar];
		ListPlot[{Re[#],Im[#]}& /@ pointsSoFar, 
			fR[ListPlot,opts], 
			fR[ListPlot,Sequence @@ Options[JuliaModified]]
		],
		Message[Julia::badFormula, Julia]], Message[Julia::noInverses] 
	]
] /; (NumberQ[N[formula /. var -> Pi]]                &&
	PolynomialQ[Numerator[Together[formula]], var]    &&
	PolynomialQ[Denominator[Together[formula]], var]) || 
	Message[Julia::badFormula, Julia]

(* Special cases for Julia *)
Julia[]   := Message[JuliaSet::argu, "Julia", 
	"0 arguments", "2 arguments are"]
Julia[c_] := Message[JuliaSet::argu, "Julia", 
	"1 argument",  "2 arguments are"]



JuliaStochastic[c_, numpoints_, opts___] := Module[
	{inv1,inv2,points,z0,oneBias,chooser, 
	valid = First /@ Options[JuliaStochastic]},
	
	Scan[If[!MemberQ[valid, First[#]],
		Message[JuliaSet::optx, ToString[First[#]]]]&,
		{opts}
	];

	oneBias = OneBias /. {opts} /. Options[JuliaStochastic];
	If[!(0<oneBias && oneBias<1),
		{Message[OneBias::outOfRange];
		oneBias = OneBias /. Options[JuliaStochastic]}];
	
	inv1 = N[Sqrt[#-c]]&; inv2 = N[-Sqrt[#-c]]&;
	z0 = Nest[inv1[inv2[#]]&, 1, 5];
	points = ComposeList[Table[chooser = Random[];
		If[chooser<oneBias, inv1, inv2], 
		{numpoints}], z0];
	ListPlot[{Re[#],Im[#]}& /@ points,  
		fR[ListPlot,opts], 
		fR[ListPlot,Sequence @@ Options[JuliaStochastic]]
	]
] /;  
	( NumberQ[c] || Message[JuliaSet::badNumber, JuliaStochastic, c] ) &&
	  ( (IntegerQ[numpoints] && numpoints > 0 ) ||
	  	Message[JuliaSet::badInt, JuliaStochastic, numpoints] ) &&
	  ( And @@ Map[OptionQ, {opts}] ||
	  	Message[JuliaSet::badOpt,JuliaStochastic] )
	  	
(* Special cases for JuliaStochastic *)
JuliaStochastic[]   := Message[JuliaSet::argu, "JuliaStochastic", 
	"0 arguments", "2 arguments are"]
JuliaStochastic[c_] := Message[JuliaSet::argu, "JuliaStochastic", 
	"1 argument",  "2 arguments are"]

      
JuliaIterationsCom = Compile[{{z, _Complex}, {c, _Complex}, {iters, _Integer}}, 
   Module[{i, s}, i = 1; s = z; While[i < iters && Abs[s] < 2, i++; s = s^2 + c]; i]];

Options[FilledJuliaSet] = {Iterations->100, PlotPoints->100, JuliaSetOnlyQ->False, ColorQ->True};

FilledJuliaSet[c_, {x0_,x1_},{y0_,y1_}, opts___] :=  (
  {iters,pp, JSonlyQ, colQ} = {Iterations,PlotPoints,JuliaSetOnlyQ, ColorQ} /. {opts} /. Options[FilledJuliaSet];
  DensityPlot[
    JuliaIterationsCom[x + y I, c, iters], {x, x0, x1}, {y,y0,y1}, Evaluate[fR[DensityPlot, opts]],PlotPoints -> pp, PlotRange->All,AspectRatio->Automatic,
    ColorFunction -> If[colQ, (If[# == 1, Black, If[JSonlyQ,White,Hue[0.9 #^(3/4)]]] &),(If[# == 1, Black, GrayLevel[ #^0.25]] &)]])
      
      
End[];  (* End Private Context *)

(*Protect[JuliaSet, Julia, JuliaModified, JuliaSimple, 
	JuliaStochastic, Bound, Resolution, OneBias]*)

EndPackage[];

(* JuliaModified[-.83-.18I]*) Show[FilledJuliaSet[-.83-.18I, {-1.6,1.6}, {-.9,.9}, Frame->True, FrameTicks->{{{0,{.5,"0.5 I"},{-.5,"-0.5 I"}},{}},{{-1,0,1},{}}}, PlotRangePadding->0, Iterations->120, PlotPoints->120   ]]

Code for Figure : Feedback Screen Shot

pic=;

ImageDimensions[pic]

pic

Code for Figure : Simulation Example

scale = 2;
screen=;
{sm,sn}=ImageDimensions[screen];
screen=ImageResize[screen, scale{sm,sn}];
{sm,sn}=ImageDimensions[screen]

backSize=1000 scale;complexSize=2.4;
twoSquare=Polygon[{{{0,0},{1,0},{1,1},{0,1}},{{1,1},{2,1},{2,2},{1,2}}}];
back=Graphics[{RGBColor[1,.25,.25], Rectangle[{0,0}, {10,10}],RGBColor[.15,.15,1], Translate[twoSquare,Flatten[Table[{i,j},{i,0,8,2},{j,0,8,2}],1]]},PlotRangePadding->0,Frame->True,FrameTicks->False,FrameStyle->Thick,ImageSize->backSize];
im=Image[Rasterize[Overlay[{back,screen},Alignment->Center, ImageSize->backSize],RasterSize->backSize/scale]];
screenOnBackground=Show[im,Graphics[{Thickness[0.005],Lighter[Yellow],Circle[Scaled[{0.5,0.5}],Scaled[(1/(2 complexSize))]]}],ImageSize->200]

Quiet[LaunchKernels[]];

c=0;  complexSize=2.4;   p=0.;
limitColor = {0.15,0.95,0.95};
{xmin,xmax}={-complexSize,complexSize};
{ymin,ymax}={-complexSize,complexSize};
{m,n}=ImageDimensions[im];
imageData = ImageData[im];

color1 = With[{c=c,imageData=imageData,m=m,n=n,p=p,limitColor=limitColor,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax},Compile[{{z,_Complex}},Module[{z1},z1=(z^2)+c;If[!(xmin≤Re[z1]≤xmax&&ymin≤Im[z1]≤ymax),{0,0,0},p limitColor+(1-p) imageData[[Round[Max[1,((n (Im[z1]-ymax))/(ymin-ymax))]],Round[Max[1,((m (Re[z1]-xmin))/(xmax-xmin))]]]]]],RuntimeAttributes->{Listable},Parallelization->True]];
colors=color1[Table[x+I y,{y,ymax,ymin,((ymin-ymax)/n)},{x,xmin,xmax,((xmax-xmin)/m)}]];
im = Image[colors];
screenShows1=im=ImageTake[im,{((backSize-sn)/2),backSize-((backSize-sn)/2)},{((backSize-sm)/2),backSize-((backSize-sm)/2)}];
cameraSees1 = Show[{im=Image[Rasterize[Overlay[{back,im},Alignment->Center],RasterSize->(backSize/scale)]],Graphics[{Thickness[0.005],Lighter[Yellow],Circle[Scaled[{0.5,0.5}],Scaled[(1/(2 complexSize))]]}]},ImageSize->300]

cameraSees1

im=Image[Rasterize[Overlay[{back,screen},Alignment->Center], RasterSize->backSize/scale]];
Monitor[Do[imageData=ImageData[im];color1=With[{c=c,imageData=imageData,m=m,n=n,p=p,limitColor=limitColor,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax},Compile[{{z,_Complex}},Module[{z1},z1=(z^2)+c;If[!(xmin≤Re[z1]≤xmax&&ymin≤Im[z1]≤ymax),{0,0,0},p limitColor+(1-p) imageData[[Round[Max[1,((n (Im[z1]-ymax))/(ymin-ymax))]],Round[Max[1,((m (Re[z1]-xmin))/(xmax-xmin))]]]]]],RuntimeAttributes->{Listable},Parallelization->True]];colors=color1[Table[x+I y,{y,ymax,ymin,((ymin-ymax)/n)},{x,xmin,xmax,((xmax-xmin)/m)}]];im=Image[colors];im=ImageTake[im,{((backSize-sn)/2),backSize-((backSize-sn)/2)},{((backSize-sm)/2),backSize-((backSize-sm)/2)}];im=Image[Rasterize[Overlay[{back,im},Alignment->Center],RasterSize->(backSize/scale)]],{i,5}],i];

screenShowsFinally = Image[Show[{im,
Graphics[{Thickness[0.005],Lighter[Yellow],
Circle[Scaled[{0.5,0.5}],Scaled[1/(2*complexSize)]]}];{}}, ImageSize->backSize]];
ImageDimensions[screenShowsFinally]
Show[screenShowsFinally,ImageSize->400]

GraphicsColumn[{GraphicsRow[{Show[screenOnBackground,ImageSize->160],
Show[cameraSees1,ImageSize->160]},Spacings->20],
Show[ImageTake[screenShowsFinally,{(backSize-sn)/2+1,backSize-(backSize-sn)/2},
{(backSize-sm)/2+1,backSize-(backSize-sm)/2}],ImageSize->.744{460,460*3/4}]},
Spacings -> -32]

Code for Figure : Affine Julia

itCnt = Compile[{{a,_Real},{b,_Real},{s,_Real},
{x0,_Real},{y0,_Real}},
Module[{cnt,x,xtemp,y},
cnt = 0;
x=x0;y=y0;
While[cnt<1000&&(x^2)+(y^2)≤8,
cnt=cnt+1;
xtemp=((((x/s))^2)-(y^2)+a);
y=2(x/s)*y+b;
x=xtemp];
cnt]];
Clear[a,b,s];
f[a_,b_,s_,x_,y_] = {((((x/s))^2)-(y^2)+a),2(x/s)*y+b};
p0 = {x,y} /. First[
NSolve[f[-0.9,0.3,1.22,x,y]=={x,y},{x,y},Reals]];
Rasterize[ArrayPlot[data1=Table[itCnt[-0.9,0.3,1.22,x,y],
{y,1.25,-1.25,-0.005},{x,-2.15,2.15,0.005}],
ColorFunction->(If[#>.99, RGBColor[0.3,0.3,1],GrayLevel[ #^.15]]&),
DataRange -> {{-2.15,2.15},{-1.25,1.25}},PlotRange->{{-2.15,2.15}, {-1.25,1.25}},
PlotRangePadding->0,
FrameTicks -> {{{-1,0,1}, None},{{-2,-1,0,1,2}, None}}, Epilog -> {Yellow,PointSize[.01],Point[p0],
Line[{{-2,-2},{2,2}}];{}}],
RasterSize->1000,ImageSize->380]

Clear[f,a,b,s];
f[a_,b_,s_,x_,y_] = {((((x/s))^2)-(y^2)+a),2(x/s)*y+b};
a=-1;b=0.05; s=1.3;
interestingPoints = {x,y}/.
NSolve[f[a,b,s,Sequence@@f[a,b,s,x,y]]=={x,y},{x,y},Reals];//Quiet
p0=interestingPoints[[3]];
{p1,p2}=interestingPoints[[{2,4}]];
itCnt = Compile[{{a,_Real},{b,_Real},{s,_Real},
{x0,_Real},{y0,_Real}},
Module[{cnt,x,xtemp,y},
cnt = 0;
x=x0;y=y0;
While[cnt<500&&(x^2)+(y^2)≤8,
cnt=cnt+1;
xtemp=((((x/s))^2)-(y^2)+a);
y=2(x/s)*y+b;
x=xtemp];
cnt]];
Rasterize[ArrayPlot[data2=Table[itCnt[-1,0.05,1.3,x,y],
{y,1.25,-1.25,-0.005},{x,-2.45,2.45,0.005}],
DataRange -> {{-2.45,2.45 },{-1.25,1.25}},ColorFunction->(If[#>(.35),RGBColor[0.3,0.3,1];Yellow,GrayLevel[#^.3]]&),
FrameTicks -> {{{-1,0,1}, None},{{-1,0,1}, None}}, Epilog -> {{Yellow;RGBColor[0.3,0.3,1],PointSize[.01],Point[p0]},
{Red,PointSize[.008],Point[{p1,p2}]}}, PlotRangePadding->0],
RasterSize->1000,ImageSize->380]

Union[data1f=Flatten[data1]]
{Count[data1f, 100], Length[data1f]}

Union[data2f=Flatten[data2]]
{Count[data2f, 165], Length[data2f]}

Code for Cover Figure: Eiffel Tower

pic=Import[FileNameJoin[{NotebookDirectory[],"EiffelTowerOriginal.jpg"}]];
Show[pic, ImageSize->200]
ImageDimensions[pic]
pic=ImageTake[pic,{40, 2217}, {0, 1440}] ;
ImageDimensions[pic]

f[z_] := (z^2);
Clear[n];

ImageDimensions[pic]
invImage[im_Image] := Image[cInvImage[ImageData[im]],ImageSize -> 150];

{xmin,xmax}={-2,2}; {ymin,ymax}={-4.,2};
DynamicModule[{a=0,b=0,da=2178,db=1440},pic1=ImageTake[pic,{a,da},{b,db}];{m,n}=ImageDimensions[pic1];si[x_]=((x (xmax-xmin))/m)+xmin;sx[i_]=x/.First[Solve[si[x]==i,x]];sj[y_]=((y (ymax-ymin))/n)+ymin;sy[j_]=y/.First[Solve[sj[y]==j,y]];fScaled=ComplexExpand[{Re[f[x+I y]],Im[f[x+I y]]}]/.{x->si[x],y->sj[y]};fScaledBack={sx[fScaled[[1]]],sy[fScaled[[2]]]};intInv=With[{fScaledBack=fScaledBack},Compile[{{x,_Integer},{y,_Integer}},Round[fScaledBack]]];rgb2=With[{intInv=intInv,m=m,n=n},Compile[{{x,_Integer},{y,_Integer},{imageData,_Real,3}},Module[{x2,y2},{x2,y2}=intInv[x,y];If[1≤x2<m&&1≤y2<n,imageData[[y2,x2]],{0.`,0.`,0.`}]]]];cInvImage=With[{rgb2=rgb2,m=m,n=n},Compile[{{imageData,_Real,3}},Table[rgb2[x,y,imageData],{y,1,n},{x,1,m}]]]];
final=invImage[pic1];

Show[ImageTake[final, {750, 2180},{30,1420}], ImageSize->600]