MATLAB非线性可视化之Mandelbrot集与分形

%Mandelbrot
%?res是目标分辨率，iter是循环次数，（xc,yc）是图像中心，xoom是放大倍数
%res分辨率越多，在一张图上能够计算的点越多。iter越多，细节越丰富。
clear;
clc
close all
% res=1024;iter=100;xc=-0.748766710846959;yc=-0.123640847970064;xoom=6.503637847981780e+08;
% res=2048;iter=10000;xc=-0.912;yc=-0.2611;xoom=150;%为了生成图形好看，计算时间非常长
res=1024;iter=100;xc=-0.5;yc=0;xoom=1;
% res=2500;iter=200;xc=0;yc=0;xoom=1.7;

x0=xc-2/xoom;x1=xc+2/xoom;
y0=yc-2/xoom;y1=yc+2/xoom;
x=linspace(x0,x1,res);
y=linspace(y0,y1,res);
[xx,yy]=meshgrid(x,y);
z=xx+yy*1i;
C=z;
% C=0.279;
tic

%matlab官方推荐
N = ones(res,res);
%N=uint32(N);
for k = 0:iter
    z = z.*z + C;
    inside = abs( z )<=2;
    %N = N + uint32(inside);
    N = N + inside;
end
toc



%蓝色-橘色颜色条
mycolorpoint=[
    [235 255 255];...
    [243 253 242];...
    [248 247 193];...
    [255 211 49];...
    [236 133 0];...
    [134 36 21];...
    [75 6 34];...
    [15 5 64];...
    [3 8 96];...
    [7 35 154];...
    [24 84 205];...
    [87 169 229];...
    [165 223 247];...
    [235 255 255];...
    ];
mycolorposition=[1 3 7 14 27 45 57 66 74 86 96 109 119 128];
mycolormap_r=interp1(mycolorposition,mycolorpoint(:,1),1:128,'pchip','extrap');
mycolormap_g=interp1(mycolorposition,mycolorpoint(:,2),1:128,'pchip','extrap');
mycolormap_b=interp1(mycolorposition,mycolorpoint(:,3),1:128,'pchip','extrap');
mycolor=[mycolormap_r',mycolormap_g',mycolormap_b']/255;


figure()
pcolor(x,y,log(N));

mycolor=interp1(linspace(1,256,128),mycolor,1:256);
colormap([mycolor;mycolor;mycolor])
shading interp

Nq=LowPrecision2High(N);
figure()
pcolor(x,y,log(Nq));
colormap([mycolor;mycolor;mycolor])
shading interp

function Zq=LowPrecision2High(Z)
%对阶梯数据进行平滑
[Nx,Ny]=size(Z);
[X,Y]=meshgrid(1:Nx,1:Ny);
XP=[];YP=[];ZP=[];
for kx=1:(Nx-1)
    for ky=1:(Ny-1)
        if Z(kx,ky)~=Z(kx,ky+1)
            x_t=0.5*(X(kx,ky)+X(kx,ky+1));
            y_t=0.5*(Y(kx,ky)+Y(kx,ky+1));
            z_t=0.5*(Z(kx,ky)+Z(kx,ky+1));
            XP=[XP;x_t];YP=[YP;y_t];ZP=[ZP;z_t];
        end
        if Z(kx,ky)~=Z(kx+1,ky)
            x_t=0.5*(X(kx,ky)+X(kx+1,ky));
            y_t=0.5*(Y(kx,ky)+Y(kx+1,ky));
            z_t=0.5*(Z(kx,ky)+Z(kx+1,ky));
            XP=[XP;x_t];YP=[YP;y_t];ZP=[ZP;z_t];
            %pause(0.1)
        end
    end
end
%边界值
for kx=1:Nx
    XP=[XP;X(kx,1)];  YP=[YP;Y(kx,1)];  ZP=[ZP;Z(kx,1)];
    XP=[XP;X(kx,end)];YP=[YP;Y(kx,end)];ZP=[ZP;Z(kx,end)];
end
for ky=2:Ny-1
    XP=[XP;X(1,ky)];  YP=[YP;Y(1,ky)];  ZP=[ZP;Z(1,ky)];
    XP=[XP;X(end,ky)];YP=[YP;Y(end,ky)];ZP=[ZP;Z(end,ky)];
end
Zq=griddata(XP,YP,ZP,X,Y,'natural');%'v4'最好，实在不行用linear
end