Merge pull request #666 from defThu/offSource

added "offSource" to adjust the position of X-ray source

Author: AnderBiguri

Common/CUDA/Siddon_projection.cu

@@ -668,9 +668,11 @@ void computeDeltas_Siddon(Geometry geo,int i, Point3D* uvorigin, Point3D* deltaU
 
     Point3D S;
     S.x=geo.DSO[i];
-    S.y=0;
-    S.z=0;
-    
+    //S.y=0;
+    //S.z=0;
+    S.y=geo.offSourceY[i];
+    S.z=geo.offSourceZ[i];
+
     //End point
     Point3D P,Pu0,Pv0;
 

Common/CUDA/ray_interpolated_projection.cu

@@ -632,9 +632,11 @@ void splitImageInterp(unsigned int splits,Geometry geo,Geometry* geoArray, unsig
 void computeDeltas(Geometry geo,unsigned int i, Point3D* uvorigin, Point3D* deltaU, Point3D* deltaV, Point3D* source){
     Point3D S;
     S.x=geo.DSO[i];
-    S.y=0;
-    S.z=0;
-    
+    //S.y=0;
+    //S.z=0;
+    S.y=geo.offSourceY[i];
+    S.z=geo.offSourceZ[i];
+
     //End point
     Point3D P,Pu0,Pv0;
 

Common/CUDA/types_TIGRE.hpp

@@ -63,6 +63,7 @@ struct  Geometry {
     float sDetecU, sDetecV;
     float dDetecU, dDetecV;
     float *offDetecU, *offDetecV;
+    float *offSourceY,*offSourceZ;
     float* DSD;
     float* dRoll;
     float* dPitch;
@@ -106,4 +107,4 @@ struct Point3Ddouble{
     }
 };
 
-#endif
+#endif

Common/CUDA/voxel_backprojection.cu

@@ -871,8 +871,8 @@ void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX,
     //Done for P, now source
     Point3Ddouble source;
     source.x=geo.DSD[i]; //already offseted for rotation
-    source.y=-geo.offDetecU[i];
-    source.z=-geo.offDetecV[i];
+    source.y=geo.offSourceY[i]-geo.offDetecU[i];
+    source.z=geo.offSourceZ[i]-geo.offDetecV[i];
     rollPitchYawT(geo,i,&source);
 
     source.x=source.x-(geo.DSD[i]-geo.DSO[i]);//   source.y=source.y-auxOff.y;    source.z=source.z-auxOff.z;

Common/CUDA/voxel_backprojection2.cu

@@ -793,8 +793,8 @@ void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX,
     //Done for P, now source
     Point3Ddouble source;
     source.x=geo.DSD[i]; //already offseted for rotation
-    source.y=-geo.offDetecU[i];
-    source.z=-geo.offDetecV[i];
+    source.y=geo.offSourceY[i]-geo.offDetecU[i];
+    source.z=geo.offSourceZ[i]-geo.offDetecV[i];
     rollPitchYawT(geo,i,&source);
 
 

MATLAB/Demos/d01_CreateGeometry.m

@@ -77,6 +77,7 @@
 geo.sVoxel=[256;256;256];                   % total size of the image       (mm)
 geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)
 % Offsets
+geo.offSource = [0;0];
 geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)              
 geo.offDetector=[0; 0];                     % Offset of Detector            (mm)
                                             % These two can be also defined

MATLAB/Demos/d025_OffsetSource.m

@@ -0,0 +1,45 @@
+%% DEMO 03: Generate sample data and add realistic CT noise to it.
+%
+% This demo will show how to generate sample data for image reconstruction
+%
+%
+%--------------------------------------------------------------------------
+%--------------------------------------------------------------------------
+% This file is part of the TIGRE Toolbox
+% 
+% Copyright (c) 2015, University of Bath and 
+%                     CERN-European Organization for Nuclear Research
+%                     All rights reserved.
+%
+% License:            Open Source under BSD. 
+%                     See the full license at
+%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
+%
+% Contact:            [email protected]
+% Codes:              https://github.com/CERN/TIGRE/
+% Coded by:           Ander Biguri 
+%--------------------------------------------------------------------------
+%% Initialize
+clear;
+close all;
+%% Geometry
+geo=defaultGeometry();
+head=headPhantom(geo.nVoxel);
+angles = linspace(0,2*pi,500);
+P0=Ax(head,geo,angles,'interpolated');
+
+img_ref = FDK(P0,geo,angles);
+imtool(img_ref(:,:,128),[])
+
+offset = 60;
+geo1 = geo;
+geo1.offSource = [offset;0];
+geo1.offDetector = [0;0];
+geo1.offOrigin = [0;0;0];
+P1=Ax(head,geo1,angles,'interpolated');
+
+img = FDK(P1,geo1,angles);
+% show the difference caused by the offset of source
+imtool([P0(:,:,100)-P1(:,:,100)],[])
+% show the similar reconstruction images
+imtool([img_ref(:,:,128),img(:,:,128)],[])

MATLAB/Utilities/checkGeo.m

@@ -6,7 +6,7 @@
     'DSO','DSD'};
 
 geofields_optional={'offOrigin','offDetector','rotDetector','COR',...
-    'mode','accuracy'};
+    'mode','accuracy','offSource'};
 allfields=horzcat(geofields_mandatory,geofields_optional);
 fnames=fieldnames(geo);
 % Find if geo has fields we do not recongize
