filter droste_effect (image in, bool No_Transparency, bool External_Transparency, bool Mirror_Effect, bool untwist, bool Do_Not_Flatten_Transparency)
##
# Droste Effect Code
# by breic ( www.flickr.com/photos/breic )
# Modified for Windows and updated by Pisco Bandito (www.flickr.com/photos/joshsommers)
# see also the flickr group: www.flickr.com/groups/escherdroste
# Version 7 - For Windows and Linux
# last modified: 4/20/2007
##

r1=1; #Inner Radius (must be less than r2)
r2=1; #Outer Radius

p1=1; #Periodicity
p2=1; #Number of strands

# procedural scaling and rotation (Whole numbers)
zoom=1;
rotate=0;

#Procedural post-effect shifting (Whole numbers)
xShift=0;
yShift=0;

#Center Shift
#Use these to calibrate an offcenter transparent region. Altering these values will move the center of the spiral.
#Note that you will need to also use the xShift and yShift parameters in conjunction with these to get your image framed correctly.

xCenterShift = 0;
yCenterShift = 0;

#Level Control
StartingLevel=1;
NumberOfLevels=10;
LevelFrequency=1; #Frequency of levels. If 1, then every level. If 2 then every other level. If 3 then every third level, and so on.

#For the mirror effect
#Fiddle with this value depending on the values you have for p1 and p2 to get an alternating reflected effect
MirrorCoefficient=1;

##############################################
# No need to change anything below this line #
##############################################

# Miscellaneous useful variables
true=1;
false=0;
epsilon=.01; # used to avoid hard comparisons
imageX=W;
imageY=H;
minDimension=min(imageX, imageY);
rotate=pi/180*rotate;

if zoom > 0 then
zoom=(zoom*p1)/10;
end;

#User Values Defined
tileBasedOnTransparency=!(No_Transparency);
transparentPointsIn=!(External_Transparency);#Default is Internal
mirrorEffect = Mirror_Effect;
retwist=!(untwist);
flattenTransparency = !(Do_Not_Flatten_Transparency);

##
# 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
angle = 2*pi*p1;

if mirrorEffect then
angle=angle/MirrorCoefficient;
end;

if (0>p2) then
angle = -angle;
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];
z = z - (ri:[xShift,yShift] / 10);

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;

z=p1*z2/beta;
z=r1*exp(z);

## End Droste effect math

if (tileBasedOnTransparency && StartingLevel > 0) then
if (!transparentPointsIn) then
ratio=r2/r1*exp( I*angle);
end;
if ( transparentPointsIn) then
ratio=r1/r2*exp(-I*angle);
end;
z=z*exp(StartingLevel*log(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]+xCenterShift);
iy=minDimension/2*(z[1]+yCenterShift);
IXY = xy:[ix,iy];
ReturnPixel = in(IXY);
colorSoFar = colorSoFar + (ReturnPixel*(alpha(ReturnPixel)*alphaRemaining));
alphaRemaining=alphaRemaining*(1-alpha(ReturnPixel));
sign=0;

# do we need to look inward from the current point, or outward?

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 (r1 > radius) then
sign=-1;
end;
if (radius > r2) then
sign= 1;
end;
end;

if (0 > sign) then ratio=r2/r1*exp( I*angle); end;
if (sign > 0) then ratio=r1/r2*exp(-I*angle); end;

if (LevelFrequency > 1) then
ratio = exp(log(ratio)*LevelFrequency);
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=StartingLevel;
maxiteration=NumberOfLevels+StartingLevel-1;

while (sign != 0 && maxiteration > iteration) do
z2=z*ratio;
z=z2;
ix=minDimension/2*(z[0]+xCenterShift);
iy=minDimension/2*(z[1]+yCenterShift);
IXY = xy:[ix,iy];
sign=0;
if (tileBasedOnTransparency) then
ReturnPixel = in(IXY);
colorSoFar = colorSoFar + (ReturnPixel*(alpha(ReturnPixel)*alphaRemaining));
alphaRemaining=alphaRemaining*(1-alpha(ReturnPixel));
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]);
colorSoFar = in(IXY);
if (r1 > radius) then sign=-1; end;
if (radius > r2) then sign= 1; end;
end;
iteration=iteration+1;
end;

ReturnPixel=colorSoFar;

if (flattenTransparency) then
ReturnPixel = rgba:[red(ReturnPixel),green(ReturnPixel), blue(ReturnPixel),1];
end;

#The return value for the current pixel is held in the variable ReturnPixel
ReturnPixel

end