分形/复平面上的迭代/斜率
< 分形
这是一种将 2D 图像转换为 3D 图像的算法
"每个像素迭代三个点，然后用于创建三维曲面。垂直于该曲面的向量用作 z 变量，并传递给着色例程。当使用斜率公式进行光照着色时，结果类似于从上方和侧面照亮的分形曲面。" Kerry Mitchell[1]
- 提升域着色
- 曼德勃罗集的阴影山体[2]
- 如何使用 GDAL 快速计算整个世界的阴影浮雕，作者：Anna Thieme
- 来自 ArcGIS Pro 3.1 的阴影山体工作原理
差异
- 阴影浮雕或阴影山体：通过考虑光源角和阴影，从表面栅格创建阴影浮雕。
- 斜率 = 识别栅格每个像元的斜率（梯度或陡峭度）。
在数学中，直线的斜率或梯度是一个描述直线方向和陡峭程度的数字。[3] 斜率通常用字母m表示；关于为什么用字母m表示斜率没有明确的答案。
y = mx + b
经过两点：(x1,y1) 和 (x2,y2) 的直线的斜率是
直线的斜率m与其倾斜角θ的关系由正切函数给出
- 由 Damien Jones 在 dmj.ufm 和 standard.ufm 中提供的 Ultrafractal 公式
-
朱利亚集的斜率 -
-
带有混合梯度平方根步长和斜率的距离曼德勃罗集的斜率。png -
hue = fmod(1 + carg(de) / (2 * pi), 1); // slope of complex number de
由 nickspiker 提供的 Matlab 代码
%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 [email protected]
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
由 LionHeart 提供的 C++ 代码[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;
}
由 Softology 提供的 OpenGLSL 代码着色器中的斜率光照[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));
}
- Ultrafractal[10]
- 方向性 DE，可用于斜率着色
- ↑ kerrymitchellart : uf1-5
- ↑ Eamonn O'Brien-Strain 的《探索曼德勃罗集》
- ↑ Clapham, C.; Nicholson, J. (2009). "牛津简明数学词典，梯度" (PDF). Addison-Wesley. p. 348. 从 原文 (PDF) 存档于 2013 年 10 月 29 日. 检索于 2013 年 9 月 1 日.
- ↑ hiddendimension 斜率教程 2
- ↑ hiddendimension 斜率教程 3
- ↑ hiddendimension 教程：斜率教程 9
- ↑ hiddendimension 斜率教程 7
- ↑ fractalforums.org : 在分形中渲染斜率的算法
- ↑ fractalforums.org : 在分形中渲染斜率的算法
- ↑ ultrafractal 标准公式 : 斜率