From d0401c49b6b7511cfdaa0534b78bd3c5b2b0637a Mon Sep 17 00:00:00 2001 From: Loriane Weber Date: Wed, 20 Apr 2016 10:49:24 +0200 Subject: [PATCH] Useful functions for simulations (created by LW, free to use) --- test_lw.txt | 2 - utilities_LW/Analytical_1D.m | 38 ++++++ utilities_LW/CTFPropagation_1D.m | 65 ++++++++++ utilities_LW/CTFPropagation_2D.m | 82 +++++++++++++ utilities_LW/CircleMat_Analytical.m | 63 ++++++++++ utilities_LW/CreateSphere3D.m | 41 +++++++ utilities_LW/FrequencySpace.m | 9 ++ utilities_LW/FrequencySpace1.m | 8 ++ utilities_LW/FresnelTransform_1D.m | 53 ++++++++ utilities_LW/FresnelTransform_2D.m | 68 +++++++++++ utilities_LW/PhantTibo_Analytical.m | 99 +++++++++++++++ utilities_LW/astra_iradon_3d.m | 32 +++++ utilities_LW/astra_radon_3d.m | 37 ++++++ utilities_LW/attenuation.edf | Bin 0 -> 5760512 bytes utilities_LW/csquare1.m | 40 ++++++ utilities_LW/csquare2.m | 34 ++++++ utilities_LW/delta_beta_map_1200.edf | Bin 0 -> 5760512 bytes utilities_LW/im_pad.m | 175 +++++++++++++++++++++++++++ utilities_LW/mha_read_header.m | 96 +++++++++++++++ utilities_LW/mha_read_slice.m | 93 ++++++++++++++ utilities_LW/mha_read_volume.m | 88 ++++++++++++++ utilities_LW/ssquare1.m | 28 +++++ utilities_LW/ssquare2.m | 34 ++++++ 23 files changed, 1183 insertions(+), 2 deletions(-) delete mode 100644 test_lw.txt create mode 100644 utilities_LW/Analytical_1D.m create mode 100644 utilities_LW/CTFPropagation_1D.m create mode 100644 utilities_LW/CTFPropagation_2D.m create mode 100644 utilities_LW/CircleMat_Analytical.m create mode 100644 utilities_LW/CreateSphere3D.m create mode 100644 utilities_LW/FrequencySpace.m create mode 100644 utilities_LW/FrequencySpace1.m create mode 100644 utilities_LW/FresnelTransform_1D.m create mode 100644 utilities_LW/FresnelTransform_2D.m create mode 100644 utilities_LW/PhantTibo_Analytical.m create mode 100644 utilities_LW/astra_iradon_3d.m create mode 100644 utilities_LW/astra_radon_3d.m create mode 100644 utilities_LW/attenuation.edf create mode 100644 utilities_LW/csquare1.m create mode 100644 utilities_LW/csquare2.m create mode 100644 utilities_LW/delta_beta_map_1200.edf create mode 100644 utilities_LW/im_pad.m create mode 100644 utilities_LW/mha_read_header.m create mode 100644 utilities_LW/mha_read_slice.m create mode 100644 utilities_LW/mha_read_volume.m create mode 100644 utilities_LW/ssquare1.m create mode 100644 utilities_LW/ssquare2.m 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 0000000000000000000000000000000000000000..26823adeef2ece5fa4538922df4b0533c17f0e6f GIT binary patch literal 5760512 zcmeF&&x%|{6bIm}^AuUQRG0}wK@tBDC9VW@F9Qy7KrzH@qWA>9gAd`u_&iRhy+mei zcURq8b?S8YH!$3p+g0bB`YNj*liwdcyT1D2`sG(oAAS7ji)SBRY`!kHkBk2vk3M?% z)z4SYuWwge-v0dg>sQxbzr5XZ`}?z-=Rabuu3HcP+!>6fcF zuWo++{&KVR@x{ev#qF<;-~RgT&2O8VLPk+{a zAIZ^-Nq_(W0t5&UAV7cs0RjXF5Fl`b0<~FH$@(K?FfsuG1PBlyK!5-N0t5&UAV7dX zB7u?&s@0lA$2JN90t5&UAV7cs0RjXF5FkK+KyrarnbeXx$!Uy7fB*pk1PBlyK!5-N z0t5&UATUs1N#4|J&p-yP2@oJafB*pk1PBlyK!5-N0tAu?^va-G)J;ZZ90CLg5FkK+ z009C72oNAZfB=DDffgB4?;61*+7Tc?fB*pk1PBlyK!5-N0t5&o5a^vb^{AhK%?Jbt z5FkK+009C72oNAZfB*pk0Rla;rKp+#JX#PSK!5-N0t5&UAV7cs0RjXF3>1jUnR@ma z$fq>{0t5&UAV7cs0RjXF5FkK+Ku3X|nNoD^j!a4uAV7cs0RjXF5FkK+009C72n-O2 z&X*$k4Pexg009C72oNAZfB*pk1PBlyK%j#_M3ywFM+Yt?2@oJafB*pk1PBlyK!5-N z0t5mDM&(M8eFHhQB0zuu0RjXF5FkK+009C72oPu}5Sby3?$wY@DFOrt5FkK+009C7 z2oNAZfB=CYfzg>#Wd9&iZ3qw`K!5-N0t5&UAV7cs0RjXX2}EW`S$Z|1Q-%Nm0t5&U zAV7cs0RjXF5FkLHvp|+SX~Z0zX_Y5HfB*pk1PBlyK!5-N0t5&UXdo~mH_Fnp0iO~C z2oNAZfB*pk1PBlyK!5-N0v!dig76Jqa5FkK+009C72oNAZfB=EU0@*U7O!G9RR*nDx z0t5&UAV7cs0RjXF5FkKcfk39bDBCOxIBg+7fB*pk1PBlyK!5-N0t5&UXe^K|Gs-ki zV`}9H5FkK+009C72oNAZfB*pk1S$nGWkp$MsU)?W009C72oNAZfB*pk1PBlyK%k*O z*32l=Tn)LEB0zuu0RjXF5FkK+009C72oR_g$dnUhou!i0asmVh5FkK+009C72oNAZ zfB=Dp0$DSoOmj8lR*C=t0t5&UAV7cs0RjXF5FkLHQXo@Kly#O$Qp*VtAV7cs0RjXF z5FkK+009C78VY31j55vDkXtDN1PBlyK!5-N0t5&UAV7csfl7f)IZ@VGDoHITK!5-N z0t5&UAV7cs0RjXF5NIfnH8aXIS3_>42oNAZfB*pk1PBlyK!5-N0t6}rGUY^BXQ?E$ zoB#m=1PBlyK!5-N0t5&UAV8p@K-SDC(_9U?l_EfZ009C72oNAZfB*pk1PBnQ6v&hl zWu2vx)N%p@2oNAZfB*pk1PBlyK!5;&h5}hLqfB!(a-S5FkK+009C72oNAZfB*pk1R4uu&5bh5)tFj20t5&UAV7cs0RjXF z5FkK+0D%Pp8FHemGc4e=g#ZBp1PBlyK!5-N0t5&UAV8q8K-TOi!(5H2l_Nla009C7 z2oNAZfB*pk1PBmVAdn#^$~waWPFn~NAV7cs0RjXF5FkK+009C7ItXOVk4Dbbfmulc z1PBlyK!5-N0t5&UAV7csfrSDibE2%h7gE|rfB*pk1PBlyK!5-N0t5&UAkaY|Ylbv( zt`5vf5+Fc;009C72oNAZfB*pk1PH7U7?~Ai>%D@}RssYF5FkK+009C72oNAZfB=C` z0@*U8k@IvSR+az(0t5&UAV7cs0RjXF5FkKcrNGF{C`<2^gtikPK!5-N0t5&UAV7cs z0RjXFbQH*vA&s1)Bkx(4{^W;m{weuBYb1v!K!5-N0t5&UAV7cs0RjYy1V&~@qk9(t zSyUuYl2w(gDQadB0RjXF5FkK+009C72oNA}Rsy3lq>=s4>Y_Uwwj_^Qv!=tyN)jMI zfB*pk1PBlyK!5-N0!Jz^GD8~G`$(@U#n@{ysU>?;9NG{B2oNAZfB*pk1PBlyK!Ct- zfl)cqh`z(Ik9dG3`BRHMM@(jP0t5&UAV7cs0RjXF5FkLHk-&%?X;jZfWM@^TMF!QP z##x~|ECB)p2oNAZfB*pk1PBlya8F=Vo)p>l-l4rXD9|EzYFXppz!nf7K!5-N0t5&U zAV7cs0Rks25Sb-K_d4l|>vXu5xl_wposO$40RjXF5FkK+009C72oNA}Bm&X7QbfNa zp-(kt%e<*atyEZsAwYlt0RjXF5FkK+009C7!UZBSrKld^cqh=XN8Z$<<_SV0zGr4p0(>wY6k%V1PBlyK!5-N0t5&UAaF_oy>q4>^-sw^>MT9;rHI;5 zPWliaK!5-N0t5&UAV7cs0RnjidgM&4Yv$FM1lHtGwR=y^ZuEH~v!#e$qdoN}K!5-N0t5&UAV7cs0RjZ_3)JROCF}E_ z`ZyQ6I zDWYd^Vw3zk=as2@1-vG5g2e; zXULW!W+<_+h5!Kq1PBlyK!5-N0t5&UAV47QZ?C2-?2vgxc9fyVAqSrmAV7cs0RjXF z5FkK+009C72*hVsB^g$RY$;*}|NaONAV7cs0RjXF5FkK+009C74hu|qQuK4*WD(g> zh93J2d`5r(0RjXF5FkK+009C72oNAJ?r*QAtZdgEWXP5xX4rLp?+FkfK!5-N0t5&U zAV7cs0RjZJ0#lAO>bdJ4BC?|lJ$9YndjbRq5FkK+009C72oNAZfB=EmwL`HcV_Ph5FkK+009C7 z2oNAZfB*pk1hxWMa;s^MxqIoE9c8FJcV_Ph5FkK+009C72oNAZfB*pk1hxWGZj|LY zcP|+-rk*p*o!L7A1PBlyK!5-N0t5&UAV7csfvrH6+-jO*?p}IkM;U6*o!L7A1PBly zK!5-N0t5&UAV7csfvv!l8)bR!x|a+&Q;!*To!5H;1PBlyK!5-N0t5&UAV7csfvrH6 z+-jO**S++}jxyBTbzbiY5FkK+009C72oNAZfB*pk1hxWGZj|M@?_M(GO)Y2GcUGSf zAV7cs0RjXF5FkK+009C72y6wic9fykzO(v_009C72oNAZfB*pk1PBly zK;W&wlpT$JzMW|Prwo}>%NgcR=`8^Q1PBlyK!5-N0t5&UAV7e?U4hZL)zts)Tywv- z%#TLanme0!1PBlyK!5-N0t5&UAV7cs0Rr;^Q+_o1IX~U~w7QOGE$QuF#2oNAZ zfB*pk1PBlyK!5;&LjuQ=U2V?MB10Nk0z%k`l zo9Rd9PiuPLm-mJM0RjXF5FkK+009C72oNAZpj2Q~hE+P>-F0hnq!D}X&hK{u1PBly zK!5-N0t5&UAV7csfjWU>%da-48<9gT=~;IZRudpVfB*pk1PBlyK!5-N0t5)`5{S;Q zcFj@qeo3Abxu@nV))F8SBh9$>TMkX0t5&UAV7cs0RjXF5FkK+K)pck z46A;kvKHd<8smGeq`K%*A zfB*pk1PBlyK!5-N0t5&UXeZDr!)iBCS*0HN)V`W!Q&~lT009C72oNAZfB*pk1PBly z&`V%xj@4_P(z^RHrxu@;&S4z^0t5&UAV7cs0RjXF5FkK+K%Bs#EGur3x=tYwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA Yz<>b*1`HT5V8DO@0|pEjFks*?5C>CD<^TWy literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..0c941d5db4ed85177434015b56bb10096712873e GIT binary patch literal 5760512 zcmeF&!HQf(6b9f~=PB|6fk_|+6gL_}vT_mkGT;yc21Cpyir~tZ@u7Sbr_){%W^Q*^ z-CK3)boVzf+?m@|=bZW~tDngq_nu#0{dE2A#nby=-~aCUmzSHbhug>H|Bw4$-Fxxo z>ecnliibBpfAit}_4n^?_T2pb{PnA!H@`o7{qDE-H=jMddVlr9hqs%h&wly!>f@`| zZ+?8Z+4|`6aW{0FM>~2oNAZfB*pk1PBlyK!5-N0s{r2a;Bbr2J&f5 zfB*pk1PBlyK!5-N0t5&UAka~uXQmWgyCaj*1PBlyK!5-N0t5&UAV7cs0RjUAqVuJQ zeghb_BtU=w0RjXF5FkK+009C72oUHX5RoN~>d}EqNdg225FkK+009C72oNAZfB=C& zfl;|qWZytetq2exK!5-N0t5&UAV7cs0RjXX3Pff|qkA=EQ;Glq0t5&UAV7cs0RjXF z5FkJxNMLlP6xlzBR2u>W2oNAZfB*pk1PBlyK!5;&Mgo!9QI=kf=#(KqfB*pk1PBly zK!5-N0t5&U=q!*WPZ}{tXIkY65FkK+009C72oNAZfB*pk1R4m8$c?h}Y`~`k0RjXF z5FkK+009C72oNAZfIvrqEIHE1IXdzxO@IIa0t5&UAV7cs0RjXF5FoHpU}SETrT0og z+X)aLK!5-N0t5&UAV7cs0RjX%31rESGR)D5SXlxD2oNAZfB*pk1PBlyK!5;&6#^M@ zqbxJ5V6>G00RjXF5FkK+009C72oNAZpo2h`+$hr=9hj9QK!5-N0t5&UAV7cs0RjXF z5LhUXDL2Y8%R)-q2oNAZfB*pk1PBlyK!5-N0t7k;WXX&&&C!8bNdg225FkK+009C7 z2oNAZfB=C70+}+SY_lxjw1ofx0t5&UAV7cs0RjXF5FkLHu|T%WDAPQRsg)x@fB*pk z1PBlyK!5-N0t5&USRjxoFUmH{0!~{95FkK+009C72oNAZfB*pk1R4uu%ZxJ3)0kR0 z0t5&UAV7cs0RjXF5FkK+0D($@Oj%LZSt>~_CqRGz0RjXF5FkK+009C72oPu}kTo;P zG*?4zr3erpK!5-N0t5&UAV7cs0RjXn1v2GCS!byvwVVI}0t5&UAV7cs0RjXF5FkLH zp+MHmDAQaGxs@V7fB*pk1PBlyK!5-N0t5&Us1(SQ6J?#HlGJhn1PBlyK!5-N0t5&U zAV7csfrbKEGows%HRM)`009C72oNAZfB*pk1PBlyK%i0}Q%;n1mP%5~2@oJafB*pk z1PBlyK!5-N0t6ZgWX+5+&DD@wDFOrt5FkK+009C72oNAZfB=C?flN73)>$e^Ehj*L z009C72oNAZfB*pk1PBmlD3CQX$~0F)ZlwqiAV7cs0RjXF5FkK+009C7Dg`p-L|JF4 zB(kJDxZ6QE_009C72oNAZfB*pk1PBmlERZ!j$}m@BYUKzJAV7cs0RjXF z5FkK+009C776@d>iL%bHfYTNN1PBlyK!5-N0t5&UAV7csfer#$^P`b-bzoMK009C7 z2oNAZfB*pk1PBlyKwzQ3$ebu^?}e1M5g*?O;Fw3Pq>0t5&UAV7cs0RjXF5FkLHlR&l% zY2-Yeh?ONkfB*pk1PBlyK!5-N0t5&USSc_vGs@C?C86yE2oNAZfB*pk1PBlyK!5-N z0v!diWJn|D=*WB4rJw$LaZ&Pp)<_OdfB*pk1PBlyK!5-N0t5&Y35?8+M)xiPvZzR) zB&#Y}Q`F2N0t5&UAV7cs0RjXF5FkL{tOQ19NF)27)kSwWY)KxqW=)5Yl_Wrb009C7 z2oNAZfB*pk1ddc-WQH`V_mN&zim}&ZQcL!xIJ6-M5FkK+009C72oNAZfB=Ev0;6)I z5q*baAMpT7@~0Mij+o5o1PBlyK!5-N0t5&UAV7dXBY_b)(x{$|$j+)viwvqojk7{` zSONqH5FkK+009C72oNAZ;EuqkJSnp8okM$ZP@qNb)Uw9Gfh{0FfB*pk1PBlyK!5-N z0t8N4ATmpe?sd`^*XeLAbElTIIvrP80t5&UAV7cs0RjXF5FkL{NCcvDrHFnM$+@lK#&kG!cz%@d$@1OfyI5FkK+ z009C72oNAZV1+YtK-)LDAwOA)oBob(|; zfB*pk1PBlyK!5-N0tE63^vIc7*UYPzf3k>-sb`P;!X_j@fB*pk1PBlyK!5-N0t6xi zTIWqI>qVUZRD1Nyl_G1O>iswu0RjXF5FkK+009C72oNApFVHe~T3f4rdOJQBh|HHF zdVPL2KNBE8fB*pk1PBlyK!5-N0tD7(PiywmZzB+qDMj{ZqoN7{0t5&UAV7cs0RjXF z5Fl`B0&8-p+P$Y{H~Ku0*-}KW(VqGfAV7cs0RjXF5FkK+009E|1#0uClJ)sdeVh|U zU1*W{I*R#fVWVRI1yT-y=0t5&UAV7cs0RjXF z5FkK+0D<;@bG0i2E3uBqk23Tqv9N{!0RjXF5FkK+009C72oNAZU@NdIZ)*L%#5+T_ z6fr}Ig*5~S5FkK+009C72oNAZfB*pkaesR?WnqWRBeJ6mJq|heoB#m=1PBlyK!5-N z0t5&UAV45KyDG`BGGt2;Gx+yMfB*pk1PBlyK!5-N0t5&UAaGb<%9Em>`zDLXjxzMv zXW%me1PBlyK!5-N0t5&UAV7csfpLF(HDzVH?jS?96fwiD^LtN#009C72oNAZfB*pk z1PBlyuoakcq*2dZ_YjdCW$3Z%{N58FK!5-N0t5&UAV7cs0RjXFjQ`uKDKne9iwxOP z#0+z1_l^Jo0t5&UAV7cs0RjXF5FkKcD=_6pqn~s45s@8b=rMP8?+6eeK!5-N0t5&U zAV7cs0RjZJ0;6-QssG%4WXP5xW|%v>cLWF!AV7cs0RjXF5FkK+009Dd{_o(f3@ai# z%Fx5VMFIo}5FkK+009C72oNAZfB=Dcfhjl2@|>G8L%!5=hPg9)M}PnU0t5&UAV7cs z0RjXF5FoG>$dX%4bIjdK&+I5e?YT31M}PnU0t5&UAV7cs0RjXF5FoG>m~x{m&$)Za zkTLa~VeZV{5g5Gds#qd+yBM5gm~x{m&wcljA#Z9q!@jfni~s=w z1PBlyK!5-N0t5&UAV6R%kR`X8=Gb>HEwiHxwf3FWX9NfkAV7cs0RjXF5FkK+009D@ z1g7k0^z+k1^FL+CoLbH>e@br&5FkK+009C72oNAZfB*pk1a1qA&aI~Yx96Jsy=8th zvew+$ydyw>009C72oNAZfB*pk1PBnA7nt&+(a-tm?!FzFJGJP2_e9wn<1PBlyK!5-N0t5&UAV7cs0Rp80qcW`0`EIXUlOv7TdwYJr6Cgl<009C7 z2oNAZfB*pk1PIg#99w?1Io*gHYDv$!o3NSy0RjXF5FkK+009C72oNAZV3$C2hP7*s zn)gfcq{uxrXR($50RjXF5FkK+009C72oNAZ;IKechIM$7@)ePJR86n)6Ie-r009C7 z2oNAZfB*pk1PBlyP%3bY`PJrZHMvs6+EQ=p2oNAZfB*pk1PBlyK!5-N0tD&>dS_Vm z6P4|V$fZhpluc_D0RjXF5FkK+009C72oNAZfWTUT);ZSNc}n+}WK2EQl+I@z0RjXF z5FkK+009C72oNAZfIvHeRvA{iiOMSV$fx$zESt(I0t5&UAV7cs0RjXF5FkK+0D)ct zOLMGV^OV-zmpQfgtaJ|R2oNAZfB*pk1PBlyK!5-N0tDg&4rN(!lhk!;kymYM_#!}n z009C72oNAZfB*pk1PBly5Ft>SXGP4S2LS>E2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5Fqeh zJ9I1n0000m$lrPt`wAfg1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r z3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@ z0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VK zfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5 zV8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM z7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b* z1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd z0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwA zz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEj zFkrxd0RsjM7%*VKfB^#r3>YwAz<>b*1`HT5V8DO@0|pEjFkrxd0RsjM7%*VKfB^#r T3>YwAz<>b*1`HT5VBi@Dc|%Vj literal 0 HcmV?d00001 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 -- 2.45.0