Remove thorax and add head

CERN · Nov 30, 2016 · a72140d · a72140d
1 parent 210adc6
commit a72140d
Show file tree

Hide file tree

Showing 20 changed files with 73 additions and 84 deletions.
diff --git a/Demos/d02_SampleData.m b/Demos/d02_SampleData.m
@@ -74,12 +74,12 @@
 %--------------------------------------------------------------------------
 
 %--------------------------------------------------------------------------
-% Thorax phantom
+% head phantom
 %
 %
-thorax=thoraxPhantom(geo.nVoxel); %default is 128^3
+head=headPhantom(geo.nVoxel); %default is 128^3
 % show it
-plotImg(thorax,'Dim','Z');
-
+plotImg(head,'Dim','Z');
+citeme('headPhantom')
 
 
diff --git a/Demos/d03_generateData.m b/Demos/d03_generateData.m
@@ -48,12 +48,12 @@
 % define projection angles (in radians)
 angles=linspace(0,2*pi,100);
 % load phatnom image
-thorax=thoraxPhantom(geo.nVoxel);
+head=headPhantom(geo.nVoxel);
 
 % Simulate forward projection.
 % Strongly suggested to use 'iterpolated' option for more accurate
 % projections. reduce geo.accuracy for better results
-projections=Ax(thorax,geo,angles,'interpolated');
+projections=Ax(head,geo,angles,'interpolated');
 
 % Add realistic noise. Adds photon scattering noise ('Poisson') and
 % electronic noise of the detector ('Gaussian').

diff --git a/Demos/d04_SimpleReconstruction.m b/Demos/d04_SimpleReconstruction.m
@@ -50,9 +50,9 @@
 % define angles
 angles=linspace(0,2*pi,100);
 % Load thorax phatom data
-thorax=thoraxPhantom(geo.nVoxel);
+head=headPhantom(geo.nVoxel);
 % generate projections
-projections=Ax(thorax,geo,angles,'interpolated');
+projections=Ax(head,geo,angles,'interpolated');
 % add noise
 noise_projections=addCTnoise(projections);
 

diff --git a/Demos/d05_plottingFunctions.m b/Demos/d05_plottingFunctions.m
@@ -48,8 +48,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 %% Reconstruct image using OS-SART and FDK
 imgFDK=FDK(noise_projections,geo,angles);
@@ -149,6 +149,6 @@
 plotImg(imgFDK,'Dim',dimension,'Step',step,'CLims',clims,'Colormap',colormap,'Savegif',giffilename);
 
 % Remember: You can always plot more than 1 image together!
-plotImg([thorax imgFDK imgOSSART],'Dim','z')
+plotImg([head imgFDK imgOSSART],'Dim','z')
 % Or even the errors!
-plotImg([abs(thorax-imgFDK) abs(thorax-imgOSSART)],'Dim',3);
+plotImg([abs(head-imgFDK) abs(head-imgOSSART)],'Dim',3);
diff --git a/Demos/d06_Algorithms01.m b/Demos/d06_Algorithms01.m
@@ -47,8 +47,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,200);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 
 %% Usage of FDK
@@ -80,4 +80,4 @@
 
 % but it can be seen that one has bigger errors in the whole image, while
 % hte other just in the boundaries
-plotImg([abs(thorax-imgFDK1) abs(thorax-imgFDK2)],'Dim','Z');
+plotImg([abs(head-imgFDK1) abs(head-imgFDK2)],'Dim','Z');
diff --git a/Demos/d07_Algorithms02.m b/Demos/d07_Algorithms02.m
@@ -48,8 +48,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 %% SART family of algorithms
 %
@@ -118,9 +118,9 @@
 
 % SIRT and SART both have no extra input parameters.
 % =========================================================================
