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
這不是你不是還在使用32位操作系統或者只有一點點RAM提供了一個特別大的矩陣。如果您向我們展示有問題的代碼,我們可能會提供幫助。更好的是,花一些時間並創建一個[小例子](http://stackoverflow.com/help/mcve),它仍然表明了這個問題。 – horchler
嗨,我已經發布了代碼..它用於3D數據的分割。你有什麼想法,爲什麼它不適合大量的數據? – lee891031
這裏有很多可怕的事情發生。你的代碼不能運行,所以我不能幫忙。我無法看到正在增長或未預先分配的大'for'循環中的任何內容(儘管頂部的'u0_color'不是 - 這有多大?)。你正在運行多少次迭代?代碼中何時何地會陷入停滯?在什麼迭代?開始調試並打印出變量。嘗試使用[profiler]運行(http://www.mathworks.com/help/matlab/ref/profile.html)('profview'是一個GUI界面,但您可能想使用'profile')。 – horchler