0 Stimmen

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 Kommentare
Keine anzeigen Keine ausblenden

Matt J
Matt J am 13 Jul. 2021
Bearbeitet: Matt J am 13 Jul. 2021
The part of the code you've shown does not reveal where the out of memory is coming from. The size of the input data sets alone is unlikely to be the problem. The dimensions you've indicated correspond to arrays of about 300 MB in double floats or 150 MB in single floats. Not really that large...
SourceCoords=rand(16777216,3);
whos SourceCoords
Name Size Kilobytes Class Attributes
SourceCoords 16777216x3 393216 double
Amit Chakraborty
Amit Chakraborty am 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

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Matt J
Matt J am 13 Jul. 2021

0 Stimmen

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.,
https://github.com/CERN/TIGRE#features

Kategorien

Mehr zu Spline Construction finden Sie in Hilfe-Center und File Exchange

Tags

Gefragt:

Amit Chakraborty
Amit Chakraborty
am 13 Jul. 2021

Beantwortet:

Matt J
Matt J
am 13 Jul. 2021

Community Treasure Hunt

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

Start Hunting!

Translated by