首先看一下这个相关的问答:
主要的想法是使用直方图来更有效地将颜色渐变分配给已使用的索引,而不是在未使用的索引上统一浪费许多颜色。此外,它还使用了特定的视觉上令人愉悦的渐变函数:
其他人提出的动态最大迭代次数只会影响整体性能和缩放中的细节。然而,如果你希望在没有缩放的情况下获得漂亮的颜色,那么你需要计算浮点型迭代次数,也称为Mandelbrot Escape。有一种数学方法可以从方程式的最后子结果计算出迭代次数的小数部分。更多信息请参见:
然而,我从未尝试过,所以带着偏见阅读此内容:如果我理解正确,你想计算这个方程式:
mu = m + frac = n + 1 - log (log |Z(n)|) / log 2
当你迭代次数为n
时,Z(n)
是你正在迭代的方程的复数域子结果。现在从mu
计算颜色,它现在是浮点型而不是从n
计算...
[编辑2]基于上述链接的GLSL mandelbrot与分数逃逸
我添加了分数逃逸并修改了直方图多通道重新着色以匹配新输出...
顶点:
#version 420 core
layout(location=0) in vec2 pos;
out smooth vec2 p;
void main()
{
p=pos;
gl_Position=vec4(pos,0.0,1.0);
}
片段:
#version 420 core
uniform vec2 p0=vec2(0.0,0.0);
uniform float zoom=1.000;
uniform int n=100;
uniform int sh=7;
uniform int multipass=0;
in smooth vec2 p;
out vec4 col;
const int n0=1;
vec3 spectral_color(float l)
{
float t; vec3 c=vec3(0.0,0.0,0.0);
if ((l>=400.0)&&(l<410.0)) { t=(l-400.0)/(410.0-400.0); c.r= +(0.33*t)-(0.20*t*t); }
else if ((l>=410.0)&&(l<475.0)) { t=(l-410.0)/(475.0-410.0); c.r=0.14 -(0.13*t*t); }
else if ((l>=545.0)&&(l<595.0)) { t=(l-545.0)/(595.0-545.0); c.r= +(1.98*t)-( t*t); }
else if ((l>=595.0)&&(l<650.0)) { t=(l-595.0)/(650.0-595.0); c.r=0.98+(0.06*t)-(0.40*t*t); }
else if ((l>=650.0)&&(l<700.0)) { t=(l-650.0)/(700.0-650.0); c.r=0.65-(0.84*t)+(0.20*t*t); }
if ((l>=415.0)&&(l<475.0)) { t=(l-415.0)/(475.0-415.0); c.g= +(0.80*t*t); }
else if ((l>=475.0)&&(l<590.0)) { t=(l-475.0)/(590.0-475.0); c.g=0.8 +(0.76*t)-(0.80*t*t); }
else if ((l>=585.0)&&(l<639.0)) { t=(l-585.0)/(639.0-585.0); c.g=0.84-(0.84*t) ; }
if ((l>=400.0)&&(l<475.0)) { t=(l-400.0)/(475.0-400.0); c.b= +(2.20*t)-(1.50*t*t); }
else if ((l>=475.0)&&(l<560.0)) { t=(l-475.0)/(560.0-475.0); c.b=0.7 -( t)+(0.30*t*t); }
return c;
}
void main()
{
int i,j,N;
vec2 pp;
float x,y,q,xx,yy,mu;
pp=(p/zoom)-p0;
pp.x-=0.5;
for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n-n0)&&(xx+yy<4.0);i++)
{
q=xx-yy+pp.x;
y=(2.0*x*y)+pp.y;
x=q;
xx=x*x;
yy=y*y;
}
for (j=0;j<n0;j++,i++)
{
q=xx-yy+pp.x;
y=(2.0*x*y)+pp.y;
x=q;
xx=x*x;
yy=y*y;
}
mu=float(i)-log(log(sqrt(xx+yy))/log(2.0));
mu*=float(1<<sh); i=int(mu);
N=n<<sh;
if (i>N) i=N;
if (i<0) i=0;
if (multipass!=0)
{
float r,g,b;
r= i &255; r/=255.0;
g=(i>> 8)&255; g/=255.0;
b=(i>>16)&255; b/=255.0;
col=vec4(r,g,b,255);
}
else{
q=float(i)/float(N);
q=pow(q,0.2);
col=vec4(spectral_color(400.0+(300.0*q)),1.0);
}
}
CPU侧C++/VCL代码:
#include <vcl.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl\\OpenGL3D_double.cpp"
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
OpenGLscreen scr;
GLSLprogram shd;
float mx=0.0,my=0.0,mx0=0.0,my0=0.0,mx1=0.0,my1=0.0;
TShiftState sh0,sh1;
int xs=1,ys=1;
float zoom=1.000;
int sh=7;
int N=256;
int _multi=0;
unsigned int queryID[2];
#define multi_pass
OpenGLtexture txr;
DWORD spectral_color(float l)
{
float t; float r,g,b; DWORD c,x; r=0.0; g=0.0; b=0.0;
if ((l>=400.0)&&(l<410.0)) { t=(l-400.0)/(410.0-400.0); r= +(0.33*t)-(0.20*t*t); }
else if ((l>=410.0)&&(l<475.0)) { t=(l-410.0)/(475.0-410.0); r=0.14 -(0.13*t*t); }
else if ((l>=545.0)&&(l<595.0)) { t=(l-545.0)/(595.0-545.0); r= +(1.98*t)-( t*t); }
else if ((l>=595.0)&&(l<650.0)) { t=(l-595.0)/(650.0-595.0); r=0.98+(0.06*t)-(0.40*t*t); }
else if ((l>=650.0)&&(l<700.0)) { t=(l-650.0)/(700.0-650.0); r=0.65-(0.84*t)+(0.20*t*t); }
if ((l>=415.0)&&(l<475.0)) { t=(l-415.0)/(475.0-415.0); g= +(0.80*t*t); }
else if ((l>=475.0)&&(l<590.0)) { t=(l-475.0)/(590.0-475.0); g=0.8 +(0.76*t)-(0.80*t*t); }
else if ((l>=585.0)&&(l<639.0)) { t=(l-585.0)/(639.0-585.0); g=0.84-(0.84*t) ; }
if ((l>=400.0)&&(l<475.0)) { t=(l-400.0)/(475.0-400.0); b= +(2.20*t)-(1.50*t*t); }
else if ((l>=475.0)&&(l<560.0)) { t=(l-475.0)/(560.0-475.0); b=0.7 -( t)+(0.30*t*t); }
r*=255.0; g*=255.0; b*=255.0;
x=r; c =x;
x=g; c|=x<<8;
x=b; c|=x<<16;
return c;
}
void gl_draw()
{
scr.cls();
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glMatrixMode(GL_TEXTURE);
glLoadIdentity();
shd.bind();
shd.set2f("p0",mx,my);
shd.set1f("zoom",zoom);
shd.set1i("n",N);
shd.set1i("sh",sh);
shd.set1i("multipass",_multi);
glQueryCounter(queryID[0], GL_TIMESTAMP);
glColor3f(1.0,1.0,1.0);
glBegin(GL_QUADS);
glVertex2f(-1.0,+1.0);
glVertex2f(-1.0,-1.0);
glVertex2f(+1.0,-1.0);
glVertex2f(+1.0,+1.0);
glEnd();
shd.unbind();
if (_multi)
{
float t,m,n=N<<sh;
DWORD *hist=new DWORD[n+1];
int sz=txr.xs*txr.ys,i,j;
glReadPixels(0,0,txr.xs,txr.ys,GL_RGBA,GL_UNSIGNED_BYTE,txr.txr);
for (i=0;i<=n;i++) hist[i]=0;
for (i=0;i<sz;i++) hist[txr.txr[i]&0x00FFFFFF]++;
for (i=1,j=1;i<=n;i++)
if (hist[i]){ hist[i]=j; j++; }
m=1.0/float(j); hist[0]=0x00000000;
for (i=1;i<=n;i++)
if (hist[i]){ t=hist[i]; t*=m; hist[i]=spectral_color(400.0+(300.0*t)); }
else hist[i]=0x00000000;
for (i=0;i<sz;i++) txr.txr[i]=hist[txr.txr[i]&0x00FFFFFF];
scr.cls();
txr.bind();
glColor3f(1.0,1.0,1.0);
glBegin(GL_QUADS);
glTexCoord2f(0.0,1.0); glVertex2f(-1.0,+1.0);
glTexCoord2f(0.0,0.0); glVertex2f(-1.0,-1.0);
glTexCoord2f(1.0,0.0); glVertex2f(+1.0,-1.0);
glTexCoord2f(1.0,1.0); glVertex2f(+1.0,+1.0);
glEnd();
txr.unbind();
glDisable(GL_TEXTURE_2D);
delete[] hist;
}
glQueryCounter(queryID[1], GL_TIMESTAMP);
scr.text_init_pix(0.75);
glColor4f(1.0,1.0,1.0,0.9);
scr.text(glGetAnsiString(GL_VENDOR));
scr.text(glGetAnsiString(GL_RENDERER));
scr.text("OpenGL ver: "+glGetAnsiString(GL_VERSION));
if (_multi) scr.text("Multi pass");
else scr.text("Single pass");
if (shd.log.Length()!=41)
for (int i=1;i<=shd.log.Length();) scr.text(str_load_lin(shd.log,i,true));
scr.text_exit();
scr.exe();
scr.rfs();
int e;
unsigned __int64 t0,t1;
for (e=0;!e;) glGetQueryObjectiv(queryID[0],GL_QUERY_RESULT_AVAILABLE,&e);
for (e=0;!e;) glGetQueryObjectiv(queryID[1],GL_QUERY_RESULT_AVAILABLE,&e);
glGetQueryObjectui64v(queryID[0], GL_QUERY_RESULT, &t0);
glGetQueryObjectui64v(queryID[1], GL_QUERY_RESULT, &t1);
Form1->Caption=AnsiString().sprintf("dt: %f ms p0:%.3fx%.3f zoom: %.1lf N:%i<<%i\n",(t1-t0)/1000000.0,mx,my,zoom,N,sh);
}
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
scr.init(this);
shd.set_source_file("","","","Mandelbrot_set.glsl_vert","Mandelbrot_set.glsl_frag");
glGenQueries(2, queryID);
_multi=1;
zoom=300.0;
mx = 0.268;
my =-0.102;
}
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
scr.exit();
}
void __fastcall TForm1::FormResize(TObject *Sender)
{
scr.resize();
xs=ClientWidth;
ys=ClientHeight;
txr.resize(xs,ys);
gl_draw();
}
void __fastcall TForm1::FormPaint(TObject *Sender)
{
gl_draw();
}
void __fastcall TForm1::FormMouseMove(TObject *Sender, TShiftState Shift, int X,int Y)
{
bool q0,q1;
mx1=1.0-divide(X+X,xs-1);
my1=divide(Y+Y,ys-1)-1.0;
sh1=Shift;
q0=sh0.Contains(ssLeft);
q1=sh1.Contains(ssLeft);
if (q1)
{
mx-=(mx1-mx0)/zoom;
my-=(my1-my0)/zoom;
}
mx0=mx1; my0=my1; sh0=sh1;
gl_draw();
}
void __fastcall TForm1::FormMouseDown(TObject *Sender, TMouseButton Button,TShiftState Shift, int X, int Y)
{
FormMouseMove(Sender,Shift,X,Y);
}
void __fastcall TForm1::FormMouseUp(TObject *Sender, TMouseButton Button,TShiftState Shift, int X, int Y)
{
FormMouseMove(Sender,Shift,X,Y);
}
void __fastcall TForm1::FormMouseWheel(TObject *Sender, TShiftState Shift, int WheelDelta, TPoint &MousePos, bool &Handled)
{
if (WheelDelta>0) zoom*=1.2;
if (WheelDelta<0) zoom/=1.2;
Handled=true;
gl_draw();
}
void __fastcall TForm1::FormKeyDown(TObject *Sender, WORD &Key, TShiftState Shift)
{
Caption=Key;
if (Key==32){ _multi=!_multi; gl_draw(); }
if (Key==33){ if (N<8192) N<<=1; gl_draw(); }
if (Key==34){ if (N> 128) N>>=1; gl_draw(); }
}
这是单次分数逃逸 n=100*32
:
![fractional escape](https://istack.dev59.com/t1GMg.webp)
这是单次整数逃逸 n=100
:
![integer escape](https://istack.dev59.com/zz2mU.webp)
可以看到,相同迭代次数(100
)的分数逃逸更好。
最后是漂亮的多次逃逸(炫耀一下),只有256次迭代和~300倍缩放:
![multi pass](https://istack.dev59.com/FGV9R.webp)
与单次逃逸相比:
![single pass](https://istack.dev59.com/4wXnQ.webp)
对修改进行一些解释:
我在计数器(定点数)中添加了sh
个分数位。因此,最大计数现在是n<<sh
而不仅仅是n
。我还添加了n0
常量,降低了逃逸的分数部分的误差。该链接建议使用2次迭代,但我认为1次看起来更好(它也从对数方程中删除了i+1
增量)。迭代循环未更改,我只是在其中添加了相同的n0
次迭代,然后计算分数逃逸mu
并将其转换为定点数(因为我的着色器输出整数)。
多次逃逸仅在CPU端代码上进行了更改。它简单地重新索引已使用的索引,以便其中没有空洞,并使用可见光谱颜色重新着色。
这里有一个演示:
Color.getHSBColor(i / 255F, 1, 1)
,其中i是迭代到发散的次数。我建议在缩放时动态增加迭代次数。 - Daniel Williams