### Algorithm for Rendering Slope in a Fractal

• 111 Replies
• 1497 Views

0 Members and 1 Guest are viewing this topic.

#### C0ryMcG #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #90 on: March 22, 2020, 11:05:54 PM »
I see that someone posted as I was typing, and said the opposite of what I said, and posted nicer looking images than the ones I'm about to post... so now my confidence is a litlle shaken. but IF I understood the question properly, here's a rather rough illustration of what the process I described can do... this is using a heightmap generated through log potential calculations, and the following image is the slope-blur result.

#### xenodreambuie #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #91 on: March 23, 2020, 05:27:22 AM »
Hi C0ryMcG,
not really saying the opposite. I was answering Adam's question specifically about field lines for that field, while you were addressing the general question of turning it into a slope.
However, you're overcomplicating it a bit. Your first step gives you dx and dy at each pixel. Then the simplest step is to calculate a normal from (dx,dy,1) by normalizing it, then the z component is a scalar value of slope you can map directly to color brightness.
What I use the normal for instead is calculating the effect of light sources and material properties.

#### LionHeart #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #92 on: March 23, 2020, 05:32:29 AM »
Hi sergioct,

I'm struggling a bit. I can get the Nova fractal okay, but the derivative seems to be difficult to find. I tried the Quotient Rule for derivative: (f/g)' = (f'g - g'f) / g*g. But it gives incorrect images. This is definitely one disadvantage of the derivative method Paul the LionHeart

#### xenodreambuie #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #93 on: March 23, 2020, 09:08:57 AM »
Hi LionHeart,

the Nova formula is convergent. The usual form of distance estimate is correct for exterior regions. It can still be used, as it is simpler than the proper distance estimate for interior regions, but it is less correct and doesn't always work well, especially for Mandelbrots. Potential can a better source for slope in such cases.
Checking it with forward differencing, it does actually work ok, though with some discontinuities in the circular area on the left. This image uses potential for coloring and simple distance estimate for the height field. I do use analytic derivatives for calculating the distance estimate whenever it's reasonable, just not for calculating the slopes.

#### LionHeart #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #94 on: March 23, 2020, 09:40:33 AM »
Thanks xenodreambuie,

Quote
Potential can a better source for slope in such cases.

Do you have any pseudo code to show how to do this?

Many thanks.

#### xenodreambuie #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #95 on: March 23, 2020, 10:22:52 AM »
Hi LionHeart,

if you're already using continuous potential, you know how the formula depends on the type of critical point. If not, well, it's a bit complicated. I've never thought about how to get an analytic derivative of it, as I've always used forward differencing for the slope, and I've only read about others doing that too.

#### LionHeart #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #96 on: March 23, 2020, 11:05:34 AM »
Okay, thanks xenodreambuie.

I'll keep experimenting and see what I can come up with.

#### gerrit #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #97 on: March 27, 2020, 03:45:47 AM »
$$e^{1/z} z^2+c$$ using signed distance estimate directly (zlog(|z|)/z' instead of z/z') to make bumpmap.

#### nickspiker #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #98 on: March 29, 2020, 12:31:35 AM »
Interior and exterior slope utilized for shading and reflections.  All using derivatives mentioned earlier in this post.  And the interior slope was determined by differentiating over the period from Z=0
The force behind it is stronger than magic, it’s engineering.
~Ido

#### Adam Majewski #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #99 on: March 29, 2020, 09:50:15 AM »
Interior and exterior slope utilized for shading and reflections.  All using derivatives mentioned earlier in this post.  And the interior slope was determined by differentiating over the period from Z=0

cool

( I allways ask for it (:-))

Can you post the code ?
Can you upload it to commons ?

#### nickspiker #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #100 on: March 29, 2020, 06:45:19 PM »
cool

( I allways ask for it (:-))

Can you post the code ?
Can you upload it to commons ?

It's kind of a mess, but here you go!
It's Matlab code and terribly documented.  I'll clean it up at some point and put it on the commons.
Code: [Select]
%Mandelbrot rendering algorithms utilizing interior and exterior slope finding steps copyright 2020 Nick Spiker nick@nickspiker.commaxn=2^32;steps=2^5;%anti aliasing stepsmns=2^8;%max newton stepsdpi=90;bleed=0;%change to 6 when final render, 3 per sidewidth=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;endrez=[(width+bleed).*dpi,(height+bleed).*dpi];zoom=zoom./height.*(height+bleed);localrez=ceil(rez./tiles);%finds closest multiple equal to or greater than rezlocalrezb=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.']                    endendclear('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']);    endend

#### gerson #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #101 on: March 29, 2020, 07:27:53 PM »
@nickspiker All sharing code, ideas, etc are welcome here, thanks.

My gallery:
www.deviantart.com/gs13

#### C0ryMcG #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #102 on: March 29, 2020, 08:33:55 PM »
Interior and exterior slope utilized for shading and reflections.  All using derivatives mentioned earlier in this post.  And the interior slope was determined by differentiating over the period from Z=0
This is beautiful, and somewhat demoralizing... my experiments with the same idea have not turned out nearly so clean or attractive. I'll see if I can understand your code, though, and keep on trying. One thing I am not sure how to tackle is that when I try to get an interior slope, the orbits with odd-numbered periods render inverted.
Either way, well done, it's a fantastic render. I hope mine will approach that eventually.

#### nickspiker #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #103 on: March 30, 2020, 02:30:59 PM »
...odd-numbered periods render inverted.
Yup, definitely doing something wrong.  I think it may be how you are figuring the derivative maybe?
Here's a link from the post earlier detailing the derivative portion:
https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Normal_map_effect

I do the shading bit a little different but the calculations to find the derivative are the same
Code: [Select]
%for each pixel do the stuff belowc=somethingreal+somethingimaginary;z=0;d=1;while (real(z)^2+imaginary(z)^2)<realbignumberz=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 oneendwhileslope=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 onsunangle=45;brightness=slope*e^(i*sunangle*pi/180);%this rotates the sign by the angle abovebrightness=real(brightness)/2+.5;Hopefully the re-coded stuff above clears it up some?
You can always send over what you have and I can look it over too.

#### LionHeart #### Re: Algorithm for Rendering Slope in a Fractal

« Reply #104 on: March 30, 2020, 09:27:57 PM »
Dear Friends,

I finally got there. I went back to scratch and rewrote the forward differencing method for slope calculation. And it works. I wasn't able to apply the derivative method to burning ship because it has discontinuities and I was unable to find a derivative. So I persevered with the alternative method. Many thanks to you all, especially superheal who hung in there and even debugged my code.

Wow nickspiker,

Those are amazing images. Thanks for sharing your code with us.

### Similar Topics ###### Algorithm for Full Exploration of 3D Fractal in Realtime

Started by AranLeer on Fractal Mathematics And New Theories

13 Replies
521 Views April 26, 2019, 07:34:08 AM ###### How would quantum computing affect fractal rendering?

Started by greentexas on Programming

1 Replies
584 Views December 12, 2017, 01:40:20 AM
by eiffie ###### Triangle Inequality Average Algorithm

Started by mikoval on Fractal Mathematics And New Theories

10 Replies
623 Views August 05, 2018, 02:36:30 PM
by FractalDave ###### Chaos to speed up numerical algorithm

Started by gerrit on Physics

1 Replies
424 Views February 15, 2018, 12:40:28 PM
by knighty ###### ultrafractal - translationg coloring algorithm ( ucl) to c

Started by Adam Majewski on Programming

4 Replies
490 Views January 05, 2018, 05:39:20 PM