diff --git a/Common/CUDA/Siddon_projection.cu b/Common/CUDA/Siddon_projection.cu index 2a025f8c..9939cbd6 100644 --- a/Common/CUDA/Siddon_projection.cu +++ b/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; diff --git a/Common/CUDA/ray_interpolated_projection.cu b/Common/CUDA/ray_interpolated_projection.cu index e71c5b59..bc974cb7 100644 --- a/Common/CUDA/ray_interpolated_projection.cu +++ b/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; diff --git a/Common/CUDA/types_TIGRE.hpp b/Common/CUDA/types_TIGRE.hpp index 0a3abc4d..a9cceb93 100644 --- a/Common/CUDA/types_TIGRE.hpp +++ b/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 \ No newline at end of file +#endif diff --git a/Common/CUDA/voxel_backprojection.cu b/Common/CUDA/voxel_backprojection.cu index bec4d909..8f3b8798 100644 --- a/Common/CUDA/voxel_backprojection.cu +++ b/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; diff --git a/Common/CUDA/voxel_backprojection2.cu b/Common/CUDA/voxel_backprojection2.cu index c9dcc957..11dc8ae9 100644 --- a/Common/CUDA/voxel_backprojection2.cu +++ b/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); diff --git a/MATLAB/Demos/d01_CreateGeometry.m b/MATLAB/Demos/d01_CreateGeometry.m index 9753e210..b49d8e4d 100644 --- a/MATLAB/Demos/d01_CreateGeometry.m +++ b/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 diff --git a/MATLAB/Demos/d025_OffsetSource.m b/MATLAB/Demos/d025_OffsetSource.m new file mode 100644 index 00000000..190d4fff --- /dev/null +++ b/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: tigre.toolbox@gmail.com +% 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)],[]) diff --git a/MATLAB/Utilities/checkGeo.m b/MATLAB/Utilities/checkGeo.m index 442b76a4..b98e5989 100644 --- a/MATLAB/Utilities/checkGeo.m +++ b/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 diff --git a/MATLAB/Utilities/computeV.m b/MATLAB/Utilities/computeV.m index 2150c703..b002cbc1 100644 --- a/MATLAB/Utilities/computeV.m +++ b/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 - diff --git a/MATLAB/Utilities/cuda_interface/Atb_mex.cpp b/MATLAB/Utilities/cuda_interface/Atb_mex.cpp index da78bfce..b61185ff 100644 --- a/MATLAB/Utilities/cuda_interface/Atb_mex.cpp +++ b/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 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"; @@ -156,12 +156,13 @@ void mexFunction(int nlhs , mxArray *plhs[], fieldnames[11]= "mode"; fieldnames[12]= "COR"; fieldnames[13]= "rotDetector"; - + fieldnames[14]= "offSource"; + // Now we know that all the input struct is good! Parse it from mxArrays to // C structures that MEX can understand. double * nVoxel, *nDetec; //we need to cast these to int double * sVoxel, *dVoxel,*sDetec,*dDetec, *DSO, *DSD; - double *offOrig,*offDetec,*rotDetector; + double *offOrig,*offDetec,*rotDetector,*offSource; double * acc, *COR; const char* mode; int c; @@ -170,7 +171,7 @@ void mexFunction(int nlhs , mxArray *plhs[], geo.unitX=1;geo.unitY=1;geo.unitZ=1; bool coneBeam=true; // mexPrintf("%d \n",nfields); - for(int ifield=0; ifield<14; ifield++) { + for(int ifield=0; ifield<15; ifield++) { tmp=mxGetField(geometryMex,0,fieldnames[ifield]); if(tmp==NULL){ //tofix @@ -290,6 +291,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