OUT OF MEMORY PROBLEM

4 views (last 30 days)
Amit Chakraborty
Amit Chakraborty on 13 Jul 2021
Answered: Matt J on 13 Jul 2021
I have 3 data sets such as , SourceCords(16777216*3), DetectorCords(16777216*3),Voxel(16777216*3). Apart from the Voxel data rest of the 2 datasets were generated during excuting the other line of the codes. Now I want to build a function called attenuation lengths .
[len_ijk, ind_ijk] = ( attenuationLengts(Nx, Ny, Nz, d, dz, bx, by, bz, sourceCords, detectorCords));
Nx, Ny, Nz, d, dz, bx, by, bz= Constant Value.
As you can see that my data size is large it is giving me the error : "Out of memory". At this point I tried clear the unnecessary variables to save some memory but I also got the error : "Out of memory attempting to serialize data for transmission".
I tried another way, I put :::::: tall ( attenuationLengts(Nx, Ny, Nz, d, dz, bx, by, bz, sourceCords, detectorCords));
But still not working as I found later that to use the "Tall Array" youd need to store the data first, but for my case I am not getting my data from the function so I can not store it anyway.
Could anyone can give me some probable solution???????????
  2 Comments
Amit Chakraborty
Amit Chakraborty on 13 Jul 2021
%%%% This is my code actually .
function [length_ijk, index_ijk, source_point, detector_point] = attenuationLengts(Nx, Ny, Nz, d, dz, bx, by, bz, sourceCords, detectorCords)
% IMPROVED_SIDDON_ATTEN
% Input source coordinates
xsource = sourceCords(:,1);
ysource = sourceCords(:,3);
zsource = sourceCords(:,2);
% Input detector coordinates
xdetector = detectorCords(:,1);
ydetector = detectorCords(:,3);
zdetector = detectorCords(:,2);
yy = xsource + d/2;
xx = ysource + d/2;
bxb = bx + double(Nx) * d;
byb = by + double(Ny) * d;
bzb = bz + double(Nz) * dz;
pituus = length(xsource);
NSlices = length(zsource);
length_ijk = cell((pituus),1);
index_ijk = cell((pituus),1);
source_point = cell((pituus),1);
detector_point = cell((pituus),1);
zmax = max(max([zsource;zdetector]));
if zmax==0
zmax = 1;
end
parfor ll = int32(1:pituus)
imax = int32(0);
imin = int32(0);
jmax = int32(0);
jmin = int32(0);
ju = 0;
iu = 0;
ku = 0;
xs = xsource(ll);
xd = xdetector(ll);
ys = ysource(ll);
yd = ydetector(ll);
zs = zsource(ll);
zd = zdetector(ll);
z_loop = int32((zs/zmax)*(NSlices-1)+1);
x_diff = xd-xs;
y_diff = yd-ys;
z_diff = zd-zs;
source_point{ll} = [xs ys zs];
detector_point{ll} = [xd yd zd];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if zs==zd %% Checking if the ray is parallel to Z axis
if abs(y_diff)<1e-8 %% Checking if the ray is parallel to Y axis
if yd <= max(yy) && yd>=min(yy)
tempi = (1:Ny)';
[~,apu] = min(abs(yy(1:Ny)-yd));
tempj = ones(Ny,1,'int32')*(int32(apu));
templ_ijk = ones(Ny,1,'double')*d;
tempk = (z_loop)*ones(size(tempi),'int32');
[tempk, tempi] = sort(((tempj-1).*Ny+tempi+(Nx.*Ny.*(tempk-1))));
index_ijk{ll} = tempk;
length_ijk{ll} = templ_ijk(tempi);
end
elseif abs(x_diff)<1e-8 %% Checking if the ray is parallel to X axis
if xd <= max(xx) && xd>=min(xx)
tempj = (1:Nx)';
[~,apu] = min(abs(xx(1:Nx)-xd));
tempi = ones(Nx,1,'int32')*(int32(apu));
templ_ijk = ones(Nx,1,'double')*d;
tempk = (z_loop)*ones(size(tempj),'int32');
[tempk, tempi] = sort(((tempj-1).*Ny+tempi+(Nx.*Ny.*(tempk-1))));
index_ijk{ll} = tempk;
length_ijk{ll} = templ_ijk(tempi);
end
end
tx0=(bx-xs)/x_diff; % Parametric Values of first X Plane
txback=(bxb-xs)/x_diff; % Parametric Values of last X Plane
ty0 = (by-ys)/y_diff; % Parametric Values of first Y Plane
tyback=(byb-ys)/y_diff; % Parametric Values of last Y Plane
txmin=min(tx0,txback); % minimum parametric value in X plane
txmax=max(tx0,txback); % maximum parametric value in X plane
tymin=min(ty0,tyback); % minimum parametric value in Y plane
tymax=max(ty0,tyback); % maximum parametric value in Y plane
tmin=max(txmin,tymin); % Calculating Overall Minimum parametric value
tmax=min(txmax,tymax); % Calculating Overall Maximum parametric value
if tmin<tmax
if xs<xd
if tmin==txmin
imin=int32(1);
else
pxt=xs+tmin*x_diff;
imin=int32(ceil((pxt-bx)/d));
end
if tmax==txmax
imax=int32(Nx);
else
pxt=xs+tmax*x_diff;
imax=int32(floor((pxt-bx)/d));
end
iu = 1;
tx0=((bx+double(imin)*d)-xs)/x_diff;
elseif xs>xd
if tmin==txmin
imax=int32(Nx-1);
else
pxt=xs+tmin*x_diff;
imax=int32(floor((pxt-bx)/d));
end
if tmax==txmax
imin=int32(0);
else
pxt=xs+tmax*x_diff;
imin=int32(ceil((pxt-bx)/d));
end
iu = -1;
tx0=((bx+double(imax)*d)-xs)/x_diff;
end
if ys<yd
if tmin==tymin
jmin=int32(1);
else
pyt=ys+tmin*y_diff;
jmin=int32(ceil((pyt-by)/d));
end
if tmax==tymax
jmax=int32(Ny);
else
pyt=ys+tmax*y_diff;
jmax=int32(floor((pyt-by)/d));
end
ju = 1;
ty0=((by+double(jmin)*d)-ys)/y_diff;
elseif ys>yd
if tmin==tymin
jmax=int32(Ny-1);
else
pyt=ys+tmin*y_diff;
jmax=int32(floor((pyt-by)/d));
end
if tmax==tymax
jmin=int32(0);
else
pyt=ys+tmax*y_diff;
jmin=int32(ceil((pyt-by)/d));
end
ju = -1;
ty0=((by+double(jmax)*d)-ys)/y_diff;
end
end
Np = (imax - imin + 1) + (jmax - jmin + 1);% Total Number of intersected Planes
L = sqrt(x_diff^2 + y_diff^2); % Distance between the two points
pt = (min(tx0, ty0) + tmin) /2;
tempi = int32(floor(((xs + pt * x_diff) - bx) / d) + 1); % Intersected pixel index in X axis
tempj = int32(floor(((ys + pt * y_diff) - by) / d)); % Intersected pixel index in Y axis
txu = d / abs(x_diff); % Updating Parameters for the pixel index in X axis
tyu = d / abs(y_diff); % Updating Parameters for the pixel index in Y axis
tempk=(z_loop - 1) * Nx * Ny; % Intersected pixel index in Z axis
tc = tmin;
templ_ijk = zeros(Np,1,'double');
idx=zeros(Np,1,'int32');
hpk = 0;
for ii = 1 : Np
if tx0 < ty0
apu = (tx0 - tc);
idx(ii) = (tempj)*Nx+tempi+tempk;
templ_ijk(ii) = apu * L;
tempi = tempi + iu;
tc = tx0;
tx0 = tx0 + txu;
elseif ty0 <= tx0
apu = (ty0 - tc);
idx(ii) = (tempj)*Nx+tempi+tempk;
templ_ijk(ii) = apu * L;
tempj = tempj + ju;
tc = ty0;
ty0 = ty0 + tyu;
end
hpk = hpk + 1;
if tempj < 0 || tempi < 1 || tempj >= Ny || tempi > Nx
break
end
end
templ_ijk = templ_ijk(1:hpk);
idx = idx(1:hpk);
length_ijk{ll} = templ_ijk;
index_ijk{ll} = idx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zs != zd %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
imax=int32(0);
imin=int32(0);
jmax=int32(0);
jmin=int32(0);
kmax=int32(0);
kmin=int32(0);
if abs(y_diff) < 1e-8 %% Checking if the ray is parallel to Y axis
if yd <= max(yy) && yd >= min(yy)
tx0=(bx-xs)/x_diff; % Parametric Values of first X Plane
txback=(bxb-xs)/x_diff; % Parametric Values of last X Plane
tz0=(bz-zs)/z_diff; % Parametric Values of first Z Plane
tzback=(bzb-zs)/z_diff; % Parametric Values of last Z Plane
txmin=min(tx0,txback); % Minimum Parametric Values of X Plane
txmax=max(tx0,txback); % Maximum Parametric Values of X Plane
tzmin=min(tz0,tzback); % Minimum Parametric Values of Z Plane
tzmax=max(tz0,tzback); % Maximum Parametric Values of Z Plane
% Calculating overall minimum and maximum parametric values
% of X and Z plane
tmin=max(txmin,tzmin);
tmax=min(txmax,tzmax);
if tmin<tmax
if xs<xd
if tmin==txmin
imin=int32(1);
else
pxt=xs+tmin*x_diff;
imin=int32(ceil((pxt-bx)/d));
end
if tmax==txmax
imax=int32(Nx);
else
pxt=xs+tmax*x_diff;
imax=int32(floor((pxt-bx)/d));
end
tx0=((bx+double(imin)*d)-xs)/x_diff;
iu = 1;
elseif xs>xd
if tmin==txmin
imax=int32(Nx-1);
else
pxt=xs+tmin*x_diff;
imax=int32(floor((pxt-bx)/d));
end
if tmax==txmax
imin=int32(0);
else
pxt=xs+tmax*x_diff;
imin=int32(ceil((pxt-bx)/d));
end
tx0=((bx+double(imax)*d)-xs)/x_diff;
iu = -1;
end
if zs<zd
if tmin==tzmin
kmin=int32(1);
else
pzt=zs+tmin*z_diff;
kmin=int32(ceil((pzt-bz)/dz));
end
if tmax==tzmax
kmax=int32(Nz);
else
pzt=zs+tmax*z_diff;
kmax=int32(floor((pzt-bz)/dz));
end
tz0=((bz+double(kmin)*dz)-zs)/z_diff;
ku = 1;
elseif zs>zd
if tmin==tzmin
kmax=int32(Nz-1);
else
pzt=zs+tmin*z_diff;
kmax=int32(floor((pzt-bz)/dz));
end
if tmax==tzmax
kmin=int32(0);
else
pzt=zs+tmax*z_diff;
kmin=int32(ceil((pzt-bz)/dz));
end
tz0=((bz+double(kmax)*dz)-zs)/z_diff;
ku = -1;
end
end
Np = (imax - imin + 1) + (kmax - kmin + 1); % Total number of intersected planes
L = sqrt(x_diff^2+z_diff^2); % Distance between two points
pt = (min(tx0, tz0) + tmin) /2;
tempi = int32(floor(((xs + pt * x_diff) - bx) / d) + 1); % Intersected pixel index in X axis
tempk = int32(floor(((zs + pt * z_diff) - bz) / dz)); % Intersected pixel index in Z axis
txu = d / abs(x_diff); % Updating Parameters for the parametric values in X axis
tzu = dz / abs(z_diff); % Updating Parameters for the parametric values in Z axis
tc = tmin;
[~,apu] = min(abs(yy(1:Ny)-yd));
tempj=int32((apu - 1)*Ny);
templ_ijk = zeros(Np,1,'double');
idx = zeros(Np,1,'int32');
hpk = 0;
for ii = 1 : Np
if tx0 < tz0
idx(ii) = (tempj) + tempi + (Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (tx0 - tc) * L; %Intersected voxel length
tempi = tempi + iu; % Updating voxel parammeter in X axis
tc = tx0; % Updating Parametric value
tx0 = tx0 + txu; % Updating Parametric value
else
idx(ii) = (tempj) + tempi + (Nx*Ny*(tempk));%Voxel Index
templ_ijk(ii) = (tz0 - tc) * L; %Intersected voxel length
tempk = tempk + ku; %Updating voxel parammeter in Z axis
tc = tz0; %Updating Parametric value
tz0 = tz0 + tzu; %Updating Parametric value
end
hpk = hpk + 1;
if tempk < 0 || tempi < 1 || tempk >= Nz || tempi > Nx
break
end
end
templ_ijk = templ_ijk(1:hpk);
idx = idx(1:hpk);
length_ijk{ll} = templ_ijk;
index_ijk{ll} = idx;
end
elseif abs(x_diff) < 1e-8 %% Checking if the ray is parallel to X axis
if xd <= max(xx) && xd>=min(xx)
ty0 = (by-ys)/y_diff;
tyback = (byb-ys)/y_diff;
tz0 = (bz-zs)/z_diff;
tzback = (bzb-zs)/z_diff;
tymin = min(ty0,tyback);
tymax = max(ty0,tyback);
tzmin = min(tz0,tzback);
tzmax = max(tz0,tzback);
% Yhteiset arvot
% tmin = max(tymin,tzmin);
tmin = 0; % since source inside voxels
tmax = min(tymax,tzmax);
if tmin < tmax
if ys < yd
if tmin == tymin
jmin = int32(1);
else
pyt = ys + tmin*y_diff;
jmin = int32(ceil((pyt-by)/d));
end
if tmax == tymax % Sama poistumistielle
jmax = int32(Ny);
else
pyt = ys+tmax*y_diff;
jmax = int32(floor((pyt-by)/d));
end
ty0=((by+double(jmin)*d)-ys)/y_diff;
ju = 1;
elseif ys>yd % Jos l���hde on ylemp���n���
if tmin==tymin
jmax=int32(Ny-1);
else
pyt=ys+tmin*y_diff;
jmax=int32(floor((pyt-by)/d));
end
if tmax==tymax
jmin=int32(0);
else
pyt=ys+tmax*y_diff;
jmin=int32(ceil((pyt-by)/d));
end
ty0=((by+double(jmax)*d)-ys)/y_diff;
ju = -1;
end
if zs<zd
if tmin==tzmin
kmin=int32(1);
else
pzt=zs+tmin*z_diff;
kmin=int32(ceil((pzt-bz)/dz));
end
if tmax==tzmax
kmax=int32(Nz);
else
pzt=zs+tmax*z_diff;
kmax=int32(floor((pzt-bz)/dz));
end
tz0=((bz+double(kmin)*dz)-zs)/z_diff;
ku = 1;
elseif zs>zd
if tmin==tzmin
kmax=int32(Nz-1);
else
pzt=zs+tmin*z_diff;
kmax=int32(floor((pzt-bz)/dz));
end
if tmax==tzmax
kmin=int32(0);
else
pzt=zs+tmax*z_diff;
kmin=int32(ceil((pzt-bz)/dz));
end
tz0=((bz+double(kmax)*dz)-zs)/z_diff;
ku = -1;
end
end
Np = (kmax - kmin + 1) + (jmax - jmin + 1); % Number of total intersected planes
L=sqrt(y_diff^2+z_diff^2); % Distance between two points
pt = (min(ty0, tz0) + tmin) /2;
% Intersected
tempj = int32(floor(((ys + pt * y_diff) - by) / d)); % Intersected pixel index in Y axis
tempk = int32(floor(((zs + pt * z_diff) - bz) / dz)); % Intersected pixel index in Z axis
tyu = d / abs(y_diff); % Updating Parameters for the parametric values in Y axis
tzu = dz / abs(z_diff); % Updating Parameters for the parametric values in Z axis
tc = tmin;
[~,apu] = min(abs(xx(1:Nx)-xd));
tempi = int32(apu);
templ_ijk = zeros(Np,1,'double');
idx = zeros(Np,1,'int32');
hpk = 0;
for ii = 1 : Np
if ty0 < tz0
idx(ii) = (tempj)*Ny+tempi+(Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (ty0 - tc) * L; %Intersected voxel length
tempj = tempj + ju; % Updating voxel parammeter in Y axis
tc = ty0; % Updating Parametric value
ty0 = ty0 + tyu; % Updating Parametric value
else
idx(ii) = (tempj)*Ny+tempi+(Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (tz0 - tc) * L; %Intersected voxel length
tempk = tempk + ku; % Updating voxel parammeter in Z axis
tc = tz0; % Updating Parametric value
tz0 = tz0 + tzu; % Updating Parametric value
end
hpk = hpk + 1;
if tempi < 1 || tempj < 0 || tempk < 0 || tempj >= Ny || tempi > Nx || tempk >= Nz
break
end
end
templ_ijk = templ_ijk(1:hpk);
idx = idx(1:hpk);
length_ijk{ll} = templ_ijk;
index_ijk{ll} = idx;
end
end
% If the ray is fully Non-Generic
ty0 = (by-ys)/y_diff; % Parametric Values of first Y Plane
tyback = (byb-ys)/y_diff; % Parametric Values of last Y Plane
tx0 = (bx-xs)/x_diff; % Parametric Values of first X Plane
txback = (bxb-xs)/x_diff; % Parametric Values of first X Plane
tz0 = (bz-zs)/z_diff; % Parametric Values of first Z Plane
tzback = (bzb-zs)/z_diff; % Parametric Values of last Z Plane
txmin = min(tx0,txback); %Minimum parametric values in X plane
txmax = max(tx0,txback); %Maximum parametric values in X plane
tymin = min(ty0,tyback); %Minimum parametric values in Y plane
tymax = max(ty0,tyback); %Maximum parametric values in Y plane
tzmin = min(tz0,tzback); %Minimum parametric values in Z plane
tzmax = max(tz0,tzback); %Maximum parametric values in Z plane
% tmin=max([txmin;tymin;tzmin]); %Overall minimum parametric values
tmin = 0;
tmax=min([txmax;tymax;tzmax]); %Overall maximum parametric values
if tmin < tmax
if xs < xd
if tmin == txmin
imin = int32(1);
else
pxt = xs + tmin*x_diff;
imin = int32(ceil((pxt-bx)/d));
end
if tmax == txmax % Sama poistumistielle
imax = int32(Nx);
else
pxt = xs+tmax*x_diff;
imax = int32(floor((pxt-bx)/d));
end
tx0 = ((bx + double(imin)*d)-xs)/x_diff;
iu = 1;
elseif xs > xd
if tmin == txmin
imax = int32(Nx-1);
else
pxt = xs+tmin*x_diff;
imax = int32(floor((pxt-bx)/d));
end
if tmax == txmax
imin = int32(0);
else
pxt = xs + tmax*x_diff;
imin = int32(ceil((pxt-bx)/d));
end
tx0 = ((bx+double(imax)*d)-xs)/x_diff;
iu = -1;
end
% Same for Y plane
if ys < yd
if tmin == tymin
jmin = int32(1);
else
pyt = ys+tmin*y_diff;
jmin = int32(ceil((pyt-by)/d));
end
if tmax == tymax
jmax = int32(Ny);
else
pyt = ys + tmax*y_diff;
jmax = int32(floor((pyt-by)/d));
end
ty0 = ((by+double(jmin)*d)-ys)/y_diff;
ju = 1;
elseif ys > yd
if tmin == tymin
jmax = int32(Ny-1);
else
pyt = ys + tmin*y_diff;
jmax = int32(floor((pyt-by)/d));
end
if tmax == tymax
jmin = int32(0);
else
pyt = ys + tmax*y_diff;
jmin = int32(ceil((pyt-by)/d));
end
ty0 = ((by+double(jmax)*d)-ys)/y_diff;
ju = -1;
end
% Same for Z plane
if zs < zd
if tmin == tzmin
kmin = int32(1);
else
pzt = zs + tmin*z_diff;
kmin = int32(ceil((pzt-bz)/dz));
end
if tmax == tzmax
kmax = int32(Nz);
else
pzt = zs+tmax*z_diff;
kmax = int32(floor((pzt-bz)/dz));
end
tz0 = ((bz+double(kmin)*dz)-zs)/z_diff;
ku = 1;
elseif zs > zd
if tmin == tzmin
kmax = int32(Nz-1);
else
pzt = zs + tmin*z_diff;
kmax = int32(floor((pzt-bz)/dz));
end
if tmax == tzmax
kmin = int32(0);
else
pzt = zs + tmax*z_diff;
kmin = int32(ceil((pzt-bz)/dz));
end
tz0 = ((bz + double(kmax)*dz) - zs)/z_diff;
ku = -1;
end
end
Np = (imax - imin + 1) + (jmax - jmin + 1) + (kmax - kmin + 1); % Total Number of intersected planes
L = sqrt(x_diff^2 + y_diff^2 + z_diff^2); % Distance between source and detectors
pt = (min(min(ty0, tz0), tx0) + tmin) /2;
tempi = int32(floor(((xs + pt * x_diff) - bx) / d) + 1); % Intersected pixel index in X axis
tempj = int32(floor(((ys + pt * y_diff) - by) / d)); % Intersected pixel index in Y axis
tempk = int32(floor(((zs + pt * z_diff) - bz) / dz)); % Intersected pixel index in Z axis
txu = d / abs(x_diff); % Updating Parameters for the parametric values in X axis
tyu = d / abs(y_diff); % Updating Parameters for the parametric values in Y axis
tzu = dz / abs(z_diff); % Updating Parameters for the parametric values in Z axis
tc = tmin;
templ_ijk = zeros(Np,1,'double');
idx = zeros(Np,1,'int32');
hpk = 0;
for ii = 1 : Np
if ty0 < tz0 && ty0 < tx0
idx(ii) = (tempj)*Ny + tempi + (Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (ty0 - tc) * L; %Intersected voxel length
tempj = tempj + ju; % Updating Parameters for the pixel index in Y axis
tc = ty0; % Updating Parametric value
ty0 = ty0 + tyu; % Updating Parametric value
elseif tz0 < tx0
idx(ii)= (tempj)*Ny+tempi + (Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (tz0 - tc) * L; %Intersected voxel length
tempk = tempk + ku; % Updating Parameters for the pixel index in Z axis
tc = tz0; % Updating Parametric value
tz0 = tz0 + tzu; % Updating Parametric value
else
idx(ii) = (tempj)*Ny + tempi + (Nx*Ny*(tempk)); %Voxel Index
templ_ijk(ii) = (tx0 - tc) * L; %Intersected voxel length
tempi = tempi + iu; % Updating Parameters for the pixel index in X axis
tc = tx0; % Updating Parametric value
tx0 = tx0 + txu; % Updating Parametric value
end
hpk = hpk + 1;
if tempi < 1 || tempj < 0 || tempk < 0 || tempj >= Ny || tempi > Nx || tempk >= Nz
break
end
end
templ_ijk = templ_ijk(1:hpk);
idx = idx(1:hpk);
length_ijk{ll} = templ_ijk;
index_ijk{ll} = idx;
end
end
end

Sign in to comment.

Answers (1)

Matt J
Matt J on 13 Jul 2021
I have the vague impression that you are trying to compute a projection matrix for the purpose of tomographic forward and back projection. It is unrealistic to do that for realistic system sizes. You should consider existing forward/back projection algorithm libraries, e.g.,

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!