DLA分形的MATLAB仿真

直接使用MATLAB跑就是慢些,但跑百把万的粒子数目无压力,混合编程在100000点以下跑起来很快,差不多效率提升了12倍,但150000以上时明显变慢,顿时亚历山大了,估计与c的随机数有关,暂记于此吧,上一幅1000000个粒子的DLA分形图。

DLA(1000000)

function dla(n,size,r)
%DLA分形生长模型
%n为模拟粒子的数目,size为随机产生粒子的半径,宜大,r为团簇的半径
tic;
switch nargin
    case 0
        n=10000; r=5; size=512;
    case 1
        size=512; r=5;
    case 2
        r=512;
end
rmax=0;  %粒子团的最大半径
%矩阵m用于储存点的位置
m=zeros(size*2+1,size*2+1,'uint8');
m(size+1,size+1)=1;
%生成一个随机数序列
ran=single(rand(n,1)*2*pi);
t=1;
while t<n
    xp=round(cos(ran(t,1))*(rmax+r));
    yp=round(sin(ran(t,1))*(rmax+r));
    disp(t);    %显示计算的点
    %(2*rmax+r)*pi为粒子行走的最大步数,随半径增大
    %若走完步数还未接触到团簇或接触到边界,则会进入while循环重新模拟
    for step=1:int32((2*rmax+r)*pi)
        rans=rand;
        if rans<0.25
            xp=xp+1;
        elseif rans<0.50
            yp=yp+1;
        elseif rans<0.75
            xp=xp-1;
        else
            yp=yp-1;
        end
        %判断粒子是否触界,判断的话影响速率
        %if xp == size || yp == size
        %    break;
        %end
        %判断当前粒子是否接触到团簇
        if m(xp+size+2,yp+size+1)==1 || m(xp+size,yp+size+1)==1 || m(xp+size+1,yp+size+2)==1 || m(xp+size+1,yp+size)==1
            m(xp+size+1,yp+size+1)=1;
            rmax=max(sqrt(xp^2+yp^2),rmax);
            t=t+1;
            break;
        end
    end
end
%画图
colormap(gray);
image(255-255*m);
%axis off;
save(['dla(',num2str(n),').mat'],'m','size','r');
toc;
end
#include "mex.h" // 使用MEX文件必须包含的头文件
//#include<ctime>
#include<math.h>

//获取最大值
inline double max( double x , double y)
{
    if( x > y )
        return x;
    else
        return y;
}

void dla_c(int n, int size, int r, double *m)
{
    //定义团簇的生长点
    int xp = 0, yp = 0;
    *(m + (xp + size) * (size * 2 + 1) + size - yp) = 1;
    double rmax = 0;
    double ran;
    double rans;
    double pi = 3.1415926;
    int step, stepmax;
    for(int t = 1; t < n; t++)
    {
        //显示正在计算的粒子个数
        mexPrintf("%d\n", t);
        //srand(time(NULL));
        ran = double(rand() % 100);
        ran = ran / 50 * pi;
        xp = int(cos(ran) * (rmax + r) + 0.5);
        yp = int(sin(ran) * (rmax + r) + 0.5);
        //(rmax*2+r)*pi为粒子行走的最大步数,视具体模型而定
        //若走完步数还未接触到团簇,则会进入while循环重新模拟
        step = 0, stepmax = int((r + 2 * rmax) * pi);
        do
        {
            //srand(time(NULL));
            rans = rand() % 100;
            if(rans < 25)
                xp++;
            else if(rans < 50)
                yp++;
            else if(rans < 75)
                xp--;
            else
                yp--;
            step++;
            if(step == stepmax || abs(xp) >= size / 2 || abs(yp) >= size / 2)
            {
                xp = int(cos(ran) * (rmax + r) + 0.5);
                yp = int(sin(ran) * (rmax + r) + 0.5);
                step = 0;
            }
        }while(*(m + (xp + 1 + size) * (size * 2 + 1) + size - yp) == 0
                    && *(m + (xp - 1 + size) * (size * 2 + 1) + size - yp) == 0
                    && *(m + (xp + size) * (size * 2 + 1) + size - yp - 1) == 0
                    && *(m + (xp + size) * (size * 2 + 1) + size - yp + 1) == 0);
        *(m + (xp + size) * (size * 2 + 1) + size - yp) = 1;
        rmax = max(sqrt(double(xp * xp + yp * yp)), rmax);
    }
    return;
}
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
    //m用于储存分形图形
    double *m;
    //n为粒子的数量,size为预想DLA分形图形的半径(宜大),r为随机产生的粒子与团簇间的距离
    int n, size, r;
    n = int(*(mxGetPr(prhs[0])));
    size = int(*(mxGetPr(prhs[1])));
    r = int(*(mxGetPr(prhs[2])));
    
    plhs[0] = mxCreateDoubleMatrix(2 * size + 1, 2 * size + 1, mxREAL);
    m = mxGetPr(plhs[0]);
    
    dla_c(n, size, r, m);
}

Leave a Reply

Your email address will not be published. Required fields are marked *