-[imgSIRT,errL2SIRT,qualitySIRT]=SIRT(noise_projections,geo,angles,30,...
+[imgSIRT,errL2SIRT,qualitySIRT]=SIRT(projections,geo,angles,30,...
                             'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas);
-[imgSART,errL2SART,qualitySART]=SART(noise_projections,geo,angles,30,...
+[imgSART,errL2SART,qualitySART]=SART(projections,geo,angles,30,...
                             'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas);
 % OS-SART
 % ========================================================================
@@ -138,7 +138,7 @@
 %                                     biggest angular distance with the
 %                                     ones used.  (default)
 order='angularDistance';
-[imgOSSART,errL2OSSART,qualityOSSART]=OS_SART(noise_projections,geo,angles,30,...
+[imgOSSART,errL2OSSART,qualityOSSART]=OS_SART(projections,geo,angles,30,...
                             'lambda',lambda,'lambda_red',lambdared,'verbose',verbose,'QualMeas',qualmeas,...
                              'BlockSize',blcks,'OrderStrategy',order);
 %% Lets have a brief show of the results
@@ -164,7 +164,7 @@
 plotImg([imgSIRT;  imgOSSART; imgSART;],'Dim','Z','Savegif','sarts.gif');
 
 % plot error
-plotImg(abs([thorax-imgSIRT thorax-imgSART thorax-imgOSSART]),'Dim','Z');
+plotImg(abs([head-imgSIRT; head-imgOSSART; head-imgSART; ]),'Dim','Z');
 
 
 

diff --git a/Demos/d08_Algorithms03.m b/Demos/d08_Algorithms03.m
@@ -52,9 +52,9 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
-noise_projections=addCTnoise(projections,'Poisson',1e4,'Gaussian',[0 50]);
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
+noise_projections=addCTnoise(projections);
 
 %% Usage CGLS
 %
@@ -95,4 +95,4 @@
 % plot images
 plotImg([imgCGLS imgSIRT],'Dim','Z','Step',2)
 %plot errors
-plotImg(abs([thorax-imgCGLS thorax-imgSIRT]),'Dim','Z','Slice',64)
+plotImg(abs([head-imgCGLS head-imgSIRT]),'Dim','Z','Slice',64)
diff --git a/Demos/d09_Algorithms04.m b/Demos/d09_Algorithms04.m
@@ -54,8 +54,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi-2*pi/30,30);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 
 %% Lets create a OS-SART test for comparison
@@ -205,10 +205,10 @@
 %    
 %     OSC-TV             B-ASD-POCS-beta   SART-TV
 
-plotImg([ imgOSCTV imgBASDPOCSbeta imgSARTTV; thorax imgOSSART  imgASDPOCS ] ,'Dim','Z','Step',2)
+plotImg([ imgOSCTV imgBASDPOCSbeta imgSARTTV; head imgOSSART  imgASDPOCS ] ,'Dim','Z','Step',2)
  % error
 
-plotImg(abs([ thorax-imgOSCTV thorax-imgBASDPOCSbeta thorax-imgSARTTV;thorax-thorax thorax-imgOSSART  thorax-imgASDPOCS ]) ,'Dim','Z','Slice',64)
+plotImg(abs([ head-imgOSCTV head-imgBASDPOCSbeta head-imgSARTTV;head-head head-imgOSSART  head-imgASDPOCS ]) ,'Dim','Z','Slice',64)
 
 
 

diff --git a/Demos/d10_QualityMeasurements.m b/Demos/d10_QualityMeasurements.m
@@ -47,8 +47,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 
 %% "Measure_Quality.m" function

diff --git a/Demos/d11_PostProcessing.m b/Demos/d11_PostProcessing.m
@@ -45,8 +45,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 %% Lets just use FDK
 
@@ -75,6 +75,6 @@
 % however, teh denoising has no knoledge of the original data (projections)
 % this it doesnt reduce the error. The error increases, specially in small
 % areas
-plotImg(abs([thorax-imgFDK thorax-imgdenoised thorax-imcroped]),'Dim','Z')
+plotImg(abs([head-imgFDK head-imgdenoised head-imcroped]),'Dim','Z')
 
 
diff --git a/Demos/d13_HelicalGeometry.m b/Demos/d13_HelicalGeometry.m
@@ -47,15 +47,15 @@
 angles=linspace(00.0,2*pi,100);
 angles=[angles angles angles];
 % Load thorax phatom data
-thorax=thoraxPhantom(geo.nVoxel); % yes, not the best example data, but It will do.
+head=headPhantom(geo.nVoxel); % yes, not the best example data, but It will do.
 
 % Thsi makes it helical
 geo.offOrigin(3,:)=linspace(-1024/2+128,1024/2-128,length(angles)).';  % about 256^3 images fit int he detector with this size. 
 geo.offOrigin(1,:)=0;
 geo.offOrigin(2,:)=0;
 
 % project data
-data=Ax(thorax,geo,angles);
+data=Ax(head,geo,angles);
 
 % Uncomment if you want to see the data
 % plotProj(data,angles);
@@ -67,6 +67,6 @@
 %% Plot results
 
 % CGLS and SIRT
-plotImg([thorax, OSSARTimg ,CGLSimg],'Dim',3,'Step',3);
+plotImg([head, OSSARTimg ,CGLSimg],'Dim',3,'Step',3);
 
 
diff --git a/Demos/d14_Offsets.m b/Demos/d14_Offsets.m
@@ -51,8 +51,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 
 
 %% lets see it
@@ -66,7 +66,7 @@
 %% Second test: lets test variying offsets:
 
 geo.offDetector=[10*sin(angles); 20*cos(angles)];                     % Offset of Detector            (mm)
-projections2=Ax(thorax,geo,angles,'interpolated');
+projections2=Ax(head,geo,angles,'interpolated');
 %% lets see it
 plotProj(projections2,angles);
 %% reconstruction
@@ -83,7 +83,7 @@
 geo.offDetector=[10*sin(angles); 10*cos(angles)];                     % Offset of Detector            (mm)
 geo.offOrigin =[40*sin(angles);linspace(-30,30,100);40*cos(angles)];                     % Offset of image from origin   (mm)              
 
-projections3=Ax(thorax,geo,angles,'interpolated');
+projections3=Ax(head,geo,angles,'interpolated');
 %% lets see it
 plotProj(projections3,angles);
 %% reconstruction

diff --git a/Demos/d15_3Dparallel.m b/Demos/d15_3Dparallel.m
@@ -53,8 +53,8 @@
 %% Load data and generate projections 
 % see previous demo for explanation
 angles=linspace(0,2*pi,100);
-thorax=thoraxPhantom(geo.nVoxel);
-projections=Ax(thorax,geo,angles,'interpolated');
+head=headPhantom(geo.nVoxel);
+projections=Ax(head,geo,angles,'interpolated');
 noise_projections=addCTnoise(projections);
 
 %% Reconsturction 

diff --git a/Demos/d16_2Dtomography.m b/Demos/d16_2Dtomography.m
@@ -66,9 +66,9 @@
 %% Define angles of projection and load phatom image
 
 angles=linspace(0,2*pi,100);
-thorax=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
-thorax=cat(3,thorax,thorax);
-projections=Ax(thorax,geo,angles,'interpolated');
+phatom=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
+phatom=cat(3,phatom,phatom);
+projections=Ax(phatom,geo,angles,'interpolated');
 %% recosntruct
 
 imgOSSART=OS_SART(projections,geo,angles,40);
@@ -124,9 +124,9 @@
 %% Define angles of projection and load phatom image
 
 angles=linspace(0,2*pi,100);
-thorax=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
-thorax=cat(3,thorax,thorax);
-projections=Ax(thorax,geo,angles,'interpolated');
+phatom=single(phantom('Modified Shepp-Logan',geo.nVoxel(1)));
+phatom=cat(3,phatom,phatom);
+projections=Ax(phatom,geo,angles,'interpolated');
 %% recosntruct
 
 imgOSSART=OS_SART(projections,geo,angles,40);

diff --git a/Test_data/MRheadbrain/head.mat b/Test_data/MRheadbrain/head.mat
diff --git a/Test_data/Thorax-phantom/thoraxPhantom.m → Test_data/MRheadbrain/headPhantom.m b/Test_data/Thorax-phantom/thoraxPhantom.m → Test_data/MRheadbrain/headPhantom.m
@@ -1,9 +1,9 @@
-function img=thoraxPhantom(varargin)
-%THORAXPHANTOM returns the thorax phantom 
+function img=headPhantom(varargin)
+%HEADXPHANTOM returns the head phantom 
 %
-%  IMG=THORAXPHANTOM() returns 128^3 image IMG
+%  IMG=HEADXPHANTOM() returns 128^3 image IMG
 %
-%  IMG=THORAXPHANTOM(SZ) returns a SZ^3 image IMG if SZ is scalar, or a [SZ(1)
+%  IMG=HEADXPHANTOM(SZ) returns a SZ^3 image IMG if SZ is scalar, or a [SZ(1)
 %   SZ(2) SZ(3] image, of SZ is a vector. 
 %--------------------------------------------------------------------------
 %--------------------------------------------------------------------------
@@ -43,7 +43,7 @@
     end
 end
 % load data
-data=load('img128.mat');
+data=load('head.mat');
 img=data.img;
 
 % interpolate data to get desired size

diff --git a/Test_data/Thorax-phantom/img128.mat b/Test_data/Thorax-phantom/img128.mat
diff --git a/Test_data/Thorax-phantom/license.txt b/Test_data/Thorax-phantom/license.txt
diff --git a/Utilities/TIGRE.bib b/Utilities/TIGRE.bib
@@ -361,6 +361,19 @@ @article{ollinger1994maximum
   publisher={IEEE}
 }
 
-
+headPhantom
+
+A. Aichert, M. Manhart, B. Navalpakkam, R. Grimm, J. Hutter, A. Maier, J. Hornegger, A. Doerfler
+A Realistic Digital Phantom for Perfusion C-Arm CT Based on MRI Data
+2013 IEEE Nuclear Science Symposium and Medical Imaging Conference Record (NSS/MIC), Seoul, Korea, 2013.
+
+@inproceedings{aichert2013realistic,
+  title={A realistic digital phantom for perfusion C-arm CT based on MRI data},
+  author={Aichert, Andr{\'e} and Manhart, Michael T and Navalpakkam, Bharath K and Grimm, Robert and Hutter, Jana and Maier, Andreas and Hornegger, Joachim and Doerfler, Arnd},
+  booktitle={2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC)},
+  pages={1--2},
+  year={2013},
+  organization={IEEE}
+}
 
 End of Bib (leave this here)
diff --git a/Utilities/addCTnoise.m b/Utilities/addCTnoise.m
@@ -66,7 +66,7 @@
     switch opt
         case 'Poisson'
             if default
-                I0=1e5;
+                I0=65535;
             else
                if ~isscalar(val);error('CBCT:addnoise:WorngInput','Input to Poisson should be scalar');end
                I0=val;
@@ -75,7 +75,7 @@
         case 'Gaussian'
             if default
                 m=0;
-                sigma=10;
+                sigma=0.5;
             else
                 if (size(val,2)~=2);error('CBCT:addnoise:WorngInput','Input to Gaussian should be 1x2');end
                 m=val(1);
@@ -88,7 +88,7 @@
 
 %% Add the noise
 %//Lambert-Beer
-Im=I0*exp(-proj);
+Im=I0*exp(-proj/max(proj(:)));
 
 % Photon noise + electronic noise
 if areTheseToolboxesInstalled({'MATLAB','Statistics Toolbox'}) || areTheseToolboxesInstalled({'MATLAB','Statistics and Machine Learning Toolbox'})
@@ -103,6 +103,6 @@
          'With I0 ~ 10000'])
      Im=randn(size(Im)).*sigma + m; % this one is slower
 end
-Im(Im<0)=1e-6;
-proj=single(log(I0./Im));
+Im(Im<=0)=1e-6;
+proj=single(log(I0./Im))*max(proj(:));
 end