@@ -128,5 +128,14 @@
 end
 
 
+if isfield(geo,'offSource')
+    assert(isequal(size(geo.offSource),[2 1]) | isequal(size(geo.offSource),[2 size(angles,2)]),'TIGRE:checkGeo:BadGeometry','geo.offSource Should be 2x1 or 2xsize(angles,2)')
+    assert(isa(geo.offSource,'double'),'TIGRE:checkGeo:BadGeometry','Field geo.offSource is not double type.' )
+    if isequal(size(geo.offSource),[2 1])
+        geo.offSource=repmat(geo.offSource,[1, size(angles,2)]);
+    end
+else
+    geo.offSource=zeros(2,size(angles,2));
+end
 end
 

MATLAB/Utilities/computeV.m

@@ -24,6 +24,7 @@
         auxgeo.DSD = geo.DSD(auxindex);
         auxgeo.DSO = geo.DSO(auxindex);
         auxgeo.offOrigin = geo.offOrigin(:,auxindex);
+        auxgeo.offSource = geo.offSource(:,auxindex);
         auxgeo.offDetector = geo.offDetector(:,auxindex);
         auxgeo.rotDetector = geo.rotDetector(:,auxindex);
         auxgeo.COR = geo.COR(auxindex);
@@ -53,4 +54,3 @@
         gpuids=p.Results.gpuids;
     end
 end
-

MATLAB/Utilities/cuda_interface/Atb_mex.cpp

@@ -165,7 +165,7 @@ void mexFunction(int  nlhs , mxArray *plhs[],
     mxArray * geometryMex=(mxArray*)prhs[1];
 
     // IMPORTANT-> Make sure Matlab creates the struct in this order.
-    const char *fieldnames[14];
+    const char *fieldnames[15];
     fieldnames[0] = "nVoxel";
     fieldnames[1] = "sVoxel";
     fieldnames[2] = "dVoxel";
@@ -180,6 +180,7 @@ void mexFunction(int  nlhs , mxArray *plhs[],
     fieldnames[11]= "mode";
     fieldnames[12]= "COR";
     fieldnames[13]= "rotDetector";
+    fieldnames[14]= "offSource";
     // Make sure input is structure
 
     mxArray    *tmp;
@@ -189,13 +190,13 @@ void mexFunction(int  nlhs , mxArray *plhs[],
 
     double * nVoxel, *nDetec; //we need to cast these to int
     double * sVoxel, *dVoxel,*sDetec,*dDetec, *DSO, *DSD,*offOrig,*offDetec;
-    double *acc, *COR,*rotDetector;
+    double *acc, *COR,*rotDetector,*offSource;
     const char* mode;
     bool coneBeam=true;
     Geometry geo;
     int c;
     geo.unitX=1;geo.unitY=1;geo.unitZ=1;
-    for(int ifield=0; ifield<14; ifield++) {
+    for(int ifield=0; ifield<15; ifield++) {
         tmp=mxGetField(geometryMex,0,fieldnames[ifield]);
         if(tmp==NULL){
             //tofix
@@ -314,6 +315,17 @@ void mexFunction(int  nlhs , mxArray *plhs[],
 
                 }
                 break;
+            case 14:
+                geo.offSourceY=(float*)malloc(nangles * sizeof(float));
+                geo.offSourceZ=(float*)malloc(nangles * sizeof(float));
+                
+                offSource=(double *)mxGetData(tmp);
+                for (int i=0;i<nangles;i++){
+                    c=i;
+                    geo.offSourceY[i]=(float)offSource[0+2*c];
+                    geo.offSourceZ[i]=(float)offSource[1+2*c];
+                }
+                break;
             default:
                 mexErrMsgIdAndTxt( "CBCT:MEX:Atb:unknown","This should not happen. Weird");
                 break;