#### By: breic

Here's some code which uses alpha values properly for blending. It also largely fixes the framing bug with a `levelsToLookOut` parameter. Higher values slow the code down, but allow for parts of images from multiple annular layers beyond the outside of the frame to contribute to the image — useful for more complicated alpha masks.

There is still a minor framing "bug": if you have a nonzero alpha value at the center [or attractive point] of the image, this will cause a discontinuity across the -x axis.

I haven't tested this code carefully, and I also removed some minor functionality (alpha values can only point inward *or* outward on any given image), so I'm leaving the previous version up for now. I've also added a few comments. Your comments are always welcome, you can leave them on the bottom of the page. :)

Thanks, Seb, for pointing out that vector-multiplying trick. This would be impossible without it. ;)

One missing feature which might be useful would be to have it output just one annular layer at a time, so you can stack them yourself in Photoshop.

*Edit: Previous code would not work when* `tileBasedOnTransparency=false`. *It should probably work now. Added comments.*

```
##
# Droste effect code
# by breic ( www.flickr.com/photos/breic )
# see also the flickr group: www.flickr.com/groups/escherdroste
# last modified: 1/11/2007
##
# Miscellaneous useful variables
true=1;
false=0;
epsilon=.01; # used to avoid hard comparisons
##
# Input parameters
##
# input should be a square image, *at 72dpi resolution*
imageX=1024; # image size, in pixels
imageY=1024;
minDimension=min(imageX, imageY);
# 0 < r1 < r2 < 1 determines the inner and outer fractional radii of the annulus which will be wrapped
#r1=user_float("inner radius", 0,1);
#r2=user_float("outer radius",0,1);
r1=.115;
r2=.7;
##
# See gadl's chessboard set (http://flickr.com/photos/gadl/sets/72157594305663827/) for an
# explanation of the (p1, p2) parameters -- although he may have switched them!
# Mathematically, we rotate point (p1,p2) in the annular lattice in log coordinates to (0,1).
# The labels at escherdroste.math.leidenuniv.nl/index.php?menu=im&sub...
# agree with our notation here.
# There is a good explanation of the Droste-effect math at
# www.josleys.com/articles/printgallery.htm .
##
p1=1; # periodicity
p2=1; # number of strands
# procedural scaling and rotation
zoom=user_float("zoom",0.1,10);
rotate=pi/180*user_float("rotate",0,360);
retwist=user_bool("retwist");
##
# Droste-effect code starts here
# Set Droste effect parameters
##
alpha=atan(p2/p1*log(r2/r1)/(2*pi));
f=cos(alpha);
beta=f*exp(I*alpha);
# the angle of rotation between adjacent annular levels
if (p2 > 0)
then angle = 2*pi*p1;
else angle =-2*pi*p1;
end;
##
# Code to set up the viewport properly
##
if (retwist) then
xbounds=[-r2,r2];
ybounds=[-r2,r2];
else
ybounds=[0,2.1*pi];
xbounds=[-log(r2/r1), log(r2/r1)];
end;
xymiddle=ri:[0.5*(xbounds[0]+xbounds[1]),0.5*(ybounds[0]+ybounds[1])];
xyrange=xy:[xbounds[1]-xbounds[0], ybounds[1]-ybounds[0]];
aspectRatio=W/H;
xyrange[0]=xyrange[1]*aspectRatio;
xbounds=[xymiddle[0]-0.5*xyrange[0],xymiddle[0]+0.5*xyrange[0]];
z=ri:[xbounds[0]+(xbounds[1]-xbounds[0])*(x+W/2)/W,ybounds[0]+(ybounds[1]-ybounds[0])*(y+H/2)/H];
if (retwist) then # only allow for procedural zooming/scaling in the standard coordinates
zinitial=z;
z=xymiddle+(z-xymiddle)/zoom*exp(-I*rotate);
else
zinitial=r1*exp(z); # save these coordinates for drawing a frame later
zinitial=zinitial*zoom*exp(I*rotate);
end;
##
# The Droste effect math all takes place over the next six lines.
# All the rest of the code is for niceties.
##
if (retwist) then z2=log(z/r1);
else z2 = z; end;
logz=z2; # save these coordinates for drawing a grid later
z=p1*z2/beta;
rotatedscaledlogz=z; # save these coordinates for drawing a grid later
z=r1*exp(z);
## End Droste effect math
##
# Tiling can be based on transparency (if the input image is a tiff), or simply based on the
# radius. Using transparency, there can be protrusions between different annular layers.
# Tiling based on transparency, you can decide whether you want to look inward or
# outward from a transparent pixel. For example, with a frame you'll want to look inward,
# while for a flower you'll want to look outward.
##
tileBasedOnTransparency=false;
transparentPointsIn=true;
##
# To avoid framing problems on the largest annulus when tiling based on transparency, look
# outside (levelsToLookOut) levels to see if something farther out should cover up this pixel
# Try setting to 0 to see framing errors; 1 should be sufficient unless you have three or more
# image layers contributing to some pixel (in which case set it to 2 or more). Larger values
# slow the code down, and may lead to floating point errors.
##
levelsToLookOut=3;
if (tileBasedOnTransparency && levelsToLookOut > 0) then
if ( transparentPointsIn) then ratio=r1/r2*exp(-I*angle); end;
if (!transparentPointsIn) then ratio=r2/r1*exp( I*angle); end;
z=z*exp(levelsToLookOut*log(ratio));
#z=z*ratio;
end;
##
# When tiling based on transparency, color is accumulated into the colorSoFar variable,
# while alphaRemaining tells how much remains for lower layers to contribute (initially 1,
# finally 0).
##
colorSoFar=rgba:[0,0,0,0];
alphaRemaining=1;
ix=minDimension/2*z[0];
iy=minDimension/2*z[1];
color=origValXY(ix,iy);
colorSoFar=colorSoFar+color*alphaRemaining;
alphaRemaining=alphaRemaining*(1-alpha(color));
# do we need to look inward from the current point, or outward?
sign=0;
if (tileBasedOnTransparency) then
if ( transparentPointsIn && alphaRemaining > epsilon) then sign=-1; end;
if (!transparentPointsIn && alphaRemaining > epsilon) then sign= 1; end;
else
radius=sqrt(z[0]*z[0]+z[1]*z[1]);
if (radius < r1) then sign=-1; end;
if (radius > r2) then sign= 1; end;
end;
if (sign < 0) then ratio=r2/r1*exp( I*angle); end;
if (sign > 0) then ratio=r1/r2*exp(-I*angle); end;
##
# Iteratively move inward or outward, until
# the point has radius r in [r1, r2), if tileBasedOnTransparency=false
# or until alphaRemaining=0, if tileBasedOnTransparency=true
# In the latter case, we accumulate color at each step
##
iteration=0; maxiteration=10;
while (sign != 0 && iteration < maxiteration) do
z2=z*ratio;
z=z2;
rotatedscaledlogz=rotatedscaledlogz+ri:[0,-sign*angle];
ix=minDimension/2*(z[0]);
iy=minDimension/2*(z[1]);
color=origValXY(ix,iy);
colorSoFar=colorSoFar+color*alphaRemaining;
alphaRemaining=alphaRemaining*(1-alpha(color));
radius=sqrt(z[0]*z[0]+z[1]*z[1]);
sign=0;
if (tileBasedOnTransparency) then
if ( transparentPointsIn && alphaRemaining > epsilon) then sign=-1; end;
if (!transparentPointsIn && alphaRemaining > epsilon) then sign= 1; end;
else
radius=sqrt(z[0]*z[0]+z[1]*z[1]);
if (radius < r1) then sign=-1; end;
if (radius > r2) then sign= 1; end;
end;
iteration=iteration+1;
end;
color=colorSoFar;
color=rgba:[color[0], color[1], color[2], 1]; # set the alpha value to 1 (it could be <1 if the loop terminated at iteration maxiteration)
##
# The green grid shows each basic annulus, with (outer radius)/(inner radius) = r2/r1.
# The blue grid shows the green grid mapped into log coordinates (in which it truly is a grid),
# rotated, and exponentiated forward again. This allows visualization of the Droste
# transformation.
##
showGrid=user_bool("show grid"); # draw green and blue grids showing the boundaries of the wrapped annulus
if (showGrid) then
gridz=xy:[(logz[0]+10*log(r2/r1))%log(r2/r1), (logz[1]+10*2*pi)%(2*pi)];
if (gridz[0] < epsilon || gridz[0] > (log(r2/r1)-epsilon) || gridz[1] < epsilon || gridz[1] > (2*pi-epsilon)) then color=rgba:[0,1,0,1]; end;
gridz=xy:[(rotatedscaledlogz[0]+10*log(r2/r1))%log(r2/r1), (rotatedscaledlogz[1]+10*2*pi)%(2*pi)];
if (gridz[0] < epsilon || gridz[0] > (log(r2/r1)-epsilon) || gridz[1] < epsilon || gridz[1] > (2*pi-epsilon)) then color=rgba:[0,0,1,1]; end;
end;
##
# The frame is a black-and-white border in the standard view, and its preimage in the log
# coordinates. In log coordinates, the shaded region to the right of the frame will be mapped
# outside the viewport, while the region to the left will be visible.
##
showFrame=user_bool("show frame"); # draw red frame around the boundary of the image
if (showFrame) then
gridz=xy:[zinitial[0],zinitial[1]];
if (gridz[0] < (aspectRatio*r2) && gridz[0] > -(aspectRatio*r2) && gridz[1] < r2 && gridz[1] > -r2) then
dx=min((aspectRatio*r2)-gridz[0], gridz[0]+(aspectRatio*r2));
dy=min(r2-gridz[1], gridz[1]+r2);
if (dx < (4*epsilon) || dy < (4*epsilon)) then
color=rgba:[1,1,1,1];
end;
if (dx < (2*epsilon) || dy < (2*epsilon)) then
color=rgba:[0,0,0,1];
end;
else
color=rgba:[0.75*red(color),0.75*green(color),0.75*blue(color),1];
color=color;
end;
end;
color
```