########################################
## Droste Effect code for Mathmap     ##
## Original Code by Breic (Ben)       ##
## Adapted by Pisco Bandito (Josh)    ##
## Version 4.0                        ##
## This version for Windows and Linux ##
########################################

#You may need to alter the values of the following 9 variables to suit your image.

r1=.4; # r1 is the inner radius. r1 is greater than zero and less than r2
r2=1; # r2 is the outer radius. r2 is greater than r1 and less than 1
rinner=user_float("inner radius", 0,1); # to experiment with r1 and r2, uncomment these lines
router=user_float("outer radius",0,1);
r1 = min(rinner, router); # sanity check
r2 = max(rinner, router);

p2=user_int("p2 (# strands)",-1,4); # number of strands - the number of "arms" the spiral will have
p1=user_int("p2 (# strands)",-1,4); # number of strands; # periodicity - the number of times the image will repeat per cycle
setp1=user_bool("auto-set p1");
if (setp1) then
p1= abs(p2)/2 * (1+sqrt( 1 - (log(r2/r1)/pi)^2) );
else
p1=user_int("p1 (periodicity)",-1,3); # periodicity
if (p1 == 0) then p1 = 1; end;
end;
# procedural scaling and rotation
zoom=user_float("zoom",1,10);
rotate=pi/180*user_float("rotate",0,360);

# Procedural Shifting
xShift=.2; #Between -1 and 1
yShift=0;  #Between -1 and 1

##
# 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;

############################################################
############################################################
## You should not have to change anything below this line ##
############################################################
############################################################

imageX=W; # image size, in pixels
imageY=H;
minDimension=min(imageX, imageY);

## User Variables, set these on the User Settings Tab ##
retwist=user_bool("Do Not Retwist (Leave Unchecked for Droste Effect)");
retwist=!retwist;

##
# 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=user_bool("Tile Based on Transparency?");
transparentPointsIn=user_bool("Transparency Points In?");

# Miscellaneous variables
true=1;
false=0;
epsilon=.01;

##
# 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)+xShift,(ybounds[0]+(ybounds[1]-ybounds[0])*(y+H/2)/H)+yShift];

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
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));
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*(alpha(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*(alpha(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)

#This last line is important, it returns the pixel value for the current pixel
color