From 0d8a3396615176665eb11cf9756fdbcf466e8ea2 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 16:49:00 +0800 Subject: [PATCH 01/19] Update checkGeo.m Added a parameter to enable different position of the source --- MATLAB/Utilities/checkGeo.m | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/MATLAB/Utilities/checkGeo.m b/MATLAB/Utilities/checkGeo.m index 442b76a4..52777273 100644 --- a/MATLAB/Utilities/checkGeo.m +++ b/MATLAB/Utilities/checkGeo.m @@ -6,14 +6,14 @@ '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 unknown=~ismember(fnames,allfields); % there must be not unknown variables % TODO: Do we really want to enforce this? Perhaps just a warning? -assert(~sum(unknown),'TIGRE:checkGeo:BadGeometry',['The following fields are not known by TIGRE:\n' strjoin(fnames(unknown)),'\nMake sure you have not misspelled any field or introduced unnecessary fields.']) +assert(~sum(unknown),'TIGRE:checkGeo:BadGeometry',['The following fields are not known by TIGRE:\n' strjoin(fnames(unknown)),'\nMake sure you have not misspelled any field or introduced unnecesary fields.']) @@ -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 From 8ce3d2b2f93d1e49d5d49addbc1dd65a94d15dec Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 16:50:12 +0800 Subject: [PATCH 02/19] Update computeV.m added the offSource parameter --- MATLAB/Utilities/computeV.m | 57 +++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/MATLAB/Utilities/computeV.m b/MATLAB/Utilities/computeV.m index 2150c703..82fcf649 100644 --- a/MATLAB/Utilities/computeV.m +++ b/MATLAB/Utilities/computeV.m @@ -25,6 +25,63 @@ auxgeo.DSO = geo.DSO(auxindex); auxgeo.offOrigin = geo.offOrigin(:,auxindex); auxgeo.offDetector = geo.offDetector(:,auxindex); + function V=computeV(geo,angles,alphablocks,orig_index,varargin) +if nargin>=6 + [gpuids] = parse_inputs(varargin{1:length(varargin)}); +elseif nargin==4 + gpuids = GpuIds(); +else + error('Wrong amount of inputs, 4 or 6 expected'); +end + +V=zeros(geo.nVoxel(1),geo.nVoxel(2) ,length(alphablocks),'single'); +geo=checkGeo(geo,angles); + +if ~isfield(geo,'mode')||~strcmp(geo.mode,'parallel') + for ii=1:length(alphablocks) + auxang=alphablocks{ii}; + auxindex=orig_index{ii}; + auxgeo = geo; + % expand the detector to avoiding zeros in backprojection +% maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); +% auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); +% auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; + % subset of projection angles + auxgeo.DSD = geo.DSD(auxindex); + auxgeo.DSO = geo.DSO(auxindex); + auxgeo.offOrigin = geo.offOrigin(:,auxindex); + auxgeo.offDetector = geo.offDetector(:,auxindex); + auxgeo.offSource = geo.offSource(:,auxindex); + auxgeo.rotDetector = geo.rotDetector(:,auxindex); + auxgeo.COR = geo.COR(auxindex); + %auxgeo=geo; + V(:,:,ii) = mean(Atb(ones(geo.nDetector(2),geo.nDetector(1),size(auxang,2),'single'),auxgeo,auxang,'gpuids',gpuids),3)+0.000001; + end + V(V==0.0)=Inf; +else + for ii=1:length(alphablocks) + V(:,:,ii)=ones(geo.nVoxel(1:2).','single')*length(alphablocks{ii}); + end +end +end + +function [gpuids]=parse_inputs(varargin) + %fprintf('parse_inputs0(varargin (%d))\n', length(varargin)); + if isempty(varargin) + gpuids = GpuIds(); + else + % create input parser + p=inputParser; + % add optional parameters + addParameter(p,'gpuids', GpuIds()); + %execute + parse(p,varargin{:}); + %extract + gpuids=p.Results.gpuids; + end +end + + auxgeo.rotDetector = geo.rotDetector(:,auxindex); auxgeo.COR = geo.COR(auxindex); %auxgeo=geo; From 28b4cd070b4eddaaf96f55354e9ceb44d19518db Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 16:54:15 +0800 Subject: [PATCH 03/19] Update Atb_mex.cpp added offSource to enable different positions of X-ray source --- MATLAB/Utilities/cuda_interface/Atb_mex.cpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/MATLAB/Utilities/cuda_interface/Atb_mex.cpp b/MATLAB/Utilities/cuda_interface/Atb_mex.cpp index da78bfce..9bcf6af9 100644 --- a/MATLAB/Utilities/cuda_interface/Atb_mex.cpp +++ b/MATLAB/Utilities/cuda_interface/Atb_mex.cpp @@ -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,7 +190,7 @@ 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; @@ -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 Date: Wed, 4 Jun 2025 16:55:54 +0800 Subject: [PATCH 04/19] Update Ax_mex.cpp added offSource --- MATLAB/Utilities/cuda_interface/Ax_mex.cpp | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/MATLAB/Utilities/cuda_interface/Ax_mex.cpp b/MATLAB/Utilities/cuda_interface/Ax_mex.cpp index 3c6f3670..1f274501 100644 --- a/MATLAB/Utilities/cuda_interface/Ax_mex.cpp +++ b/MATLAB/Utilities/cuda_interface/Ax_mex.cpp @@ -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; @@ -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 Date: Wed, 4 Jun 2025 16:58:47 +0800 Subject: [PATCH 05/19] Update ray_interpolated_projection.cu Enable different source positions by using offSource --- Common/CUDA/ray_interpolated_projection.cu | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) 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; From 6e43bd69c32724d563f97eef728b62afcb0105e4 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 17:03:21 +0800 Subject: [PATCH 06/19] Update Siddon_projection.cu Enable different source positions by using offSource --- Common/CUDA/Siddon_projection.cu | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) 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; From e30cc34fb19a76fafe219da0a4ec9c9bd5a02d6d Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 17:05:35 +0800 Subject: [PATCH 07/19] Update d01_CreateGeometry.m added offSource --- MATLAB/Demos/d01_CreateGeometry.m | 1 + 1 file changed, 1 insertion(+) 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 From 898ac62febe9f638e2498c0de85d7e2bdbd99eb8 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 17:06:29 +0800 Subject: [PATCH 08/19] Update defaultGeometry.m added offSource --- MATLAB/Utilities/defaultGeometry.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/MATLAB/Utilities/defaultGeometry.m b/MATLAB/Utilities/defaultGeometry.m index fe566758..5dcf1f19 100644 --- a/MATLAB/Utilities/defaultGeometry.m +++ b/MATLAB/Utilities/defaultGeometry.m @@ -27,6 +27,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]; % Offset of Source (mm) geo.offOrigin =[0;0;0]; % Offset of image from origin (mm) geo.offDetector=[0; 0]; % Offset of Detector (mm) @@ -104,4 +105,4 @@ warning('In Parallel mode nVoxel(2:3) is generally equal to nDetector. Consider setting them equal for better reconstruction quality.'); end -end \ No newline at end of file +end From 290d10316a4e6f8b3fc352ca72cb97ee709500ad Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Wed, 4 Jun 2025 17:15:43 +0800 Subject: [PATCH 09/19] Update computeV.m added offSource --- MATLAB/Utilities/computeV.m | 57 ------------------------------------- 1 file changed, 57 deletions(-) diff --git a/MATLAB/Utilities/computeV.m b/MATLAB/Utilities/computeV.m index 82fcf649..2150c703 100644 --- a/MATLAB/Utilities/computeV.m +++ b/MATLAB/Utilities/computeV.m @@ -25,63 +25,6 @@ auxgeo.DSO = geo.DSO(auxindex); auxgeo.offOrigin = geo.offOrigin(:,auxindex); auxgeo.offDetector = geo.offDetector(:,auxindex); - function V=computeV(geo,angles,alphablocks,orig_index,varargin) -if nargin>=6 - [gpuids] = parse_inputs(varargin{1:length(varargin)}); -elseif nargin==4 - gpuids = GpuIds(); -else - error('Wrong amount of inputs, 4 or 6 expected'); -end - -V=zeros(geo.nVoxel(1),geo.nVoxel(2) ,length(alphablocks),'single'); -geo=checkGeo(geo,angles); - -if ~isfield(geo,'mode')||~strcmp(geo.mode,'parallel') - for ii=1:length(alphablocks) - auxang=alphablocks{ii}; - auxindex=orig_index{ii}; - auxgeo = geo; - % expand the detector to avoiding zeros in backprojection -% maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); -% auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); -% auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; - % subset of projection angles - auxgeo.DSD = geo.DSD(auxindex); - auxgeo.DSO = geo.DSO(auxindex); - auxgeo.offOrigin = geo.offOrigin(:,auxindex); - auxgeo.offDetector = geo.offDetector(:,auxindex); - auxgeo.offSource = geo.offSource(:,auxindex); - auxgeo.rotDetector = geo.rotDetector(:,auxindex); - auxgeo.COR = geo.COR(auxindex); - %auxgeo=geo; - V(:,:,ii) = mean(Atb(ones(geo.nDetector(2),geo.nDetector(1),size(auxang,2),'single'),auxgeo,auxang,'gpuids',gpuids),3)+0.000001; - end - V(V==0.0)=Inf; -else - for ii=1:length(alphablocks) - V(:,:,ii)=ones(geo.nVoxel(1:2).','single')*length(alphablocks{ii}); - end -end -end - -function [gpuids]=parse_inputs(varargin) - %fprintf('parse_inputs0(varargin (%d))\n', length(varargin)); - if isempty(varargin) - gpuids = GpuIds(); - else - % create input parser - p=inputParser; - % add optional parameters - addParameter(p,'gpuids', GpuIds()); - %execute - parse(p,varargin{:}); - %extract - gpuids=p.Results.gpuids; - end -end - - auxgeo.rotDetector = geo.rotDetector(:,auxindex); auxgeo.COR = geo.COR(auxindex); %auxgeo=geo; From e5b20a7321608f20010bae1732645cd7ecbc46fb Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 22:49:21 +0800 Subject: [PATCH 10/19] Update types_TIGRE.hpp add offSource --- Common/CUDA/types_TIGRE.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 From 4fafecaefe10d1fdbad9b483b5e0857c01d3cf32 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 22:50:49 +0800 Subject: [PATCH 11/19] Update Atb_mex.cpp fix bug caused by offSource --- MATLAB/Utilities/cuda_interface/Atb_mex.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/MATLAB/Utilities/cuda_interface/Atb_mex.cpp b/MATLAB/Utilities/cuda_interface/Atb_mex.cpp index 9bcf6af9..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"; @@ -196,7 +196,7 @@ void mexFunction(int nlhs , mxArray *plhs[], 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 From c4223b9519d903daea8c3825aed0fd302aa56ba9 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 22:51:29 +0800 Subject: [PATCH 12/19] Update Ax_mex.cpp fix bug caused by offSource --- MATLAB/Utilities/cuda_interface/Ax_mex.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/MATLAB/Utilities/cuda_interface/Ax_mex.cpp b/MATLAB/Utilities/cuda_interface/Ax_mex.cpp index 1f274501..e1b4c497 100644 --- a/MATLAB/Utilities/cuda_interface/Ax_mex.cpp +++ b/MATLAB/Utilities/cuda_interface/Ax_mex.cpp @@ -141,7 +141,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"; @@ -171,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 From a19066fb688deaae3ee4819fa4788bd96e5d56ba Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:03:24 +0800 Subject: [PATCH 13/19] Add files via upload add offSource --- Common/CUDA/voxel_backprojection.cu | 4 ++-- Common/CUDA/voxel_backprojection2.cu | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) 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); From 8df67e7b8816211f047117cfbdfc76352d49d737 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:05:54 +0800 Subject: [PATCH 14/19] add offSource --- MATLAB/Demos/d025_OffsetSource.m | 46 ++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 MATLAB/Demos/d025_OffsetSource.m diff --git a/MATLAB/Demos/d025_OffsetSource.m b/MATLAB/Demos/d025_OffsetSource.m new file mode 100644 index 00000000..6504586a --- /dev/null +++ b/MATLAB/Demos/d025_OffsetSource.m @@ -0,0 +1,46 @@ +%% 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; +run 'D:\dyf20\CT\TIGRE-offSource\MATLAB\InitTIGRE.m' +%% 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)],[]) \ No newline at end of file From 6b86348fc836801f41216da6a65a7b82a5ec59b4 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:08:32 +0800 Subject: [PATCH 15/19] add offSource --- MATLAB/Utilities/computeV.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/MATLAB/Utilities/computeV.m b/MATLAB/Utilities/computeV.m index 2150c703..5d149fc7 100644 --- a/MATLAB/Utilities/computeV.m +++ b/MATLAB/Utilities/computeV.m @@ -8,7 +8,6 @@ end V=zeros(geo.nVoxel(1),geo.nVoxel(2) ,length(alphablocks),'single'); -%geo.offDetector(1) = geo.offDetector(1) + (geo.DSD(1) / geo.DSO(1)) * geo.COR(1); geo=checkGeo(geo,angles); if ~isfield(geo,'mode')||~strcmp(geo.mode,'parallel') @@ -17,14 +16,15 @@ auxindex=orig_index{ii}; auxgeo = geo; % expand the detector to avoiding zeros in backprojection - maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); - auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); - auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; +% maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); +% auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); +% auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; % subset of projection angles auxgeo.DSD = geo.DSD(auxindex); auxgeo.DSO = geo.DSO(auxindex); auxgeo.offOrigin = geo.offOrigin(:,auxindex); auxgeo.offDetector = geo.offDetector(:,auxindex); + auxgeo.offSource = geo.offSource(:,auxindex); auxgeo.rotDetector = geo.rotDetector(:,auxindex); auxgeo.COR = geo.COR(auxindex); %auxgeo=geo; From 14f252baf0a936695ba4cf79d721f2fa30b0fab0 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:09:15 +0800 Subject: [PATCH 16/19] Add files via upload From 02f3ba3eb03c43a01175fabdf6e5115285ec697d Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:12:00 +0800 Subject: [PATCH 17/19] Update checkGeo.m --- MATLAB/Utilities/checkGeo.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MATLAB/Utilities/checkGeo.m b/MATLAB/Utilities/checkGeo.m index 52777273..b98e5989 100644 --- a/MATLAB/Utilities/checkGeo.m +++ b/MATLAB/Utilities/checkGeo.m @@ -13,7 +13,7 @@ unknown=~ismember(fnames,allfields); % there must be not unknown variables % TODO: Do we really want to enforce this? Perhaps just a warning? -assert(~sum(unknown),'TIGRE:checkGeo:BadGeometry',['The following fields are not known by TIGRE:\n' strjoin(fnames(unknown)),'\nMake sure you have not misspelled any field or introduced unnecesary fields.']) +assert(~sum(unknown),'TIGRE:checkGeo:BadGeometry',['The following fields are not known by TIGRE:\n' strjoin(fnames(unknown)),'\nMake sure you have not misspelled any field or introduced unnecessary fields.']) From 03157276e745b754e93b3f939dd097027d1563c2 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:15:46 +0800 Subject: [PATCH 18/19] Update computeV.m --- MATLAB/Utilities/computeV.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/MATLAB/Utilities/computeV.m b/MATLAB/Utilities/computeV.m index 5d149fc7..b002cbc1 100644 --- a/MATLAB/Utilities/computeV.m +++ b/MATLAB/Utilities/computeV.m @@ -8,6 +8,7 @@ end V=zeros(geo.nVoxel(1),geo.nVoxel(2) ,length(alphablocks),'single'); +%geo.offDetector(1) = geo.offDetector(1) + (geo.DSD(1) / geo.DSO(1)) * geo.COR(1); geo=checkGeo(geo,angles); if ~isfield(geo,'mode')||~strcmp(geo.mode,'parallel') @@ -16,15 +17,15 @@ auxindex=orig_index{ii}; auxgeo = geo; % expand the detector to avoiding zeros in backprojection -% maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); -% auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); -% auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; + maxsize=max(auxgeo.sVoxel+geo.offOrigin(:,auxindex),[],2); + auxgeo.sDetector=max(auxgeo.sDetector , [maxsize(1); maxsize(3)] *geo.DSD/geo.DSO); + auxgeo.dDetector = auxgeo.sDetector ./ auxgeo.nDetector; % subset of projection angles auxgeo.DSD = geo.DSD(auxindex); auxgeo.DSO = geo.DSO(auxindex); auxgeo.offOrigin = geo.offOrigin(:,auxindex); - auxgeo.offDetector = geo.offDetector(:,auxindex); auxgeo.offSource = geo.offSource(:,auxindex); + auxgeo.offDetector = geo.offDetector(:,auxindex); auxgeo.rotDetector = geo.rotDetector(:,auxindex); auxgeo.COR = geo.COR(auxindex); %auxgeo=geo; @@ -53,4 +54,3 @@ gpuids=p.Results.gpuids; end end - From e4823f85332f51b8feed44e3ed2fe75621b0bfa8 Mon Sep 17 00:00:00 2001 From: defThu <167285267+defThu@users.noreply.github.com> Date: Fri, 6 Jun 2025 23:49:39 +0800 Subject: [PATCH 19/19] Add a demo for simulations with source/focal spot offset --- MATLAB/Demos/d025_OffsetSource.m | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/MATLAB/Demos/d025_OffsetSource.m b/MATLAB/Demos/d025_OffsetSource.m index 6504586a..190d4fff 100644 --- a/MATLAB/Demos/d025_OffsetSource.m +++ b/MATLAB/Demos/d025_OffsetSource.m @@ -22,7 +22,6 @@ %% Initialize clear; close all; -run 'D:\dyf20\CT\TIGRE-offSource\MATLAB\InitTIGRE.m' %% Geometry geo=defaultGeometry(); head=headPhantom(geo.nVoxel); @@ -43,4 +42,4 @@ % 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)],[]) \ No newline at end of file +imtool([img_ref(:,:,128),img(:,:,128)],[])