The quadratic camera

ImageTake[Import[FileNameJoin[{NotebookDirectory[],"quadcam.tiff"}],
ImageSize -> 500, ImageResolution->800],{20,-40},{40,-30}]

On my Macbook Pro, I've got a Photo Booth like plug-in that transforms the image seen by the built in camera into the inverse of the image under the complex square function.  If you use an attached monitor, an external camera, or even just a mirror to create a video feedback loop, you can use this to generate Julia set like images.  The picture above shows the program in action as I point the built in camera at my bathroom mirror.  This post shows the program in action and provides an implementation in Mathematica version 8.  It also serves as a computational supplement to a paper I co-authored that's been submitted to The Mathematical Intelligencer.

Iterated inverses of an image

Here's the screen of my MacBook. I simply used Mac's Grab to generate a TIFF file and dragged the result into Mathematica.  The result is an Image with dimensions sm and sn.

Now, suppose that we attach a camera to the computer and point it at the screen.  Generally, the field of vision of the camera will be larger than just the screen.  In this simulation, we generate an artificial, checkered background representing what the computer sees behind the screen.  Each pixel in this square image is associated with a complex number; the yellow circle is the unit circle in the complex plane.

backSize=1000;
twoSquare = 
Polygon[{
{{0,0},{1,0},{1,1},{0,1}},
{{1,1},{2,1},{2,2},{1,2}}}];
back=Graphics[{Darker[Gray,0.4],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]]];
screenOnBackground = Show[im,
Graphics[{Thickness[0.005],Lighter@Yellow,
Circle[Scaled[{0.5,0.5}],Scaled[1/(2*2.4)]]}],
ImageSize->400]

The image on the screen is now passed to a digital filter.  For each pixel, we square the corresponding complex number and then color the original point according to where it lands in the field of view of the camera.  We have essentially computed the inverse image of the camera's view under the complex square function.  We then replace the screen in the previous image with that inverse image.  Experience shows that there is also a damping effect that biases the color towards a limiting cyan color.  Here is the result, taking all of this into account.

c=0; (* The circle *)
complexSize=2.4;
p=0.2;
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}]];
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]]],
Graphics[{Thickness[0.005],Lighter[Yellow],
Circle[Scaled[{0.5,0.5}],Scaled[1/(2*2.4)]]}]},
ImageSize->400]

If we iterate this procedure 15 times, the screen will show the following.

im=Image[Rasterize[Overlay[{back,screen},
Alignment->Center]]];
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}]];
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]]],{i,15}],i];
screenShowsFinally = ImageTake[Image[Show[{im,
Graphics[{Thickness[0.005],Lighter[Yellow],
Circle[Scaled[{0.5,0.5}],Scaled[1/(2*2.4)]]}]}]],
{(backSize-sn)/2,backSize-(backSize-sn)/2},
{(backSize-sm)/2,backSize-(backSize-sm)/2}]

It's well known that inverse iteration of quadratic functions produces Julia sets.  In this case, we've generated the Julia set of f(z)=(z^2).  Now, suppose that we don't the camera directly at the center of the screen, rather we translate the camera so that we are pointed directly at complex number c.  Then, this procedure effectively computes the inverse image under the quadratic function f(z)=(z^2)+c.  Thus, we might expect to see a more interesting Julia set.  Here's 32 iterations for the case where c=-0.83-0.18i.

{sm,sn}=ImageDimensions[screen];
backSize=1000;
twoSquare = 
Polygon[{
{{0,0},{1,0},{1,1},{0,1}},
{{1,1},{2,1},{2,2},{1,2}}}];
back=Graphics[{Darker[Gray,0.4],
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]]];
imageData = ImageData[im];
c=-0.83-0.18 I ;
complexSize=2.4;
p=0.125;
limitColor = {0.15,0.95,0.95};
{xmin,xmax}={-complexSize,complexSize};
{ymin,ymax}={-complexSize,complexSize};
{m,n}=ImageDimensions[im];
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}]];
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]]],{i,32}],i];
screenShowsFinally = ImageTake[im,
{(backSize-sn)/2,backSize-(backSize-sn)/2},
{(backSize-sm)/2,backSize-(backSize-sm)/2}]

The quadratic camera with Mathematica

The previous code all works with a single static image for explanatory purposes.  Mathematica V8 allows easy access to a camera attached to the computer so this can all be done dynamically.  First, we fire up the camera and determine the resolution of the corresponding image.

pic = CurrentImage[];
{m,n}=ImageDimensions[pic]

Now, we define the function we want to play with.  I suppose the classic case is f(z)=(z^2) but this should work for any polynomial.  We also define some code to translate between the complex plane and the space of the image.

f[z_] = (z^2);
{xmin,xmax}={-2,2};
{ymin,ymax}={-1.5,1.5};
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}]]];
invImage[im_Image] := Image[cInvImage[
ImageData[im]],ImageSize -> 1.5{m,n}];

Finally, we do it!

(* Do it with dynamic image! *)
Dynamic[invImage[CurrentImage[]]]

OK, point the camera at the screen or point your laptop at a mirror.  You can use the following command to shut down the camera.

IMAQ`StopCamera[]