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 Mitchellkerrymitchellart : uf1-5 names lifted domain coloring Hill-Shading the Mandelbrot SetExploring the Mandelbrot Set by Eamonn O'Brien-Strain almondbread Hill-Shading the Mandelbrot Set 29 stycznia 2021 - article 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 right|thumb|Slope: In mathematics, the slope or gradient of a line is a number that describes both the direction and the steepness of the line.http://web.cortland.edu/matresearch/OxfordDictionaryMathematics.pdf Oxford Concise Dictionary of Mathematics, Gradient C. Clapham J. Nicholson Addison-Wesley 2009 348 1 September 2013 dead https://web.archive.org/web/20131029203826/http://web.cortland.edu/matresearch/OxfordDictionaryMathematics.pdf 29 October 2013 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 Mandelbrothiddendimension Slope_Tutorial_2 Julia Newton hiddendimension Slope_Tutorial_3 Kleinian group fractalshiddendimension tutorials: Slope_Tutorial_9 strange attractors hiddendimension slope tutorial 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 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 LionHeartfractalforums.org : algorithm-for-rendering-slope-in-a-fractal /* 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 Softologyfractalforums.org : algorithm-for-rendering-slope-in-a-fractal #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 ultrafractalultrafractal standard formulas : slope directional DE, which can be used for slope colouring Compare with bumpmaping Derivative References