# Fractals/Iterations in the complex plane/slope

It is an algorithm for converting 2D image to 3D

## Description

"three points are iterated for each pixel and then are used to create a three-dimensional surface. A vector perpendicular to that surface is used as the z variable and is passed to the coloring routine. When a Slope formula is used with the Lighting coloring, the result resembles a fractal surface being illuminated from above and to the side." Kerry Mitchell[1]

## names

Differences

• Slope = Identifies the slope (gradient or steepness) from each cell of a raster.

### slope of the line

In mathematics, the slope or gradient of a line is a number that describes both the direction and the steepness of the line.[3] Slope is often denoted by the letter m; there is no clear answer to the question why the letter m is used for slope,

y = mx + b


Slope of the line through 2 points: (x1,y1) and (x2,y2) is:

${\displaystyle m={\frac {y_{2}-y_{1}}{x_{2}-x_{1}}}.}$

the slope m of a line is related to its angle of inclination θ by the tangent function:

${\displaystyle m=\tan(\theta )}$

## Apply to

• Mandelbrot[4]
• Julia
• Newton [5]
• Kleinian group fractals[6]
• strange attractors [7]

## code

• Ultrafractal formulas by Damien Jones in dmj.ufm and standard.ufm

### c

hue = fmod(1 + carg(de) / (2 * pi), 1); //  slope of complex number de


### matlab

Matlab code by nickspiker

%for each pixel do the stuff below
c=somethingreal+somethingimaginary;
z=0;
d=1;
while (real(z)^2+imaginary(z)^2)<realbignumber
z=z^2+c;
d=2*z*d+1;%derivative step and I think the plus one comes from the derivative of c, not totally sure on that one
endwhile
slope=sign(z/d)%values are strange because of where we abruptly end the step but the sign is clean

%then this is what I do for the shading bit
%you can look at the real and imaginary parts of slope before you do the steps below to better understand what's going on
sunangle=45;
brightness=slope*e^(i*sunangle*pi/180);%this rotates the sign by the angle above
brightness=real(brightness)/2+.5;


%Mandelbrot rendering algorithms utilizing interior and exterior slope finding steps copyright 2020 Nick Spiker nick@nickspiker.com

maxn=2^32;
steps=2^5;%anti aliasing steps
mns=2^8;%max newton steps
dpi=90;
bleed=0;%change to 6 when final render, 3 per side
width=42;
height=48;
zoom=3;
center=-.75;
rotation=90;
alias=2^-6;
tiles=3;
starttile=1;
preview=false;

if preview
dpi=30;
tiles=1;
bleed=0;
end

rez=[(width+bleed).*dpi,(height+bleed).*dpi];
zoom=zoom./height.*(height+bleed);
localrez=ceil(rez./tiles);%finds closest multiple equal to or greater than rez
localrezb=localrez+2;
pixelsize=zoom./rez(2);
tlx=-zoom./rez(2).*rez(1)./2;
tly=zoom./2;
localzoom=localrezb(2).*pixelsize;

tic;
for tiley=floor((starttile-1)./tiles)+1:tiles
for tilex=mod(starttile-1,tiles)+1:tiles
tn=toc;
starttile=1;
localcenter=tlx+pixelsize.*localrez(1).*(tilex-.5)+(tly-pixelsize.*localrez(2).*(tiley-.5)).*i;
localcenter=center+localcenter.*exp(1i*rotation*pi/180);

xlim = [(pixelsize-localzoom./localrezb(2).*localrezb(1))./2,(localzoom./localrezb(2).*localrezb(1)-pixelsize)./2];
ylim = [(localzoom-pixelsize)./2,(pixelsize-localzoom)./2];

x=linspace(xlim(1),xlim(2),localrezb(1));
y=linspace(ylim(1),ylim(2),localrezb(2));
[xGrid,yGrid]=meshgrid(x,y);

cleanc=xGrid+1i*yGrid;
cleanc=cleanc*exp(1i*rotation*pi/180)+localcenter;

stackalpha=zeros(localrezb(2),localrezb(1),3);
stackimage=zeros(localrezb(2),localrezb(1),3);
stackalias=true(localrezb(2),localrezb(1));
aliasstep=alias;

