+++ /dev/null
-loriane
-19 avril
\ No newline at end of file
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## Analytical_1D
+
+## Author: Loriane Weber <lweber@gpid16a-1802>
+## 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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## CTFPropagation 1D
+
+## Author: Loriane Weber <lweber@gpid16a-1801>
+## 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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## CTFPropagation_2D
+
+## Author: Loriane Weber <lweber@gpid16a-1802>
+## 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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## PhantTibo_Analytical
+
+## Author: Loriane Weber <lweber@rnice7-0102>
+## 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
--- /dev/null
+%%%%%%%%%%%%%%%%
+%%% 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)
+
--- /dev/null
+
+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
--- /dev/null
+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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## FresnelTransform
+
+## Author: Loriane Weber <lweber@hpc2-0601>
+## 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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## FresnelTransform
+
+## Author: Loriane Weber <lweber@hpc2-0601>
+## 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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## PhantTibo_Analytical
+
+## Author: Loriane Weber <lweber@rnice7-0102>
+## 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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## csquare1
+
+## Author: Loriane Weber <lweber@rnice7-0103>
+## 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
--- /dev/null
+%# 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
--- /dev/null
+% ## 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 <tpikonen@pcu.helsinki.fi>
+% ## Created: 5.5.2000
+% ## Keywords: padding image processing
+% ## 2006-11-28 P. Cloetens <cloetens@esrf.fr>
+% ## * 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
--- /dev/null
+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);
--- /dev/null
+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);
--- /dev/null
+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);
--- /dev/null
+%# 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
--- /dev/null
+%# 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
projects / CreaPhase.git / commitdiff