2014-02-13 108 views
0

我有一個1024 * 1024 * 51的矩陣。我將進行計算以改變for循環內矩陣的某個值(更改每次迭代的矩陣值)。我發現計算速度變得越來越慢,最後我的計算機陷入困境。但矩陣的大小並沒有改變。任何人都可以解釋這個問題?矩陣計算在matlab中每次迭代後變得更慢

function ActiveContours3D(method,grad,im,mu,nu,lambda1,lambda2,TimeSteps) 
epsilon = 10e-10; 
tic 
fid=fopen('Chr18_z_25of25tiles-C=0_c0_n000.raw','rb','ieee-le'); 
Xdim = 1024; 
Ydim = 1024; 
Zdim = 51; 
A = fread(fid,[Xdim Ydim*Zdim],'int16'); 
A = double(A); 
size_of_A = size(A) 
for(i=1:Zdim) 
    u0_color(:,:,i) = A(1 : Xdim , (i-1)*Ydim+1 : i*Ydim); 
end 
fclose(fid) 

time = toc 

[M,N,P,color] = size(u0_color); 
size(u0_color); 
u0_color = double(u0_color); % Convert u0_color values to double; 
u0 = u0_color(:,:,:,1);   % Define the Grayscale volumetric image. 
u0_color = uint8(u0_color);  % Necessary for color visualization 
x = 1:M; 
y = 1:N; 
z = 1:P; 
dx = 1 
dy = 1 
dz = 1 
dim_approx = 2*M*N*P/sqrt(M*N*P); 
if(method == 'Explicit') 
    dt = 0.9/((2*mu/(dx^2)) + (2*mu/(dy^2)) + (2*mu/(dz^2))) % 90% CFL 
elseif(method == 'Implicit') 
    dt = (10e7) * 0.9/((2*mu/(dx^2)) + (2*mu/(dy^2)) + (2*mu/(dz^2))) 
end 
[X,Y,Z] = meshgrid(x,y,z); 
x0 = (M+1)/2; 
y0 = (N+1)/2; 
z0 = (P+1)/2; 
r0 = min(min(M,N),P)/3; 
phi = sqrt((X-x0).^2 + (Y-y0).^2 + (Z-z0).^2) - r0; 
phi_visualize = phi;   % Use this for visualization in 3D 
phi = permute(phi,[2,1,3]);  % Use this for computations in 3D 
write_to_binary_file(phi_visualize,0);  % record initial conditions 

tic 

for(n=1:TimeSteps) 
n 

c1 = C1_3d(u0,phi); 
c2 = C2_3d(u0,phi); 

% x 
phi_xp = [phi(2:M,:,:); phi(M,:,:)];  % vertical concatenation 
phi_xm = [phi(1,:,:); phi(1:M-1,:,:)];  % (since x values are rows) 
         % cat(1,A,B) is the same as [A;B] 
Dx_m = (phi - phi_xm)/dx;   % first derivatives 
Dx_p = (phi_xp - phi)/dx; 

Dxx = (Dx_p - Dx_m)/dx;  % second derivative 

% y 
phi_yp = [phi(:,2:N,:) phi(:,N,:)];  % horizontal concatenation 
phi_ym = [phi(:,1,:) phi(:,1:N-1,:)];   % (since y values are columns) 
         % cat(2,A,B) is the same as [A,B] 
Dy_m = (phi - phi_ym)/dy; 
Dy_p = (phi_yp - phi)/dy; 

Dyy = (Dy_p - Dy_m)/dy; 

% z 
phi_zp = cat(3,phi(:,:,2:P),phi(:,:,P)); 
phi_zm = cat(3,phi(:,:,1) ,phi(:,:,1:P-1)); 

Dz_m = (phi - phi_zm)/dz; 
Dz_p = (phi_zp - phi)/dz; 

Dzz = (Dz_p - Dz_m)/dz; 

% x,y,z 
Dx_0 = (phi_xp - phi_xm)/(2*dx); 
Dy_0 = (phi_yp - phi_ym)/(2*dy); 
Dz_0 = (phi_zp - phi_zm)/(2*dz); 

phi_xp_yp = [phi_xp(:,2:N,:) phi_xp(:,N,:)]; 
phi_xp_ym = [phi_xp(:,1,:) phi_xp(:,1:N-1,:)]; 
phi_xm_yp = [phi_xm(:,2:N,:) phi_xm(:,N,:)]; 
phi_xm_ym = [phi_xm(:,1,:) phi_xm(:,1:N-1,:)]; 

phi_xp_zp = cat(3,phi_xp(:,:,2:P),phi_xp(:,:,P)); 
phi_xp_zm = cat(3,phi_xp(:,:,1) ,phi_xp(:,:,1:P-1)); 
phi_xm_zp = cat(3,phi_xm(:,:,2:P),phi_xm(:,:,P)); 
phi_xm_zm = cat(3,phi_xm(:,:,1) ,phi_xm(:,:,1:P-1)); 

phi_yp_zp = cat(3,phi_yp(:,:,2:P),phi_yp(:,:,P)); 
phi_yp_zm = cat(3,phi_yp(:,:,1) ,phi_yp(:,:,1:P-1)); 
phi_ym_zp = cat(3,phi_ym(:,:,2:P),phi_ym(:,:,P)); 
phi_ym_zm = cat(3,phi_ym(:,:,1) ,phi_ym(:,:,1:P-1));  