for step=1 :steps
tm=toc;
aliasstep=aliasstep+alias;
totalpixels=sum(sum(stackalias));
if totalpixels==0
break
end
c=(cleanc(stackalias)+rand(totalpixels,1).*(cleanc(1,1)-cleanc(1,2))+rand(totalpixels,1).*(cleanc(1,1)-cleanc(2,1)));
ac=abs(c);
acs=ac.*ac;
acs=acs.*acs;
acs=acs.*acs;
cr=real(c);
ci=imag(c);
inside=false(totalpixels,1);
a=zeros(totalpixels,1);
frame=cat(3,a,a,a);
dr=a;
di=a;
er=a;
ei=a;
slope=a;
maxcount=1;
parfor p=1:totalpixels%parfor
insidep=false;
%orbitcount=1;
crp=cr(p);
cip=ci(p);
zr=crp;
zi=cip;
ap=maxn;
drp=0;
dip=0;
dsp=0;
erp=0;
eip=0;
esp=0;
srp=1;
sip=0;
osr=1;
osi=0;
%orbit=0;
zsmallest=realmax;
for n=1:maxn

srpt=2.*(srp.*zr-sip.*zi);
sip=2.*(srp.*zi+sip.*zr);
srp=srpt+1;

zrs=zr.*zr;
zis=zi.*zi;
zs=zrs+zis;
ns=n.*n;
if zs>4
llzs=log(log(zs));
taper=(llzs-0.326634259978281)./6.238324625039508;
taper=1-taper.*taper;
taper=taper.*taper;
taper=taper.*taper.*taper.*ns;
zsq=taper./sqrt(zs);
esp=esp+taper;
erp=erp+zr.*zsq;
eip=eip+zi.*zsq;

zsi=ns./zs;
dsp=dsp+zsi;
drp=drp+zr.*zsi;
dip=dip+zi.*zsi;

if zs>1.340780792994260e+154
ap=n-(llzs-5.174750624576761).*1.442695040888963;
break
end

durt=2.*(zr.*osr-zi.*osi)+1;
osi=2.*(osr.*zi+osi.*zr);
osr=durt;

zi=2.*zr.*zi+cip;
zr=zrs-zis+crp;
else
taper=ns;
zsq=taper./sqrt(zs);
esp=esp+taper;
erp=erp+zr.*zsq;
eip=eip+zi.*zsq;
zsi=ns./zs;
dsp=dsp+zsi;
drp=drp+zr.*zsi;
dip=dip+zi.*zsi;
if zsmallest>zs
zsmallest=zs;
period=n;
wr=zr;
wi=zi;
for l=1:mns
ur=wr;
ui=wi;
dur=1;
dui=0;
for m=1:period

%du=2*du*u
durt=2.*(dur.*ur-dui.*ui);
dui=2.*(dur.*ui+dui.*ur);
dur=durt;

%u=u*u+c
uis=ui.*ui;
ui=2.*ur.*ui+cip;
ur=ur.*ur-uis+crp;

end

%nw=w-(u-w)/(du-1)
nwr=ur-wr;
nwi=ui-wi;
duro=dur-1;
durt=duro.*duro+dui.*dui;
nwrt=(nwr.*duro+nwi.*dui)./durt;
nwi=(nwi.*duro-nwr.*dui)./durt;
nwr=nwrt;
ns=nwr.*nwr+nwi.*nwi;
wr=wr-nwr;
wi=wi-nwi;

if ns<2^-32
if dur.*dur+dui.*dui<1
drp=dur;
dip=dui;
erp=wr;
eip=wi;
insidep=true;
ap=m;
break
end
end
end
if insidep
break
end
end
end

durt=2.*(zr.*osr-zi.*osi)+1;
osi=2.*(osr.*zi+osi.*zr);
osr=durt;

zi=2.*zr.*zi+cip;
zr=zrs-zis+crp;

end
if insidep
break
end
end

maxcount=max(maxcount,n);
if insidep
inside(p)=true;
dr(p)=drp;
di(p)=dip;
er(p)=erp;
ei(p)=eip;

