From: Loriane Weber Date: Wed, 20 Apr 2016 08:49:24 +0000 (+0200) Subject: Useful functions for simulations (created by LW, free to use) X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=commitdiff_plain;h=d0401c49b6b7511cfdaa0534b78bd3c5b2b0637a Useful functions for simulations (created by LW, free to use) --- diff --git a/test_lw.txt b/test_lw.txt deleted file mode 100644 index 262f9ae..0000000 --- a/test_lw.txt +++ /dev/null @@ -1,2 +0,0 @@ -loriane -19 avril \ No newline at end of file diff --git a/utilities_LW/Analytical_1D.m b/utilities_LW/Analytical_1D.m new file mode 100644 index 0000000..0db8f94 --- /dev/null +++ b/utilities_LW/Analytical_1D.m @@ -0,0 +1,38 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## Analytical_1D + +## Author: Loriane Weber +## Created: 2015-12-04 + +%% delta : delta corresponding to the wire material +%% u : abscisse, in um +%% R : radius of the wire, in um +%% u0 : wire center +%% delta : coeff delta + +function [ Phi ] = Analytical_1D(delta, u, R, u0) + +if(isscalar(delta) && isscalar(u0) && isscalar(R)) + Phi=zeros(size(u)); + Phi=2*delta*real(sqrt(R^2-(u-u0).^2)); +else + disp('delta, radius and center position should be scalars') +end + + +endfunction diff --git a/utilities_LW/CTFPropagation_1D.m b/utilities_LW/CTFPropagation_1D.m new file mode 100644 index 0000000..2801aab --- /dev/null +++ b/utilities_LW/CTFPropagation_1D.m @@ -0,0 +1,65 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## CTFPropagation 1D + +## Author: Loriane Weber +## Created: 2015-11-20 + +function [ Id ] = CTFPropagation_1D (ph, B, D, lambda, oversamp, ftot) + +[mp, np]=size(ph); + +if(np!=1) + disp('error') + return +end + +%% propagators in 1D (sin and cos) +ssq=2*ssquare1(lambda*D, ftot, 1)'; +csq=2*csquare1(lambda*D, ftot, 1)'; + +tmp_B=zeros(2*mp, 1); +tmp_ph=zeros(2*mp, 1); +tmp_B(1:mp,1)=B; +tmp_ph(1:mp,1)=ph; + +% fourier domain +tmp_B=fft(tmp_B); +tmp_ph=fft(tmp_ph); + +tmp_Id=-csq.*tmp_B + ssq.*tmp_ph; +tmp_Id=1+real(ifft(tmp_Id)); + + +Id=tmp_Id(1:mp, 1); + + if oversamp==2 + ipf = [0.5 1 0.5]; % modelise le binning du detecteur + Idd = conv(Id,ipf,'same')./conv(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + Id=Idd; + elseif oversamp==4 + ipf = [0.5 1.0 1.0 1.0 0.5]; % modelise le binning du detecteur + Idd = conv(Id,ipf,'same')./conv(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + + + Id=Idd; + end % oversamp + + +endfunction diff --git a/utilities_LW/CTFPropagation_2D.m b/utilities_LW/CTFPropagation_2D.m new file mode 100644 index 0000000..f307e63 --- /dev/null +++ b/utilities_LW/CTFPropagation_2D.m @@ -0,0 +1,82 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## CTFPropagation_2D + +## Author: Loriane Weber +## Created: 2015-12-11 + +function [ IdCTF ] = CTFPropagation_2D (ph, B, D, lambda, oversamp, ftot, gtot) + + +[mp, np]=size(ph); +[mp_a, np_a] = size(B); + +% check dimensions + if (mp_a != mp || np_a != np) + disp('error in the dimensions !! Phi and B should have the same dimension') + return; + end + +% propagators in 2D (sin and cos) +ssq=2*ssquare2(lambda*D, ftot', gtot', 1); +csq=2*csquare2(lambda*D, ftot', gtot', 1); +%keyboard +% initialisation +tmp_B=zeros(2*mp, 2*np); +tmp_Phi=zeros(2*mp, 2*np); + +%tmp_B(1:mp,1:np)=B; +%tmp_ph(1:mp,1:np)=ph; + +tmp_B=im_pad(B, 2*np, 2*mp, "zeros"); +tmp_Phi=im_pad(ph, 2*np, 2*mp, "zeros"); + +tmp_B=fft2(tmp_B); +tmp_Phi=fft2(tmp_Phi); + +tmp_Id=-csq.*tmp_B + ssq.*tmp_Phi; +tmp_Id=1+real(ifft2(tmp_Id)); + +IdCTF=cut(tmp_Id, mp, np); + + if oversamp==2 + + ipf = [ 0.25 0.5 0.25; % modelise le binning du detecteur + 0.50 1.0 0.50; + 0.25 0.5 0.25]; + Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end); + IdCTF=Idd; + + elseif oversamp==4 + ipf = [0.25 0.50 0.50 0.50 0.25; + 0.50 1.00 1.00 1.00 0.50; + 0.50 1.00 1.00 1.00 0.50; + 0.50 1.00 1.00 1.00 0.50; + 0.25 0.50 0.50 0.50 0.25]; + + Idd = conv2(IdCTF,ipf,'same')./conv2(ones(size(IdCTF)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end); + IdCTF=Idd; + + end % end oversampling + + + + + +endfunction diff --git a/utilities_LW/CircleMat_Analytical.m b/utilities_LW/CircleMat_Analytical.m new file mode 100644 index 0000000..1f28b72 --- /dev/null +++ b/utilities_LW/CircleMat_Analytical.m @@ -0,0 +1,63 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## PhantTibo_Analytical + +## Author: Loriane Weber +## Created: 2015-12-10 + +%%%%%% INPUT PARAMETERS +%% ps : is the physical pixel size, in um +%% RedAttFact : reduced attenuation factor +%% angles : angles of projections +%% save : 1 or 0 (1 if you want to save in edf, 0 sinon + +function [MuP_tot, Phi_tot, B_tot] = CircleMat_Analytical (angles, ps, R_circle, circle_center, absorption_mat, beta_mat, delta_mat, taille, save) + +du= 1; % echantillonage de la projection +angles_rad=(angles*pi/180) + pi/2; % in radians % +pi/2 to be coherent with radon transform + +circle_center=circle_center - floor(taille/2); +u=[-(taille-1)/2:du:+(taille-1)/2]; +disp('proj') + + +for i=1:length(angles) + + u0_mat= cos(angles_rad(i))*circle_center(1) + sin(angles_rad(i))*circle_center(2); + +## Proj attenuation + MuP_mat=Analytical_1D(absorption_mat, u, R_circle, u0_mat); + MuP_tot(:, i) = MuP_mat * ps; + +## Proj Phi + Phi_mat=Analytical_1D(delta_mat, u, R_circle, u0_mat); + Phi_tot(:, i) = Phi_mat * ps; + +## Proj Beta + B_mat=Analytical_1D(beta_mat, u, R_circle, u0_mat); + B_tot(:, i) = B_mat *ps; + +end + +if save + edfwrite('Sino_AnalyticalCirclePET_Mu.edf', MuP_tot, 'float32') + edfwrite('Sino_AnalyticalCirclePET_Phi.edf', Phi_tot, 'float32') + edfwrite('Sino_AnalyticalCirclePET_B.edf', B_tot, 'float32') +end + + +endfunction diff --git a/utilities_LW/CreateSphere3D.m b/utilities_LW/CreateSphere3D.m new file mode 100644 index 0000000..32c0bce --- /dev/null +++ b/utilities_LW/CreateSphere3D.m @@ -0,0 +1,41 @@ +%%%%%%%%%%%%%%%% +%%% Function that creates a 3D matrix (size dimx, dimy, dimz) containing a sphere, centered at (cx, +%%% cy, cz), of radius r +%%% Loriane, oct 2014 +%c=128 +%r=40 +%%%%%%%%%%%%%%%% + + +r=70 +cx=200; +cy=cx; +cz=50; +dimx=400; +dimy=dimx; +dimz=100; +coeff=9.917*10^-7; % mu ou autre + +Sphere=zeros(dimx, dimy, dimz); + +for(i=1:dimx) + for(j=1:dimy) + for(k=1:dimz) + if((i-cx)^2+(j-cy)^2+(k-cz)^2 <= r^2) + Sphere(i,j,k)=1; + end + end + end +end + + +Sphere=Sphere*coeff; +name='sphere_phase_Mg_19keV' + +fid=fopen(strcat(name, '.raw'), 'w') +fwrite(fid, Sphere, 'float32'); +fclose(fid); + +cmd = ['/data/id19/bones01/bones/program/linux_debian/bin/rawtomhd filein=', strcat(name, '.raw'), ' fileout=', strcat(name, '.mhd'), ' dimx=', num2str(dimx), ' dimy=', num2str(dimy), ' dimz=', num2str(dimz), ' format=FLOAT']; +unix(cmd) + diff --git a/utilities_LW/FrequencySpace.m b/utilities_LW/FrequencySpace.m new file mode 100644 index 0000000..00308eb --- /dev/null +++ b/utilities_LW/FrequencySpace.m @@ -0,0 +1,9 @@ + +function [ftot, gtot]=FrequencySpace(mf, nf, pixelsize) + +fsamplex = 1/pixelsize; +fsampley= 1/pixelsize; +%[ftot,gtot] = meshgrid([0:fsamplex/mf:fsamplex/2 -fsamplex/2+fsamplex/mf:fsamplex/mf:-fsamplex/mf],[0:fsampley/nf:fsampley/2 -fsampley/2+fsampley/nf:fsampley/nf:-fsampley/nf]); %originale version +[ftot,gtot] = meshgrid([0:fsamplex/mf:fsamplex/2 -fsamplex/2+fsamplex/mf:fsamplex/mf:-fsamplex/mf], [0:fsampley/nf:fsampley/2 -fsampley/2+fsampley/nf:fsampley/nf:-fsampley/nf]); + +end diff --git a/utilities_LW/FrequencySpace1.m b/utilities_LW/FrequencySpace1.m new file mode 100644 index 0000000..ae9e02a --- /dev/null +++ b/utilities_LW/FrequencySpace1.m @@ -0,0 +1,8 @@ +function [ftot]=FrequencySpace1(mf, pixelsize) + +% /!\ pixel size in meter, to get rid of "le" +fsamplex = 1/(pixelsize); + +%[ftot,gtot] = meshgrid([0:fsamplex/mf:fsamplex/2 -fsamplex/2+fsamplex/mf:fsamplex/mf:-fsamplex/mf],[0:fsampley/nf:fsampley/2 -fsampley/2+fsampley/nf:fsampley/nf:-fsampley/nf]); %originale version +ftot = meshgrid([0:fsamplex/mf:fsamplex/2 -fsamplex/2+fsamplex/mf:fsamplex/mf:-fsamplex/mf], 1); +end diff --git a/utilities_LW/FresnelTransform_1D.m b/utilities_LW/FresnelTransform_1D.m new file mode 100644 index 0000000..6a5d8ed --- /dev/null +++ b/utilities_LW/FresnelTransform_1D.m @@ -0,0 +1,53 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## FresnelTransform + +## Author: Loriane Weber +## Created: 2015-07-01 + + +function [FrIm] = FresnelTransform_1D(ph, B, D, lambda, oversamp, ftot) + +[m, n]=size(ph); +u = exp(-B).*exp(+i*ph); +ftot=ftot'; + +Pd = exp(-i*pi*lambda*D*(ftot.^2)); + +ud = ones(size(ftot, 1), 1); +ud(1:m,1) = u; + +ud= ifft( fft(ud).*Pd ); +ud=ud(1:m, 1); +FrIm = abs(ud).^2; +%keyboard + if oversamp==2 + ipf = [0.5 1 0.5]; + Idd = conv2(FrIm,ipf,'same')./conv2(ones(size(FrIm)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + FrIm=Idd; + elseif oversamp==4 + ipf = [0.5 1 1 1 0.5 ]; + Idd = conv2(FrIm,ipf,'same')./conv2(ones(size(FrIm)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, :); + % smooth the data + %ipf = [0.5 1 0.5]; % modelise le binning du detecteur + %Idd = conv(Idd,ipf,'same')./conv(ones(size(Idd)),ipf,'same'); + FrIm=Idd; + end +%keyboard +endfunction diff --git a/utilities_LW/FresnelTransform_2D.m b/utilities_LW/FresnelTransform_2D.m new file mode 100644 index 0000000..770ec46 --- /dev/null +++ b/utilities_LW/FresnelTransform_2D.m @@ -0,0 +1,68 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## FresnelTransform + +## Author: Loriane Weber +## Created: 2015-07-01 + + +function Id = FresnelTransform_2D(ph, B, D, lambda, oversamp, ftot, gtot) + +[m, n]=size(ph); +[m_a, n_a] = size(B); + +% check dimensions +if (m_a != m || n_a != n) + disp('error in the dimensions !! Phi and B should have the same dimension') +return; +end + +u = exp(-B).*exp(+i.*ph); + +% propagator for the distance D +Pd = exp(-i*pi*lambda*D*(ftot.^2+gtot.^2)); % ifftshift +ud = ones(2*m, 2*n); +ud(1:m,1:n) = u; + +ud=fft2(ud); + +ud = ud.*(Pd)'; +ud = ifft2(ud); + +ud = ud(1:m,1:n); +Id = abs(ud).^2; + + if oversamp==2 + ipf = [0.25 0.5 0.25; + 0.50 1.0 0.50; + 0.25 0.5 0.25]; + Idd = conv2(Id,ipf,'same')./conv2(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end); + Id=Idd; + elseif oversamp ==4 + + ipf = [0.25 0.5 0.5 0.5 0.25; + 0.50 1.0 1.0 1.0 0.50; + 0.50 1.0 1.0 1.0 0.50; + 0.50 1.0 1.0 1.0 0.50; + 0.25 0.5 0.5 0.5 0.25]; + Idd = conv2(Id,ipf,'same')./conv2(ones(size(Id)),ipf,'same'); + Idd = Idd(oversamp:oversamp:end, oversamp:oversamp:end); + Id=Idd; + end % end binning detector + +endfunction diff --git a/utilities_LW/PhantTibo_Analytical.m b/utilities_LW/PhantTibo_Analytical.m new file mode 100644 index 0000000..b603e62 --- /dev/null +++ b/utilities_LW/PhantTibo_Analytical.m @@ -0,0 +1,99 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## PhantTibo_Analytical + +## Author: Loriane Weber +## Created: 2015-12-10 + +%%%%%% INPUT PARAMETERS +%% ps : is the physical pixel size, in um +%% RedAttFact : reduced attenuation factor +%% angles : angles of projections +%% save : 1 or 0 (1 if you want to save in edf, 0 sinon + +function [MuP_tot, Phi_tot, B_tot ] = PhantTibo_Analytical (angles, ps, save) + +disp('params'); + +RedAttFact=1; + +mu_Al=10.7765 * 100/RedAttFact; % *100 to be in m-1 +mu_Mg= 5.56697 * 100/RedAttFact; +mu_PET=0.89889 * 100/RedAttFact; + +delta_Al=15.03*10^-7/RedAttFact; +delta_Mg=9.917*10^-7/RedAttFact; +delta_PET=8.265*10^-7/RedAttFact; + +beta_Al=5.59*10^-9/RedAttFact; +beta_Mg=2.89*10^-9/RedAttFact; +beta_PET=0.46*10^-9/RedAttFact; + +R_Al=73/2; +R_Mg=37/2; +R_PET=58/2; + +cen_Al=[607-600; 585-600]; % 600 to put the center at the center of the image (1200*1200) +cen_Mg=[679-600; 640-600]; +cen_PET=[730-600; 630-600]; + +du= 1; % echantillonage de la projection + +angles_rad=(angles*pi/180) + pi/2; % in radians % +pi/2 to be coherent with radon transform + +dim=1701; % taille de la proj, wrt radon function +u=[-(dim-1)/2:du:+(dim-1)/2]; +disp('proj') + + +for i=1:length(angles) + + u0_Al= cos(angles_rad(i))*cen_Al(1) + sin(angles_rad(i))*cen_Al(2); + u0_Mg= cos(angles_rad(i))*cen_Mg(1) + sin(angles_rad(i))*cen_Mg(2); + u0_PET= cos(angles_rad(i))*cen_PET(1) + sin(angles_rad(i))*cen_PET(2); + +## Proj attenuation + MuP_Al=Analytical_1D(mu_Al, u, R_Al, u0_Al); + MuP_Mg=Analytical_1D(mu_Mg, u, R_Mg, u0_Mg); + MuP_PET=Analytical_1D(mu_PET, u, R_PET, u0_PET); + + MuP_tot(:, i) = ( MuP_Al+ MuP_Mg + MuP_PET ) * ps; + +## Proj Phi + Phi_Al=Analytical_1D(delta_Al, u, R_Al, u0_Al); + Phi_Mg=Analytical_1D(delta_Mg, u, R_Mg, u0_Mg); + Phi_PET=Analytical_1D(delta_PET, u, R_PET, u0_PET); + + Phi_tot(:, i) = ( Phi_Al+ Phi_Mg + Phi_PET ) * ps; + +## Proj Beta + B_Al=Analytical_1D(beta_Al, u, R_Al, u0_Al); + B_Mg=Analytical_1D(beta_Mg, u, R_Mg, u0_Mg); + B_PET=Analytical_1D(beta_PET, u, R_PET, u0_PET); + + B_tot(:, i) = ( B_Al+ B_Mg + B_PET ) *ps; + +end + +if save + edfwrite('Sino_Analytical_Mu.edf', MuP_tot, 'float32') + edfwrite('Sino_Analytical_Phi.edf', Phi_tot, 'float32') + edfwrite('Sino_Analytical_B.edf', B_tot, 'float32') +end + + +endfunction diff --git a/utilities_LW/astra_iradon_3d.m b/utilities_LW/astra_iradon_3d.m new file mode 100644 index 0000000..ba74957 --- /dev/null +++ b/utilities_LW/astra_iradon_3d.m @@ -0,0 +1,32 @@ +function vol = astra_iradon_3d(sino, theta) +% this function calculates 3D iradon transform of the volume using ASTRA tools + + sz = size(sino); + + % initialize ASTRA projector: + proj_geom = astra_create_proj_geom('parallel3d', 1, 1, sz(3), sz(1), theta); + vol_geom = astra_create_vol_geom(sz(1), sz(1), sz(3)); + %proj_id = astra_create_projector('linear3d', proj_geom, vol_geom); + + % copy data to GPU: + volume_id = astra_mex_data3d('create','-vol', vol_geom, 0); + sino_id = astra_mex_data3d('create','-sino', proj_geom, sino); + + % initialize algorithm: + cfg = astra_struct('BP3D_CUDA'); + %cfg.ProjectorId = proj_id; + cfg.ProjectionDataId = sino_id; + cfg.ReconstructionDataId = volume_id; + %cfg.FilterType = 'none'; + + % run: + alg_id = astra_mex_algorithm('create', cfg); + astra_mex_algorithm('run', alg_id); + vol = astra_mex_data3d('get', volume_id) * 2 * size(proj_geom.ProjectionAngles,2) / pi; + + % kill stuff: + astra_mex_data3d('delete', sino_id); + astra_mex_data3d('delete', volume_id); + astra_mex_algorithm('delete', alg_id); + %astra_mex_projector('delete', proj_id); +end \ No newline at end of file diff --git a/utilities_LW/astra_radon_3d.m b/utilities_LW/astra_radon_3d.m new file mode 100644 index 0000000..5574a68 --- /dev/null +++ b/utilities_LW/astra_radon_3d.m @@ -0,0 +1,37 @@ +function sino = astra_radon_3d(volume, theta) +% this function calculates 3D radon transform of the volume using ASTRA tools + + sz = size(volume); + + % initialize ASTRA geometries: + proj_geom = astra_create_proj_geom('parallel3d', 1, 1, sz(3), sz(1), theta); + vol_geom = astra_create_vol_geom(sz(1), sz(2), sz(3)); + %proj_id = astra_create_projector('linear3d', proj_geom, vol_geom); + %proj_id = astra_create_projector('cuda3d', proj_geom, vol_geom); + + % copy data to GPU: + volume_id = astra_mex_data3d('create','-vol', vol_geom, volume); + sino_id = astra_mex_data3d('create','-sino', proj_geom, 0); + + % initialize algorithm: + cfg = astra_struct('FP3D_CUDA'); + %cfg = astra_struct('FP_CUDA'); + %cfg.ProjectorId = proj_id; + %cfg.ProjectionGeometry = proj_geom; + %cfg.VolumeGeometry = vol_geom; + cfg.ProjectionDataId = sino_id; + cfg.VolumeDataId = volume_id; + + % run algorithm: + alg_id = astra_mex_algorithm('create', cfg); + astra_mex_algorithm('iterate', alg_id); + + % get data: + sino = astra_mex_data3d('get',sino_id); + + % kill stuff: + astra_mex_data3d('delete', sino_id); + astra_mex_data3d('delete', volume_id); + astra_mex_algorithm('delete', alg_id); + % astra_mex_projector('delete', proj_id); +end \ No newline at end of file diff --git a/utilities_LW/attenuation.edf b/utilities_LW/attenuation.edf new file mode 100644 index 0000000..26823ad Binary files /dev/null and b/utilities_LW/attenuation.edf differ diff --git a/utilities_LW/csquare1.m b/utilities_LW/csquare1.m new file mode 100644 index 0000000..4783def --- /dev/null +++ b/utilities_LW/csquare1.m @@ -0,0 +1,40 @@ +## Copyright (C) 2015 Loriane Weber +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## csquare1 + +## Author: Loriane Weber +## Created: 2015-06-04 + +function csq=csquare1(betah,f,varargin) + +f_isall = 0; +switch nargin +case 3 + f_isall = varargin{1}; +end + +csq=cos((pi*betah)*f.^2); + +if (f_isall~=1) + [n m]=size(f); + n=n-1; + m=m-1; + csq=cat(2,cat(1,csq,flipdim(csq(2:n,:,:),1)),cat(1,flipdim(csq(:,2:m,:),2),flipdim(flipdim(csq(2:n,2:m,:),1),2))); +end + + +end diff --git a/utilities_LW/csquare2.m b/utilities_LW/csquare2.m new file mode 100644 index 0000000..7cd951e --- /dev/null +++ b/utilities_LW/csquare2.m @@ -0,0 +1,34 @@ +%# csq=csquare2(beta,f,g) returns cos(pi*beta*(f^2+g^2)) +%# beta is equal to lambda*D/l^2, i.e. the Fresnelnumber +%# for the complete spectrum (from 0 to fsample) +%# f_isall = 0 +%# only f and g from 0 to fsample/2 should be given, the other part is +%# calculated using the periodicity in frequency-space +%# f_isall = 1 +%# complete frequency range is given +%# allows for astigmatism (5th argument) + +function csq=csquare2(betah,f,g,varargin) + +f_isall = 0; +switch nargin +case 3 + betav = betah; +case 4 + f_isall = varargin{1}; + betav = betah; +case 5 + f_isall = varargin{1}; + betav = varargin{2}; +end + +csq=cos((pi*betah)*f.^2+(pi*betav)*g.^2); + +if (f_isall~=1) + [n m]=size(f); + n=n-1; + m=m-1; + csq=cat(2,cat(1,csq,flipdim(csq(2:n,:,:),1)),cat(1,flipdim(csq(:,2:m,:),2),flipdim(flipdim(csq(2:n,2:m,:),1),2))); +end + +end diff --git a/utilities_LW/delta_beta_map_1200.edf b/utilities_LW/delta_beta_map_1200.edf new file mode 100644 index 0000000..0c941d5 Binary files /dev/null and b/utilities_LW/delta_beta_map_1200.edf differ diff --git a/utilities_LW/im_pad.m b/utilities_LW/im_pad.m new file mode 100644 index 0000000..a460fef --- /dev/null +++ b/utilities_LW/im_pad.m @@ -0,0 +1,175 @@ +% ## Copyright (C) 2000 Teemu Ikonen +% ## +% ## This program is free software; you can redistribute it and/or +% ## modify it under the terms of the GNU General Public License +% ## as published by the Free Software Foundation; either version 2 +% ## of the License, or (at your option) any later version. +% ## +% ## This program is distributed in the hope that it will be useful, but +% ## WITHOUT ANY WARRANTY; without even the implied warranty of +% ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% ## General Public License for more details. +% ## +% ## You should have received a copy of the GNU General Public License +% ## along with this program; if not, write to the Free Software +% ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +% + +% ## -*- texinfo -*- +% ## @deftypefn {Function File} {} im_pad(@var{A}, @var{xpad}, @var{ypad}, [@var{padding}, [@var{const}]]) +% ## Pad (augment) a matrix for application of image processing algorithms. +% ## +% ## Pads the input image @var{A} with @var{xpad}(1) elements from left, +% ## @var{xpad}(2), elements from right, @var{ypad}(1) elements from above +% ## and @var{ypad}(2) elements from below. +% ## Values of padding elements are determined from the optional arguments +% ## @var{padding} and @var{const}. @var{padding} is one of +% ## +% ## @table @samp +% ## @item "zeros" +% ## pad with zeros (default) +% ## +% ## @item "ones" +% ## pad with ones +% ## +% ## @item "constant" +% ## pad with a value obtained from the optional fifth argument const +% ## +% ## @item "symmetric" +% ## pad with values obtained from @var{A} so that the padded image mirrors +% ## @var{A} starting from edges of @var{A} +% ## +% ## @item "reflect" +% ## same as symmetric, but the edge rows and columns are not used in the padding +% ## +% ## @item "replicate" +% ## pad with values obtained from A so that the padded image +% ## repeates itself in two dimensions +% ## +% ## @end table +% ## @end deftypefn +% +% ## Author: Teemu Ikonen +% ## Created: 5.5.2000 +% ## Keywords: padding image processing +% ## 2006-11-28 P. Cloetens +% ## * add padding possibility 'extend' +% ## * rename to im_pad +% ## 2006-12-03 PC +% ## * extend with average of last value over extend_avg rows/columns +% +% ## A nice test matrix for padding: +% ## A = 10*[1:5]' * ones(1,5) + ones(5,1)*[1:5] + +function retval = im_pad(A, xpad, ypad, padding, const) + +try empty_list_elements_ok_save = empty_list_elements_ok; +catch empty_list_elements_ok_save = 0; +end +try warn_empty_list_elements_save = warn_empty_list_elements; +catch warn_empty_list_elements_save = 0; +end +unwind_protect + +if nargin < 4 + padding = 'zeros'; +end +if nargin < 5 + const = 1; +end + +extend_avg = const; %# we use the same fifth argument + +origx = size(A,2); +origy = size(A,1); + +if isscalar(xpad) + retx = xpad; + xpad(1) = ceil((retx-origx)/2); + xpad(2) = retx-origx-xpad(1); +else + retx = origx + xpad(1) + xpad(2); +end + +if isscalar(ypad) + rety = ypad; + ypad(1) = ceil((rety-origy)/2); + ypad(2) = rety-origy-ypad(1); +else + rety = origy + ypad(1) + ypad(2); +end + +empty_list_elements_ok = 1; +warn_empty_list_elements = 0; + +if(strcmp(padding, 'zeros')) + retval = zeros(rety,retx); + retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A; + elseif(strcmp(padding,'ones')) + retval = ones(rety,retx); + retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A; + elseif(strcmp(padding,'constant')) + retval = const.*ones(rety,retx); + retval(ypad(1)+1 : ypad(1)+origy, xpad(1)+1 : xpad(1)+origx) = A; + elseif(strcmp(padding,'replicate')) + y1 = origy-ypad(1)+1; + x1 = origx-xpad(1)+1; + if(y1 < 1 || x1 < 1 || ypad(2) > origy || xpad(2) > origx) + error('Too large padding for this padding type'); + else + yrange1 = y1 : origy; + yrange2 = 1 : ypad(2); + xrange1 = x1 : origx; + xrange2 = 1 : xpad(2); + retval = [ A(yrange1, xrange1), A(yrange1, :), A(yrange1, xrange2); + A(:, xrange1), A, A(:, xrange2); + A(yrange2, xrange1), A(yrange2, :), A(yrange2, xrange2) ]; + end + elseif(strcmp(padding,'symmetric')) + y2 = origy-ypad(2)+1; + x2 = origx-xpad(2)+1; + if(ypad(1) > origy || xpad(1) > origx || y2 < 1 || x2 < 1) + error('Too large padding for this padding type'); + else + yrange1 = 1 : ypad(1); + yrange2 = y2 : origy; + xrange1 = 1 : xpad(1); + xrange2 = x2 : origx; + retval = [ fliplr(flipud(A(yrange1, xrange1))), flipud(A(yrange1, :)), fliplr(flipud(A(yrange1, xrange2))); + fliplr(A(:, xrange1)), A, fliplr(A(:, xrange2)); + fliplr(flipud(A(yrange2, xrange1))), flipud(A(yrange2, :)), fliplr(flipud(A(yrange2, xrange2))) ]; + end + elseif(strcmp(padding,'reflect')) + y2 = origy-ypad(2); + x2 = origx-xpad(2); + if(ypad(1)+1 > origy || xpad(1)+1 > origx || y2 < 1 || x2 < 1) + error('Too large padding for this padding type'); + else + yrange1 = 2 : ypad(1)+1; + yrange2 = y2 : origy-1; + xrange1 = 2 : xpad(1)+1; + xrange2 = x2 : origx-1; + retval = [ fliplr(flipud(A(yrange1, xrange1))), flipud(A(yrange1, :)), fliplr(flipud(A(yrange1, xrange2))); + fliplr(A(:, xrange1)), A, fliplr(A(:, xrange2)); + fliplr(flipud(A(yrange2, xrange1))), flipud(A(yrange2, :)), fliplr(flipud(A(yrange2, xrange2))) ]; + end + elseif(strcmp(padding,'extend')) + if (extend_avg > 1) + retval = [ repmat(mean2(A(1:extend_avg,1:extend_avg)),ypad(1),xpad(1)) , repmat(mean(A(1:extend_avg,:),1),ypad(1),1) , repmat(mean2(A(1:extend_avg,end-extend_avg+1:end)),ypad(1),xpad(2)) ; + repmat(mean(A(:,1:extend_avg),2),1,xpad(1)) , A , repmat(mean(A(:,end-extend_avg+1:end),2),1,xpad(2)) ; + repmat(mean2(A(end-extend_avg+1:end,1:extend_avg)),ypad(2),xpad(1)) , repmat(mean(A(end-extend_avg+1:end,:),1),ypad(2),1) , repmat(mean2(A(end-extend_avg+1:end,end-extend_avg+1:end)),ypad(2),xpad(2)) ]; + else + retval = [ repmat(A(1,1),ypad(1),xpad(1)) , repmat(A(1,:),ypad(1),1) , repmat(A(1,end),ypad(1),xpad(2)) ; + repmat(A(:,1),1,xpad(1)) , A , repmat(A(:,end),1,xpad(2)) ; + repmat(A(end,1),ypad(2),xpad(1)) , repmat(A(end,:),ypad(2),1) , repmat(A(end,end),ypad(2),xpad(2)) ]; + end + else + error('Unknown padding type'); + end + +unwind_protect_cleanup + empty_list_elements_ok = empty_list_elements_ok_save; + warn_empty_list_elements = warn_empty_list_elements_save; +end_unwind_protect + + end \ No newline at end of file diff --git a/utilities_LW/mha_read_header.m b/utilities_LW/mha_read_header.m new file mode 100644 index 0000000..21fd4b7 --- /dev/null +++ b/utilities_LW/mha_read_header.m @@ -0,0 +1,96 @@ +function info =mha_read_header(filename) +% Function for reading the header of a Insight Meta-Image (.mha,.mhd) file +% +% info = mha_read_header(filename); +% +% examples: +% 1, info=mha_read_header() +% 2, info=mha_read_header('volume.mha'); + +if(exist('filename','var')==0) + [filename, pathname] = uigetfile('*.mha', 'Read mha-file'); + filename = [pathname filename]; +end + +fid=fopen(filename,'rb'); +if(fid<0) + fprintf('could not open file %s\n',filename); + return +end + +info.Filename=filename; +info.Format='MHA'; +info.CompressedData='false'; +readelementdatafile=false; +while(~readelementdatafile) + str=fgetl(fid); + s=find(str=='=',1,'first'); + if(~isempty(s)) + type=str(1:s-1); + data=str(s+1:end); + while(type(end)==' '); type=type(1:end-1); end + while(data(1)==' '); data=data(2:end); end + else + type=''; data=str; + end + + switch(lower(type)) + case 'ndims' + info.NumberOfDimensions=sscanf(data, '%d')'; + case 'dimsize' + info.Dimensions=sscanf(data, '%d')'; + case 'elementspacing' + info.PixelDimensions=sscanf(data, '%lf')'; + case 'elementsize' + info.ElementSize=sscanf(data, '%lf')'; + if(~isfield(info,'PixelDimensions')) + info.PixelDimensions=info.ElementSize; + end + case 'elementbyteordermsb' + info.ByteOrder=lower(data); + case 'anatomicalorientation' + info.AnatomicalOrientation=data; + case 'centerofrotation' + info.CenterOfRotation=sscanf(data, '%lf')'; + case 'offset' + info.Offset=sscanf(data, '%lf')'; + case 'binarydata' + info.BinaryData=lower(data); + case 'compresseddatasize' + info.CompressedDataSize=sscanf(data, '%d')'; + case 'objecttype', + info.ObjectType=lower(data); + case 'transformmatrix' + info.TransformMatrix=sscanf(data, '%lf')'; + case 'compresseddata'; + info.CompressedData=lower(data); + case 'binarydatabyteordermsb' + info.ByteOrder=lower(data); + case 'elementdatafile' + info.DataFile=data; + readelementdatafile=true; + case 'elementtype' + info.DataType=lower(data(5:end)); + case 'headersize' + val=sscanf(data, '%d')'; + if(val(1)>0), info.HeaderSize=val(1); end + otherwise + info.(type)=data; + end +end + +switch(info.DataType) + case 'char', info.BitDepth=8; + case 'uchar', info.BitDepth=8; + case 'short', info.BitDepth=16; + case 'ushort', info.BitDepth=16; + case 'int', info.BitDepth=32; + case 'uint', info.BitDepth=32; + case 'float', info.BitDepth=32; + case 'double', info.BitDepth=64; + otherwise, info.BitDepth=0; +end +if(~isfield(info,'HeaderSize')) + info.HeaderSize=ftell(fid); +end +fclose(fid); diff --git a/utilities_LW/mha_read_slice.m b/utilities_LW/mha_read_slice.m new file mode 100644 index 0000000..e7b69c0 --- /dev/null +++ b/utilities_LW/mha_read_slice.m @@ -0,0 +1,93 @@ +function S = mha_read_slice(info, num_slice) +% Function for reading the volume of a Insight Meta-Image (.mha, .mhd) file +% +% volume = tk_read_volume(file-header) +% +% examples: +% 1: info = mha_read_header() +% V = mha_read_volume(info); +% imshow(squeeze(V(:,:,round(end/2))),[]); +% +% 2: V = mha_read_volume('test.mha'); + +if(~isstruct(info)), info=mha_read_header(info); end + + +switch(lower(info.DataFile)) + case 'local' + otherwise + % Seperate file + info.Filename=fullfile(fileparts(info.Filename),info.DataFile); +end + +% Open file +switch(info.ByteOrder(1)) + case ('true') + fid=fopen(info.Filename,'rb','ieee-be'); + otherwise + fid=fopen(info.Filename,'rb','ieee-le'); +end + +switch(lower(info.DataFile)) + case 'local' + % Skip header + fseek(fid,info.HeaderSize,'bof'); + otherwise + fseek(fid,0,'bof'); +end + + +datasize=prod(info.Dimensions)*info.BitDepth/8; % in bytes +slicesize=info.Dimensions(1)*info.Dimensions(2)*info.BitDepth/8; +OffsetToSelectedSlice=(num_slice-1)*slicesize; +fseek(fid, OffsetToSelectedSlice, 'cof'); + + +switch(info.CompressedData(1)) + case 'f' + % Read the Data + switch(info.DataType) + case 'char' + S = int8(fread(fid,slicesize,'char')); + case 'uchar' + S = uint8(fread(fid,slicesize,'uchar')); + case 'short' + S = int16(fread(fid,slicesize,'short')); + case 'ushort' + S = uint16(fread(fid,slicesize,'ushort')); + case 'int' + S = int32(fread(fid,slicesize,'int')); + case 'uint' + S = uint32(fread(fid,slicesize,'uint')); + case 'float' + S = single(fread(fid,info.Dimensions(1)*info.Dimensions(2),'float')); + case 'double' + S = double(fread(fid,slicesize,'double')); + end + case 't' + switch(info.DataType) + case 'char', DataType='int8'; + case 'uchar', DataType='uint8'; + case 'short', DataType='int16'; + case 'ushort', DataType='uint16'; + case 'int', DataType='int32'; + case 'uint', DataType='uint32'; + case 'float', DataType='single'; + case 'double', DataType='double'; + end + Z = fread(fid,inf,'uchar=>uint8'); + S = zlib_decompress(Z,DataType); + +end + +fclose(fid); +S = reshape(S,[info.Dimensions(1), info.Dimensions(2)]); + +function M = zlib_decompress(Z,DataType) +import com.mathworks.mlwidgets.io.InterruptibleStreamCopier +a=java.io.ByteArrayInputStream(Z); +b=java.util.zip.InflaterInputStream(a); +isc = InterruptibleStreamCopier.getInterruptibleStreamCopier; +c = java.io.ByteArrayOutputStream; +isc.copyStream(b,c); +M=typecast(c.toByteArray,DataType); diff --git a/utilities_LW/mha_read_volume.m b/utilities_LW/mha_read_volume.m new file mode 100644 index 0000000..e6d4a44 --- /dev/null +++ b/utilities_LW/mha_read_volume.m @@ -0,0 +1,88 @@ +function V = mha_read_volume(info) +% Function for reading the volume of a Insight Meta-Image (.mha, .mhd) file +% +% volume = tk_read_volume(file-header) +% +% examples: +% 1: info = mha_read_header() +% V = mha_read_volume(info); +% imshow(squeeze(V(:,:,round(end/2))),[]); +% +% 2: V = mha_read_volume('test.mha'); + +if(~isstruct(info)), info=mha_read_header(info); end + + +switch(lower(info.DataFile)) + case 'local' + otherwise + % Seperate file + info.Filename=fullfile(fileparts(info.Filename),info.DataFile); +end + +% Open file +switch(info.ByteOrder(1)) + case ('true') + fid=fopen(info.Filename,'rb','ieee-be'); + otherwise + fid=fopen(info.Filename,'rb','ieee-le'); +end + +switch(lower(info.DataFile)) + case 'local' + % Skip header + fseek(fid,info.HeaderSize,'bof'); + otherwise + fseek(fid,0,'bof'); +end + +datasize=prod(info.Dimensions)*info.BitDepth/8; % in bytes + +switch(info.CompressedData(1)) + case 'f' + % Read the Data + switch(info.DataType) + case 'char' + V = int8(fread(fid,datasize,'char')); + case 'uchar' + V = uint8(fread(fid,datasize,'uchar')); + case 'short' + V = int16(fread(fid,datasize,'short')); + case 'ushort' + V = uint16(fread(fid,datasize,'ushort')); + case 'int' + V = int32(fread(fid,datasize,'int')); + case 'uint' + V = uint32(fread(fid,datasize,'uint')); + case 'float' + V = single(fread(fid,datasize,'float')); + case 'double' + V = double(fread(fid,datasize,'double')); + end + case 't' + switch(info.DataType) + case 'char', DataType='int8'; + case 'uchar', DataType='uint8'; + case 'short', DataType='int16'; + case 'ushort', DataType='uint16'; + case 'int', DataType='int32'; + case 'uint', DataType='uint32'; + case 'float', DataType='single'; + case 'double', DataType='double'; + end + Z = fread(fid,inf,'uchar=>uint8'); + V = zlib_decompress(Z,DataType); + +end + +fclose(fid); +V = reshape(V,info.Dimensions); + +function M = zlib_decompress(Z,DataType) +import com.mathworks.mlwidgets.io.InterruptibleStreamCopier +a=java.io.ByteArrayInputStream(Z); +b=java.util.zip.InflaterInputStream(a); +isc = InterruptibleStreamCopier.getInterruptibleStreamCopier; +c = java.io.ByteArrayOutputStream; +isc.copyStream(b,c); +M=typecast(c.toByteArray,DataType); diff --git a/utilities_LW/ssquare1.m b/utilities_LW/ssquare1.m new file mode 100644 index 0000000..9b7a44b --- /dev/null +++ b/utilities_LW/ssquare1.m @@ -0,0 +1,28 @@ +%# ssq=ssquare2(beta,f,g) returns sin(pi*beta*(f^2+g^2)) +%# beta is equal to lambda*D/l^2, i.e. the Fresnelnumber +%# for the complete spectrum (from 0 to fsample) +%# f_isall = 0 +%# only f and g from 0 to fsample/2 should be given, the other part is +%# calculated using the periodicity in frequency-space +%# f_isall = 1 +%# complete frequency range is given +%# allows for astigmatism (5th argument) + +function ssq=ssquare1(betah,f,varargin) + +f_isall = 0; +switch nargin +case 3 + f_isall = varargin{1}; +end + +ssq=sin((pi*betah)*f.^2); + +if (f_isall~=1) + [n m]=size(f); + n=n-1; + m=m-1; + ssq=cat(2,cat(1,ssq,flipdim(ssq(2:n,:,:),1)),cat(1,flipdim(ssq(:,2:m,:),2),flipdim(flipdim(ssq(2:n,2:m,:),1),2))); +end + +end diff --git a/utilities_LW/ssquare2.m b/utilities_LW/ssquare2.m new file mode 100644 index 0000000..8e4758b --- /dev/null +++ b/utilities_LW/ssquare2.m @@ -0,0 +1,34 @@ +%# ssq=ssquare2(beta,f,g) returns sin(pi*beta*(f^2+g^2)) +%# beta is equal to lambda*D/l^2, i.e. the Fresnelnumber +%# for the complete spectrum (from 0 to fsample) +%# f_isall = 0 +%# only f and g from 0 to fsample/2 should be given, the other part is +%# calculated using the periodicity in frequency-space +%# f_isall = 1 +%# complete frequency range is given +%# allows for astigmatism (5th argument) + +function ssq=ssquare2(betah,f,g,varargin) + +f_isall = 0; +switch nargin +case 3 + betav = betah; +case 4 + f_isall = varargin{1}; + betav = betah; +case 5 + f_isall = varargin{1}; + betav = varargin{2}; +end + +ssq=sin((pi*betah)*f.^2+(pi*betav)*g.^2); + +if (f_isall~=1) + [n m]=size(f); + n=n-1; + m=m-1; + ssq=cat(2,cat(1,ssq,flipdim(ssq(2:n,:,:),1)),cat(1,flipdim(ssq(:,2:m,:),2),flipdim(flipdim(ssq(2:n,2:m,:),1),2))); +end + +end \ No newline at end of file