if(grad == 'Dirac') 
    Grad = DiracDelta(phi);     % Dirac delta  
    %Grad = 1; 
elseif(grad == 'Grad ') 
    Grad = (((Dx_0.^2)+(Dy_0.^2)+(Dz_0.^2)).^(1/2)); % |grad phi| 
end 

if(method == 'Explicit') 
    % CURVATURE: *mu*k|grad phi|* (central differences): 
    K = zeros(M,N,P); 

    Dxy = (phi_xp_yp - phi_xp_ym - phi_xm_yp + phi_xm_ym)/(4*dx*dy); 
    Dxz = (phi_xp_zp - phi_xp_zm - phi_xm_zp + phi_xm_zm)/(4*dx*dz); 
    Dyz = (phi_yp_zp - phi_yp_zm - phi_ym_zp + phi_ym_zm)/(4*dy*dz); 

    K = ((Dx_0.^2).*Dyy - 2*Dx_0.*Dy_0.*Dxy + (Dy_0.^2).*Dxx ... 
     + (Dx_0.^2).*Dzz - 2*Dx_0.*Dz_0.*Dxz + (Dz_0.^2).*Dxx ... 
     + (Dy_0.^2).*Dzz - 2*Dy_0.*Dz_0.*Dyz + (Dz_0.^2).*Dyy) ./ ((Dx_0.^2 + Dy_0.^2 + Dz_0.^2).^(3/2) + epsilon); 

    phi_temp = phi + dt * Grad .* (mu.*K + lambda1*(u0 - c1).^2 - lambda2*(u0 - c2).^2); 

elseif(method == 'Implicit') 
    C1x = 1 ./ sqrt(Dx_p.^2 + Dy_0.^2 + Dz_0.^2 + (10e-7)^2); 
    C2x = 1 ./ sqrt(Dx_m.^2 + Dy_0.^2 + Dz_0.^2 + (10e-7)^2); 
    C3y = 1 ./ sqrt(Dx_0.^2 + Dy_p.^2 + Dz_0.^2 + (10e-7)^2); 
    C4y = 1 ./ sqrt(Dx_0.^2 + Dy_m.^2 + Dz_0.^2 + (10e-7)^2); 
    C5z = 1 ./ sqrt(Dx_0.^2 + Dy_0.^2 + Dz_p.^2 + (10e-7)^2); 
    C6z = 1 ./ sqrt(Dx_0.^2 + Dy_0.^2 + Dz_m.^2 + (10e-7)^2); 

    % m = (dt/(dx*dy)) * Grad .* mu;  % 2D 
    m = (dt/(dx*dy)) * Grad .* mu; 
    C = 1 + m.*(C1x + C2x + C3y + C4y + C5z + C6z); 

    C1x_2x = C1x.*phi_xp + C2x.*phi_xm; 
    C3y_4y = C3y.*phi_yp + C4y.*phi_ym; 
    C5z_6z = C5z.*phi_zp + C6z.*phi_zm; 

    phi_temp = (1 ./ C) .* (phi + m.*(C1x_2x+C3y_4y) + (dt*Grad).*(lambda1*(u0 - c1).^2) - (dt*Grad).*(lambda2*(u0 - c2).^2)); 
end 

phi = phi_temp; 

phi_visualize = permute(phi,[2,1,3]); 
write_to_binary_file(phi_visualize,n);  % record 

end 

time = toc 

n = n 
T = dt*n 

clear 
clear all 
+6

這不是你不是還在使用32位操作系統或者只有一點點RAM提供了一個特別大的矩陣。如果您向我們展示有問題的代碼,我們可能會提供幫助。更好的是,花一些時間並創建一個[小例子](http://stackoverflow.com/help/mcve),它仍然表明了這個問題。 – horchler

+0

嗨,我已經發布了代碼..它用於3D數據的分割。你有什麼想法,爲什麼它不適合大量的數據? – lee891031

+0

這裏有很多可怕的事情發生。你的代碼不能運行,所以我不能幫忙。我無法看到正在增長或未預先分配的大'for'循環中的任何內容(儘管頂部的'u0_color'不是 - 這有多大?)。你正在運行多少次迭代?代碼中何時何地會陷入停滯?在什麼迭代?開始調試並打印出變量。嘗試使用[profiler]運行(http://www.mathworks.com/help/matlab/ref/profile.html)('profview'是一個GUI界面,但您可能想使用'profile')。 – horchler

回答

0

通常,Matlab會以矩陣的形式跟蹤所有變量。當你處理很多具有許多維度的變量時,RAM存儲器將被分配用於存儲這個變量。因此,在處理大量要運行多次迭代的變量時,最好清除內存中的變量。爲此,請使用命令

clear variable_name_1,variable_name_2,... variable_name_3;

雖然保留所有變量使代碼看起來有組織,但是當您遇到此類問題時,請嘗試清除不需要的變量。

檢查此鏈接詳細使用明確的命令:http://www.mathworks.in/help/matlab/ref/clear.html