Fractals/Iterations in the complex plane/slope
Jump to navigation
Jump to search
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
- lifted domain coloring
- Hill-Shading the Mandelbrot Set[2]
- How to calculate shaded relief for the entire world with GDAL, fast by Anna Thieme
- How Hillshade works from ArcGIS Pro 3.1
Differences
- shaded relief or hill shade: Creates a shaded relief from a surface raster by considering the illumination source angle and shadows
- 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:
the slope m of a line is related to its angle of inclination θ by the tangent function:
Apply to
code
- Ultrafractal formulas by Damien Jones in dmj.ufm and standard.ufm
c
-
slope of Julia set -
-
slope of the distanceMandelbrot set with blended gradients sqrt step and slope.png -
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
adcs=zs./zsmallest;
zsmallest=zs;
if adcs<.25
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;
shaderotation=-pi./5;
slopescale=.5;
aic=abs(ic);
isreflect=sign(is).*aic;
shading=real(isreflect.*exp(shaderotation.*i))./3+.75;
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;
shading=1-(1-aics).*shading;
cc=1-cc.*(1-shading);
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)]
finalimage((tiley-1).*localrez(2)+1:(tiley).*localrez(2),(tilex-1).*localrez(1)+1:(tilex).*localrez(1),:)=imread(['Z ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex),'.tif']);
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)]
finalimage(:,(tilex-1).*localrez(1)+1:(tilex).*localrez(1),:)=imread(['Z ',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tiley),'x',sprintf(['%0',num2str(floor(log10(tiles))+1),'.0f'],tilex),'.tif']);
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/core/ThreadDraw.java#L3640
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;
thread = threadIn;
NumThreads = NumThreadsIn;
subtype = subtypeIn;
hwnd = hwndIn;
for (i = 0; i < NUMSLOPEDERIVPARAM; i++)
param[i] = paramIn[i];
if (NumThreads == 0)
StripWidth = xdots;
else
StripWidth = xdots / NumThreads;
StripStart = StripWidth * thread;
*ThreadComplete = false;
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);
sprintf(StatusBarInfo, "Progess (%d%%), %d Threads", int(progress * 100), NumThreads);
}
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);
}
}
*ThreadComplete = true;
return 0;
}
/**************************************************************************
Slope Fractal
**************************************************************************/
int CSlope::RenderSlope(long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr)
{
double dotp, gradAbs, gradCorr, cosAngle, sizeCorr, smoothGrad, lightAngleRadians, lightx, lighty;
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;
double gradx, grady;
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;
lightx = cos(lightAngleRadians) * gradCorr;
lighty = sin(lightAngleRadians) * gradCorr;
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);
gradx = getGradientX(SlopeDataPtr, index, xdots);
grady = getGradientY(SlopeDataPtr, index, xdots, ydots);
dotp = gradx * lightx + grady * lighty;
// int original_color = modified; // not sure what this is for
if (dotp != 0)
{
gradAbs = sqrt(gradx * gradx + grady * grady);
cosAngle = dotp / gradAbs;
smoothGrad = -2.3562 / (gradAbs * sizeCorr + 1.5) + 1.57;
//smoothGrad = Math.atan(gradAbs * sizeCorr);
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
z=cadd(csqr(z),c);
//END OF FRACTAL FORMULA
magnitude=(z.x * z.x + z.y * z.y);
if(magnitude>bailout_squared) break;
// der = der*2*z + dc
der = cadd(cmul(der,cmul(ctwo,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
Compare with
References
- ↑ kerrymitchellart : uf1-5
- ↑ Exploring the Mandelbrot Set by Eamonn O'Brien-Strain
- ↑ Template:Cite web
- ↑ hiddendimension Slope_Tutorial_2
- ↑ hiddendimension Slope_Tutorial_3
- ↑ hiddendimension tutorials: Slope_Tutorial_9
- ↑ hiddendimension slope tutorial 7
- ↑ fractalforums.org : algorithm-for-rendering-slope-in-a-fractal
- ↑ fractalforums.org : algorithm-for-rendering-slope-in-a-fractal
- ↑ ultrafractal standard formulas : slope