cca=crp+cip.*i;
ccc=cca;
for l=1:mns
zzz=ccc;
ddd=1;
for z=2:period
ddd=2.*ddd.*zzz+1;
zzz=zzz.*zzz+ccc;
end
ddd=ccc-zzz./ddd+2^-256.*i;
if ddd==ccc
ccc=ddd;
break
end
ccc=ddd;
end
slope(p)=ccc-cca;
else
dsp=1./dsp;
esp=1./esp;
dr(p)=drp.*dsp;
di(p)=dip.*dsp;
er(p)=erp.*esp;
ei(p)=eip.*esp;
slope(p)=sign((osr+osi.*i)./(zr+zi.*i));
end
a(p)=ap;
end

tm=toc-tm;
hrs=floor(tm/(60*60));
mnt=floor((tm-hrs*60*60)/60);
snd=tm-hrs*60*60-mnt*60;
['Tile ',num2str(tilex+tiley.*tiles-tiles),' of ',num2str(tiles.^2),', ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex),', Step ',num2str(step),', ',num2str(hrs),' hours, ',num2str(mnt),' minutes, ',num2str(snd),' seconds per step, Max N=2^',num2str(log2(maxcount)),' alias%',num2str(totalpixels./(localrezb(1).*localrezb(2)).*100)]

ia=a(inside);%interior period count
oa=a(not(inside));%exterior iteration count

d=dr+di.*i;
ic=d(inside);%interior coordinate,
ea=d(not(inside));%exterior average

e=er+ei.*i;
iw=e(inside);%interior W
ef=e(not(inside));%Exterior flames

is=slope(inside);

inside3=cat(3,inside,inside,inside);
stackalias3=cat(3,stackalias,stackalias,stackalias);

%start color
interiorrotationoffset=1.5;
exteriorrotationoffset=pi./4;
squarerotation=-.05;
squarebrightness=.7;
fresnelbrightness=.125;
reflectionrotation=.1;
reflectionoffset=1/4;
reflectioncurve=1/16;
reflectionwavies=1/8;
reflectionbrightness=0.125;
slopescale=.5;

aic=abs(ic);
isreflect=sign(is).*aic;
iss=isreflect.*exp(i.*squarerotation);
cc=exp(i.*(angle(iw)+1./aic+interiorrotationoffset)).*aic;
waves=real(exp(i.*(angle(is)+1./aic)).*aic);
o=real(cc.*exp(i*120*pi/180));
p=real(cc);
q=real(cc.*exp(i*270*pi/180));
cc=cat(3,o,p,q);
cc=cc./2+.5;
aic=aic.*aic;
iss=iss.*aic;
aic=aic.*aic;
aics=aic.*aic;
aic=aics+fresnelbrightness;
aics=aics.*aics;

issr=real(iss);
issi=imag(iss);
cc=1-cc.*(1-(and(and(issr>.3,issr<.7),and(issi>0.1,issi<.4)).*squarebrightness+(real(isreflect.*exp(reflectionrotation.*i))>(reflectionoffset-aics.*reflectioncurve+waves.*reflectionwavies)).*reflectionbrightness).*aic);

frame(inside3)=cc;

outsideslope=slope(not(inside));
outsideslope=0.5-real(outsideslope.*exp(exteriorrotationoffset.*i))./2;

oc=erf(12./oa);
%cc=abs(ef).*exp(i.*(angle(ef)+oa.*16.*oa));
cc=ef;
o=real(cc.*exp(i*120*pi/180));
p=real(cc);
q=real(cc.*exp(i*270*pi/180));
cc=cat(3,o,p,q);
cc=cc./2+.5;
oc3=5./oa;
oc3=1-oc3;
oc3=1-cat(3,oc3,oc3,oc3).*(1-cc)./2;
cc=oc3.*oc.*(outsideslope.*slopescale+(1-slopescale));
frame(not(inside3))=cc;

if preview
stackimage(stackalias3)=reshape(frame,size(frame,1).*size(frame,2).*size(frame,3),1);
imshow(stackimage)
return
%else
%frame=frame./2+.25;%declip for post
end
%end color

stackimage(stackalias3)=stackimage(stackalias3)+reshape(frame,size(frame,1).*size(frame,2).*size(frame,3),1);
stackalpha=stackalpha+stackalias3;
frame=stackimage./stackalpha;

if step~=steps
if step==1
stackdiff=false(localrezb(2),localrezb(1));
else
previousframe=abs(frame-previousframe);
stackdiff=previousframe(:,:,1)+previousframe(:,:,2)+previousframe(:,:,3)+stackdiff./stackalpha(:,:,2);
end
r=frame(:,:,1);
g=frame(:,:,2);
b=frame(:,:,3);
r=ordfilt2(r,1,ones(3,3),'symmetric');
g=ordfilt2(g,1,ones(3,3),'symmetric');
b=ordfilt2(b,1,ones(3,3),'symmetric');
mn=cat(3,r,g,b);
r=frame(:,:,1);
g=frame(:,:,2);
b=frame(:,:,3);
r=ordfilt2(r,9,ones(3,3),'symmetric');
g=ordfilt2(g,9,ones(3,3),'symmetric');
b=ordfilt2(b,9,ones(3,3),'symmetric');
mx=cat(3,r,g,b);
mx=mx-mn;
mx=mx(:,:,1).*2+mx(:,:,2).*3+mx(:,:,3);
stackalias=or(mx>aliasstep,stackdiff>alias);
end
previousframe=frame;

end
imwrite(uint16(frame(2:localrez(2)+1,2:localrez(1)+1,:).*2^16),['Z ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex),'.tif'])

tm=toc-tn;
hrs=floor(tm/(60*60));
mnt=floor((tm-hrs*60*60)/60);
snd=tm-hrs*60*60-mnt*60;
['Tile ',num2str(tilex+tiley.*tiles-tiles),' of ',num2str(tiles.^2),', ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex),', ',num2str(hrs),' hours, ',num2str(mnt),' minutes, ',num2str(snd),' seconds per tile.']

end
end
clear('a','a3','afull','aliasstep','ang','b','c','ca','cavg','cavgi','cavgr','cavgs','ci','cleanc','count','cr','cscale','fcs','frame','g','gby','gbys','gc','gm','go','grg','grgs','gs','hrs','inside','inside3','insidef','lcolor','lmag','localcenter','localrezb','loop','loopi','loopr','loops','maxn','mn','mnt','mx','n','o','outside','p','pixelsize','previousframe','r','snd','stackalias','stackalpha','stackdiff','stackimage','step','tic1','tilex','tiley','tlx','tly','tm','totalpixels','waitstepcheck','x','xGrid','xlim','y','yGrid','ylim','zi','zr');
try
finalimage=zeros(localrez(2).*tiles,localrez(1).*tiles,3,'uint16');
for tiley=1:tiles
for tilex=1:tiles
['Z ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex)]
end
end
imwrite(finalimage,'out.tif');
catch
finalimage=zeros(localrez(2),localrez(1).*tiles,3,'uint16');
for tiley=1:tiles
for tilex=1:tiles
['Z ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex)]
end
imwrite(finalimage,['F ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'.tif']);
end
end


### C++

C++ code by LionHeart[8]

/*
FWDDIFF.CPP a module for determining the slope using forward differencing calculations of fractals.

Written in Microsoft Visual 'C++' by Paul de Leeuw.

https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/main/app_settings/BumpMapSettings.java
*/

#include	"slope.h"

/**************************************************************************
Calculate functions
**************************************************************************/

void	CSlope::DoSlopeFwdDiffFn(Complex *z, Complex *q)
{
Complex	sqr, temp;
Complex	t = { 1.0, 0.0 };
double	real_imag;
int		k;
int		degree = (int)param[5];
int		degree1, degree2;

//    PaletteStart = (int)param[4];

switch (subtype)
{
case 0:					// Mandelbrot
sqr.x = z->x * z->x;
sqr.y = z->y * z->y;
real_imag = z->x * z->y;
z->x = q->x + sqr.x - sqr.y;
z->y = q->y + real_imag + real_imag;
break;
// other cases removed
}
}

/**************************************************************************
Get the gradients in the x and y directions
**************************************************************************/

double	CSlope::getGradientX(double *wpixels, int index, int width)
{
int x = index % width;
double it = *(wpixels + index);

if (x == 0) {
return (*(wpixels + index + 1) - it) * 2;
}
else if (x == width - 1) {
return (it - *(wpixels + index - 1)) * 2;
}
else {
double diffL = it - *(wpixels + index - 1);
double diffR = it - *(wpixels + index + 1);
return diffL * diffR >= 0 ? 0 : diffL - diffR;
}
}

double	CSlope::getGradientY(double *wpixels, int index, int width, int height)
{
int y = index / width;
double it = *(wpixels + index);

if (y == 0) {
return (it - *(wpixels + index + width)) * 2;
}
else if (y == height - 1) {
return (*(wpixels + index - width) - it) * 2;
}
else {
double diffU = it - *(wpixels + index - width);
double diffD = it - *(wpixels + index + width);
return diffD * diffU >= 0 ? 0 : diffD - diffU;
}
}

/**************************************************************************
Brightness Scaling
**************************************************************************/

int	CSlope::changeBrightnessOfColorScaling(int rgb, double delta)
{
int	    new_color = 0;
double	bump_transfer_factor = param[0];		// future add to dialogue box

//    double mul = getBumpCoef(delta);
double  mul = (1.5 / (fabs(delta * bump_transfer_factor) + 1.5));

if (delta > 0) {
rgb ^= 0xFFFFFF;
int r = rgb & 0xFF0000;
int g = rgb & 0x00FF00;
int b = rgb & 0x0000FF;
int ret = (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);
new_color = 0xff000000 | (ret ^ 0xFFFFFF);
}
else {
int r = rgb & 0xFF0000;
int g = rgb & 0x00FF00;
int b = rgb & 0x0000FF;
new_color = 0xff000000 | (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);
}

return new_color;
}

/**************************************************************************
Slope Fractal
**************************************************************************/

int	CSlope::RunSlopeFwdDiff(HWND hwndIn, void(*plot)(WORD, WORD, DWORD), int user_data(HWND hwnd), char* StatusBarInfo, bool *ThreadComplete, int subtypeIn, int NumThreadsIn, int threadIn, Complex j, double mandel_width, double hor, double vert, double xgap, double ygap,
double rqlim, long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr, BYTE juliaflag)
{
Complex	z = 0.0;			// initial value for iteration Z0
Complex	c, q;
BigComplex	cBig, zBig, qBig;
BigDouble	BigTemp, BigTemp1;

double	log_zn, nu;
int		lastChecked = -1;

int		x, y, i;
DWORD	index;
double	iterations;

subtype = subtypeIn;
hwnd = hwndIn;
for (i = 0; i < NUMSLOPEDERIVPARAM; i++)
param[i] = paramIn[i];

StripWidth = xdots;
else

for (y = 0; y < ydots; y++)
{
if (BigNumFlag)
cBig.y = BigVert + BigWidth - Big_ygap * y;
else
c.y = vert + mandel_width - y * ygap;
double progress = (double)y / ydots;
if (int(progress * 100) != lastChecked)
{
lastChecked = int(progress * 10);
}
for (x = StripStart; x < StripStart + StripWidth; x++)
{
if (user_data(hwnd) < 0)	// trap user input
return -1;

c.x = hor + x * xgap;
{
if (juliaflag)
{
q = j;
z = c;
}
else
{
q = c;
z = 0.0;
}
}
iterations = 0.0;
for (;;)
{
iterations++;
if (iterations >= threshold)
break;
DoSlopeFwdDiffFn(&z, &q);
if ((z.x * z.x + z.y * z.y) >= rqlim)
break;
}
if (iterations < threshold)
{
if (BigNumFlag)
{
// sqrt of inner term removed using log simplification rules.
BigTemp = zBig.x * zBig.x + zBig.y * zBig.y;
BigTemp1 = BigTemp.BigLog();
log_zn = mpfr_get_d(BigTemp1.x, MPFR_RNDN) / 2;
nu = log(log_zn / log(2.0)) / log(2.0);
}
else
{
// sqrt of inner term removed using log simplification rules.
log_zn = log(z.x * z.x + z.y * z.y) / 2;
nu = log(log_zn / log(2.0)) / log(2.0);
}
iterations = iterations + 1 - nu;
}
index = ((DWORD)y * (DWORD)width) + (DWORD)x;
if (x >= 0 && x < xdots - 1 && y >= 0 && y < ydots - 1)
*(SlopeDataPtr + index) = iterations;
plot(x, y, (long)iterations);
}
}
return 0;
}

/**************************************************************************
Slope Fractal
**************************************************************************/

int	CSlope::RenderSlope(long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr)

{

for (int i = 0; i < NUMPARAM; i++)
param[i] = paramIn[i];

double	bumpMappingDepth = param[3];		// future add to dialogue box
double	bumpMappingStrength = param[4];		// future add to dialogue box
double	lightDirectionDegrees = param[2];	// future add to dialogue box
int		PaletteStart = (int)param[1];

double	iterations;
int		lastChecked = -1;
DWORD	index;
int		x, y;
unsigned char r, g, b;
RGBTRIPLE	colour;
int		modified;

lastChecked = -1;

sizeCorr = 0.0;
lightx = 0.0;
lighty = 0.0;

gradCorr = pow(2, (bumpMappingStrength - DEFAULT_BUMP_MAPPING_STRENGTH) * 0.05);
sizeCorr = ydots / pow(2, (MAX_BUMP_MAPPING_DEPTH - bumpMappingDepth) * 0.16);
lightAngleRadians = lightDirectionDegrees * PI / 180.0;

for (y = 0; y < ydots; y++)
{
for (x = 0; x < xdots; x++)
{
index = ((DWORD)y * (DWORD)width) + (DWORD)x;
if (SlopeDataPtr)
iterations = *(SlopeDataPtr + index);
else
return 0;			// oops, we don't actually have forward differencing

if (iterations >= threshold)
{				//  interior of Mandelbrot set = inside_color
colour.rgbtRed = (BYTE)TrueCol->InsideRed;		// M_waves
colour.rgbtGreen = (BYTE)TrueCol->InsideGreen;
colour.rgbtBlue = (BYTE)TrueCol->InsideBlue;
}
else
{
// modified = rgbs[index];
if (iterations < PaletteStart)
modified = 0x00FFFFFF;
else
modified = 0xFF000000 | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 0) << 16) | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 1) << 8) | *(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 2);
//		int	original_color = modified;		// not sure what this is for
if (dotp != 0)
{
modified = changeBrightnessOfColorScaling(modified, cosAngle * smoothGrad);
}
r = (modified >> 16) & 0xFF;
g = (modified >> 8) & 0xFF;
b = modified & 0xFF;
colour.rgbtRed = r;
colour.rgbtGreen = g;
colour.rgbtBlue = b;
}
outRGBpoint(Dib, x, y, colour);
}
}
return 0;
}


### GLSL

slope lighting in an OpenGLSL code shader by Softology[9]

#version 400

//the following uniform values are set by Visions of Chaos prior to shader execution
uniform vec2 resolution;
uniform vec3 palette[256];
uniform double xmin;
uniform double xmax;
uniform double ymin;
uniform double ymax;
uniform double bailout;
uniform int maxiters;
uniform int samplepixels;

double sqrsamplepixels=double(samplepixels*samplepixels);
double bailout_squared=double(bailout*bailout);
double magnitude,r1,r2,g1,g2,b1,b2,r,g,b,tweenval;
float realiters;
vec4 finalcol,col;
int superx,supery;
double stepx=(xmax-xmin)/resolution.x/double(samplepixels);
double stepy=(ymax-ymin)/resolution.y/double(samplepixels);
int index,colval,colval2;
dvec2 z,c,dc,der,u,v,ctwo;
double h2,angle,pi,t;

//------------------------------------------------------------
// complex number operations
dvec2 cadd( dvec2 a, dvec2 b ) { return dvec2( a.x+b.x, a.y+b.y ); }
dvec2 csub( dvec2 a, dvec2 b ) { return dvec2( a.x-b.x, a.y-b.y ); }
dvec2 cmul( dvec2 a, dvec2 b )  { return dvec2( a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x ); }
dvec2 cdiv( dvec2 a, dvec2 b )  { double d = dot(b,b); return dvec2( dot(a,b), a.y*b.x - a.x*b.y ) / d; }
dvec2 csqr( dvec2 a ) { return dvec2(a.x*a.x-a.y*a.y, 2.0*a.x*a.y ); }
dvec2 csqrt( dvec2 z ) { double m = length(z); return sqrt( 0.5*dvec2(m+z.x, m-z.x) ) * dvec2( 1.0, sign(z.y) ); }
dvec2 conj( dvec2 z ) { return dvec2(z.x,-z.y); }
dvec2 cabs( dvec2 c) { return dvec2(sqrt(c.x * c.x + c.y * c.y)); }
//------------------------------------------------------------

void main(void)
{
finalcol=vec4(0,0,0,0);
pi=3.14159265359;
h2=1.5; // height factor of the incoming light
angle=45; // incoming direction of light
v=dvec2(exp(float(1.0*angle*2*pi/360)),0.0); // unit 2D vector in this direction
ctwo=dvec2(2.0,0.0);

for (supery=0;supery<samplepixels;supery++)
{
for (superx=0;superx<samplepixels;superx++)
{
c.x = xmin+gl_FragCoord.x/resolution.x*(xmax-xmin)+(stepx*double(superx));
c.y = ymin+gl_FragCoord.y/resolution.y*(ymax-ymin)+(stepy*double(supery));
int i;
z = dvec2(0.0,0.0);
dc = dvec2(1.0,0.0);
der = dc;
for(i=0; i<maxiters; i++)
{
//START OF FRACTAL FORMULA
//END OF FRACTAL FORMULA

magnitude=(z.x * z.x + z.y * z.y);
if(magnitude>bailout_squared) break;

// der = der*2*z + dc
}

if (i==maxiters) {
col=vec4(0.0,0.0,0.0,1.0);
} else {
//CPM smooth colors
//note that double precision does not support log so it needs to be cast as float
realiters=float(i+1-((log(log(sqrt(float(magnitude))))/log(2.0))));
colval=int(mod(realiters,255));
colval2=int(mod(colval+1,255));
tweenval=realiters-int(realiters);
r1=palette[colval].r;
g1=palette[colval].g;
b1=palette[colval].b;
r2=palette[colval2].r;
g2=palette[colval2].g;
b2=palette[colval2].b;
r=r1+((r2-r1)*tweenval);
g=g1+((g2-g1)*tweenval);
b=b1+((b2-b1)*tweenval);

u = cdiv(z,der);
u = cdiv(u,cabs(u)); // normal vector: (u.re,u.im,1)
t = u.x*v.x + u.y*v.y + h2; // dot product with the incoming light
t = t/(1+h2); // rescale so that t does not get bigger than 1
if (t<0) { t=0; }

//p.color = linear_interpolation(black,white,t)
t=1.0-t;
r=r*t;
g=g*t;
b=b*t;

col=vec4(r,g,b,1.0);
}

finalcol+=col;
}
}
gl_FragColor = vec4(finalcol/double(sqrsamplepixels));
}


## programs

• ultrafractal[10]
• directional DE, which can be used for slope colouring

## References

1. kerrymitchellart : uf1-5
2. Exploring the Mandelbrot Set by Eamonn O'Brien-Strain
3. Clapham, C.; Nicholson, J. (2009). "Oxford Concise Dictionary of Mathematics, Gradient" (PDF). Addison-Wesley. p. 348. Archived from the original (PDF) on 29 October 2013. Retrieved 1 September 2013.
4. hiddendimension Slope_Tutorial_2
5. hiddendimension Slope_Tutorial_3
6. hiddendimension tutorials: Slope_Tutorial_9
7. hiddendimension slope tutorial 7
8. fractalforums.org : algorithm-for-rendering-slope-in-a-fractal
9. fractalforums.org : algorithm-for-rendering-slope-in-a-fractal
10. ultrafractal standard formulas : slope