Fractals/Iterations in the complex plane/slope

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

Description

edit

"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

edit

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

edit
 
Slope:  

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

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

code

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



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

matlab

edit

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++ 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

edit

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

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

Compare with

edit

References

edit
  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