diff --git a/Common/data/NIST_material_attenuation.json b/Common/data/NIST_material_attenuation.json new file mode 100644 index 00000000..6f6a11d2 --- /dev/null +++ b/Common/data/NIST_material_attenuation.json @@ -0,0 +1 @@ +{"Al":{"energy":[20,30,40,50,60,80,10,110,120,130,140,150],"mu":[0.9291,0.3046,0.1535,0.0994,0.075,0.0545,0.046,0.0431,0.0414,0.0398,0.0384,0.0372]},"PMMA":{"energy":[20,30,40,50,60,80,10,110,120,130,140,150],"mu":[0.0677,0.0359,0.0278,0.0246,0.0228,0.0207,0.0194,0.0188,0.0184,0.0179,0.0175,0.0172]}} \ No newline at end of file diff --git a/Common/data/dual_energy_calibrations.json b/Common/data/dual_energy_calibrations.json new file mode 100644 index 00000000..396f78fa --- /dev/null +++ b/Common/data/dual_energy_calibrations.json @@ -0,0 +1 @@ +{"VarianTrueBeam":{"V_90_140":{"Al":[80.397,-85.329,96.637,-200.79,104.18,4.2368,-21.295,30.648,-13.608],"PMMA":[-279.8,341.65,-211.44,444.98,-233.62,-14.54,65.484,-89.464,38.571]},"V_80_140":{"Al":[80,-87.778,11.875,-27.992,16.506,-1.3622,4.1877,-4.1797,1.3352],"PMMA":[-244.29,312.13,-12.957,36.246,-24.015,2.3154,-7.553,7.973,-2.7003]}}} \ No newline at end of file diff --git a/MATLAB/Demos/d23_DualEnergy.m b/MATLAB/Demos/d23_DualEnergy.m new file mode 100644 index 00000000..1922d0d0 --- /dev/null +++ b/MATLAB/Demos/d23_DualEnergy.m @@ -0,0 +1,46 @@ +%% Demo 23: Dual Energy CBCT in TIGRE. +% +% This demo will illustrate the options for performing dual-energy analysis +% using the TIGRE toolkit. +% +% Coded by: Andrew Keeler +% Contact: akeeler@luc.edu +% Date: February 28, 2024 + +%% Data loading (See d21_ScannerDataLoader.m) +% We'll use VarianDataLoader as an example. Load two sets of scan data +% at different kV settings into memory. Geo should be the same for both +% scans, and so only needs to be loaded once. Beam hardening correction +% is not needed. + +datafolder_low='~/your_data_path/varian/low_kV/2020-01-01_123456/'; +datafolder_high='~/your_data_path/varian/high_kV/2020-01-01_123456/'; +[proj_h, geo, angles_h] = VarianDataLoader(datafolder_high, 'bh', false); +[proj_l, ~, angles_l] = VarianDataLoader(datafolder, 'bh', false); + +%% Material Decomposition +% DE-CBCT decomposes the two scans into equivalent thicknesses of two basis +% materials. Generally, either Al and PMMA or Iodine and Water are used. +% The calibration in TIGRE has uses Al and PMM, but other materials are +% supported if the calibration files are updated with the data (contact +% developers or make a PR for expanding this) +% +% The decomposition procedure defaults to taking scans at 80 and 140 kV. +% Different energies can be specified using the 'kVl' and 'kVh' tags. + +[Al_proj, PMMA_proj, angles] = DeDecompose(proj_l, angles_l, proj_h, angles_h, ... % The following are defaults + 'scanner','VarianTrueBeam','material1','Al','material2', 'PMMA','kVl', 80, 'kVh', 140); + +%% Image construction +% Equivalent thickness projections can now be used to synthesize various +% DE-derived image types. The most common of these are virtual +% monoenergetic (VM) images. + +energy = 60; % keV +VM_proj = MakeVMproj(Al_proj, PMMA_proj, energy,'material1','Al','material2','PMMA'); + +VMimg = FDK(VM_proj, geo, angles); +VMimg = OS_SART(VM_proj, geo, angles, 100); + +%% Plot image +plotImg(VMimg) \ No newline at end of file diff --git a/MATLAB/InitTIGRE.m b/MATLAB/InitTIGRE.m index b101f7c5..c49304b4 100644 --- a/MATLAB/InitTIGRE.m +++ b/MATLAB/InitTIGRE.m @@ -30,6 +30,7 @@ addpath('./Utilities/IO/Diondo'); addpath('./Utilities/IO/DXChange'); addpath('./Utilities/GPU'); +addpath('./Utilities/DE'); addpath(genpath('./Test_data')); diff --git a/MATLAB/Utilities/DE/DeDecompose.m b/MATLAB/Utilities/DE/DeDecompose.m new file mode 100644 index 00000000..0b2c42df --- /dev/null +++ b/MATLAB/Utilities/DE/DeDecompose.m @@ -0,0 +1,216 @@ +function [material1_proj,material2_proj, angles_out] = DeDecompose(proj_l, angles_l , proj_h, angles_h, varargin) +%DEDECOMPSE Maps CBCT projections taken at two kVp settings to equivalent +%thicknesses of material1 and material2 using precalibrated transfer functions. +% Coefficients C and D are the coefficients for the third-degree transfer +% functions (may need to be modified for different systems) +% +% Other basis materials and kV settings can be employed by adding +% additional calibration data to "de_cal.mat" +% +% Proj_l and angles_l are the projection and angle data for the low-kV +% scan +% +% Proj_h and angles_h are the projection and angle data for the high-kV +% scan +% +% Optional Arguments +% scanner - default VarianTrueVeam +% +% kVl, kVh - the kV settings for the low- and high-kV scans, +% respectively. defaults: 80, 140 +% +% interp_model - model for sinogram interpolation +% +% 'nearest' - nearest neighbor (default). +% +% 'skip' - same as 'nearest, but skips projections without close +% matches. Requires setting 'tolerance' value. +% +% 'linear' - linear interpolation (unreliable for large +% projection separations) +% +% Coded by: Andrew Keeler +% Contact: akeeler@luc.edu +% Date: February 28, 2024 + +% Load calibrated coefficients for transfer functions + + +[scanner, calibration_pair, model, tolerance, material1, material2] = parse_inputs(varargin); + +% read NIST database (values in TIGRE) +fid = fopen('./../../../Common/data/dual_energy_calibrations.json'); +raw = fread(fid,inf); +str = char(raw'); +fclose(fid); +calibration_data = jsondecode(str); + +if ~isfield(calibration_data, scanner) + error(['Scanner: ', scanner, ' has no calibration data']) +end +if ~isfield(calibration_data.(scanner),calibration_pair) + error(['Calibration pair of energies: ', calibration_pair, ' not found in the scanner data']) +end +materials = fieldnames(calibration_data.(scanner).(calibration_pair)); + +if ~isfield(calibration_data.(scanner).(calibration_pair),material1) + disp([material1, ' not a (yet) supported material in the scanner at the given energies']) + disp("Supported materials:") + for i=1:len(materials) + disp(materials{i}) + end + error("") +end +if ~isfield(calibration_data.(scanner).(calibration_pair),material2) + disp([material2, ' not a (yet) supported material in the scanner at the given energies']) + disp("Supported materials:") + for i=1:len(materials) + disp(materials{i}) + end + error("") +end + +c = calibration_data.(scanner).(calibration_pair).(material1); +d = calibration_data.(scanner).(calibration_pair).(material2); + +%% check that the projections are the same size +for n = 1:2 + if size(proj_h,n) ~= size(proj_l,n) + fprintf('High- and low-energy projections are not the same size!\n') + return; + end +end + +% initialize skip counter for unmatched projections +skipcount = 0; + +for k = 1:size(proj_h,3) + if (~mod(k,50)) + fprintf("Decomposing: projection %d/%d\n", k, size(proj_h,3)); + end + + L_array = proj_interp(proj_l, angles_l, angles_h(k), model, tolerance); + + if(isnan(L_array)) + fprintf("No matches for projection %d. Skipping.\n", k); + skipcount = skipcount+1; + continue; + end + + % map pixels from H and L to equivalent thicknesses + for j = 1:size(proj_h,2) + for i = 1:size(proj_h,1) + H = proj_h(i,j,k); + L = L_array(i,j); + material1_proj(i,j,k-skipcount) = c(1)*L + c(2)*H + c(3)*(L.^2) + c(4)*L.*H + c(5)*(H.^2) + c(6)*(L.^3) + c(7)*(H.*(L.^2)) + c(8)*(L.*(H.^2)) + c(9)*(H.^3); + material2_proj(i,j,k-skipcount) = d(1)*L + d(2)*H + d(3)*(L.^2) + d(4)*L.*H + d(5)*(H.^2) + d(6)*(L.^3) + d(7)*(H.*(L.^2)) + d(8)*(L.*(H.^2)) + d(9)*(H.^3); + angles_out(k-skipcount) = angles_h(k); + end + end +end +end + +function [scanner,calibration_pair, model, tolerance, material1, material2] = parse_inputs(argin) +opts = {'scanner','kVl', 'kVh', 'interpolation', 'tolerance','material1','material2'}; +defaults=ones(length(opts),1); + +% check if option has been passed as input +nVarargs = length(argin); +for ii=1:2:nVarargs + ind=find(ismember(opts,lower(argin{ii}))); + if ~isempty(ind) + defaults(ind)=0; + end +end + +for ii=1:length(opts) + opt=opts{ii}; + default=defaults(ii); + % if one option is not default, then extract value from input + if default==0 + ind=double.empty(0,1);jj=1; + while isempty(ind) + ind=find(isequal(opt,lower(argin{jj}))); + jj=jj+1; + end + + if isempty(ind) + error('TIGRE:DEDecompose:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'scanner' + if default + warning("You have not given a scanner for DE calibration.\n Dual Energy Calibration is strictly dependant on scanner calibration, so do not use the results here for quantitative analysis.\n Assuming VarianTrueBeam") + scanner='VarianTrueBeam'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','Scanner must be a string') + end + scanner=val; + end + case 'kVl' + if default + kVl=80; + else + if ~isscalar(val) + error('TIGRE:DEDecompose:InvalidInput','KVl must be a scalar') + end + kVl=val; + end + case 'kVh' + if default + kVh=140; + else + if ~isscalar(val) + error('TIGRE:DEDecompose:InvalidInput','KVh must be a scalar') + end + kVh=val; + end + case 'interpolation' + if default + model='nearest'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','Interpolation must be a string') + end + model=val; + end + case 'tolerance' + if default + tolerance=1; + else + if ~isnumeric(val) + error('TIGRE:DEDecompose:InvalidInput','tolerance must be a number') + end + tolerance=val; + end + case 'material1' + if default + warning('Assuming Al for material1') + material1='Al'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','material1 must be a string') + end + material1=val; + end + case 'material2' + if default + warning('Assuming PMMA for material2') + material2='PMMA'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','material2 must be a string') + end + material2=val; + end + + otherwise + error('TIGRE:DEDecompose:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in DEDecompose()']); + end +end +calibration_pair=['V_', int2str(kVl), '_' , int2str(kVh)]; +end diff --git a/MATLAB/Utilities/DE/MakeVMproj.m b/MATLAB/Utilities/DE/MakeVMproj.m new file mode 100644 index 00000000..9efd9bff --- /dev/null +++ b/MATLAB/Utilities/DE/MakeVMproj.m @@ -0,0 +1,112 @@ +function [VM_proj] = MakeVMproj(material1_proj, material2_proj, energy, varargin) +%MAKEVMPROJ Synethesizes VM projections from basis material thickness +%projections. +% +% Coded by: Andrew Keeler +% Contact: akeeler@luc.edu +% Date: February 28, 2024 + +[material1, material2] = parse_inputs(varargin); + +% read NIST database (values in TIGRE) +fid = fopen('./../../../Common/data/NIST_material_attenuation.json'); +raw = fread(fid,inf); +str = char(raw'); +fclose(fid); +NIST = jsondecode(str); + +materials = fieldnames(NIST); + +if ~isfield(NIST,material1) + disp([material1, ' not a (yet) supported material in NIST dataset']) + disp("Supported materials:") + for i=1:len(materials) + disp(materials{i}) + end + error("") +end +if ~isfield(NIST,material2) + disp([material2, ' not a (yet) supported material in NIST dataset']) + disp("Supported materials:") + for i=1:len(materials) + disp(materials{i}) + end + error("") +end + +if (energy < max(min(NIST.(material1).energy),min(NIST.(material2).energy)) || energy > min(max(NIST.(material1).energy),max(NIST.(material2).energy))) + warning('Requested energy is outside of expected range. Results may be unreliable.') +end + +% Compute attenuation coefficients for basis materials at selected energy +material1_att = spline(NIST.(material1).energy, NIST.(material1).mu, energy); +material2_att = spline(NIST.(material2).energy, NIST.(material2).mu, energy); + +VM_proj = material1_att*material1_proj + material2_att*material2_proj; +end + + + +function [material1, material2]=parse_inputs(argin) + +opts = {'material1','material2'}; +defaults=ones(length(opts),1); + +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('TIGRE:MakeVMproj:InvalidInput','Invalid number of inputs') +end + +% check if option has been passed as input +for ii=1:2:nVarargs + ind=find(ismember(opts,lower(argin{ii}))); + if ~isempty(ind) + defaults(ind)=0; + end +end + +for ii=1:length(opts) + opt=opts{ii}; + default=defaults(ii); + % if one option is not default, then extract value from input + if default==0 + ind=double.empty(0,1);jj=1; + while isempty(ind) + ind=find(isequal(opt,lower(argin{jj}))); + jj=jj+1; + end + + if isempty(ind) + error('CBCT:MakeVMproj:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'material1' + if default + warning('Assuming Al for material1') + material1='Al'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','material1 must be a string') + end + material1=val; + end + case 'material2' + if default + warning('Assuming PMMA for material2') + material2='PMMA'; + else + if ~ischar( val) + error('TIGRE:DEDecompose:InvalidInput','material2 must be a string') + end + material2=val; + end + + otherwise + error('TIGRE:MakeVMproj:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in MakeVMproj()']); + end +end +end diff --git a/MATLAB/Utilities/DE/demo_de_cal.mat b/MATLAB/Utilities/DE/demo_de_cal.mat new file mode 100644 index 00000000..faa8dc40 Binary files /dev/null and b/MATLAB/Utilities/DE/demo_de_cal.mat differ diff --git a/MATLAB/Utilities/DE/proj_interp.m b/MATLAB/Utilities/DE/proj_interp.m new file mode 100644 index 00000000..4eb4265c --- /dev/null +++ b/MATLAB/Utilities/DE/proj_interp.m @@ -0,0 +1,74 @@ +function proj_out = proj_interp(proj_in, angles_in, target_angle, model, tolerance) +%PROJ_INTERP Performs sinogram interpolation of the input projections +% according to the specified interpolation model +% proj_out is the interpolated projection at target_angle +% +% proj_in is the list of projections to be interpolated, given +% at angles listed in angles_in +% +% model specifies the interpolation model to be used. +% 'linear' uses a linear interpolation model +% 'nearest' uses a nearest-neighbor interpolation model +% 'skip' is the same as 'nearest' but will skip projections without +% sufficiently close matches. Requires setting 'tolerance.' +% +% Coded by: Andrew Keeler +% Contact: akeeler@luc.edu +% Date: February 28, 2024 + +[delta, idx, idx_next] = angle_align(target_angle, angles_in); + +switch model + case 'linear' + dy = proj_in(:,:,idx) - proj_in(:,:,idx_next); + proj_out = proj_in(:,:,idx) + dy*delta; + case 'nearest' + proj_out = squeeze(proj_in(:,:,idx)); + case 'skip' + if delta > tolerance + proj_out = NaN; + return; + end + proj_out = squeeze(proj_in(:,:,idx)); +end + +end + +function [weight, idx, idx_next] = angle_align(target, angle_array) + +angle_diff = angle_array - target; +last = length(angle_diff); +[~, idx] = min(abs(angle_diff)); + +if (angle_diff(idx) == 0) + idx_next = idx; + weight = 1; + return +end + +if (idx == 1 && sign(angle_diff(1)) == sign(angle_diff(2))) + idx_next = idx; + weight = 1; + return +elseif (idx == 1 && sign(angle_diff(1)) ~= sign(angle_diff(2))) + idx_next = 2; + weight = angle_diff(idx)/(abs(angle_array(idx) - angle_array(idx_next))); + return +end + + + +if (idx == last && sign(angle_diff(last)) == sign(angle_diff(last-1))) + idx_next = idx; + weight = 1; + return +end + +if (sign(angle_diff(idx)) ~= sign(angle_diff(idx-1))) + idx_next = idx - 1; + weight = angle_diff(idx)/(abs(angle_array(idx) - angle_array(idx_next))); +else + idx_next = idx + 1; + weight = angle_diff(idx)/(abs(angle_array(idx) - angle_array(idx_next))); +end +end \ No newline at end of file diff --git a/MATLAB/Utilities/TIGRE.bib b/MATLAB/Utilities/TIGRE.bib index 6d87d39e..5b3fc20d 100644 --- a/MATLAB/Utilities/TIGRE.bib +++ b/MATLAB/Utilities/TIGRE.bib @@ -239,7 +239,7 @@ @notabib{ J-H Kim,J Nuyts,A Kyme,Z Kuncic and R Fulton A rigid motion correction method for helical computed tomography (CT) -Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047�2073 +Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047 2073 DOI:10.1088/0031-9155/60/5/2047 @article{kim2015rigid, @@ -257,7 +257,7 @@ @article{kim2015rigid J-H Kim,J Nuyts,A Kyme,Z Kuncic and R Fulton A rigid motion correction method for helical computed tomography (CT) -Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047�2073 +Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047 2073 DOI:10.1088/0031-9155/60/5/2047 @article{kim2015rigid, @@ -275,7 +275,7 @@ @article{kim2015rigid J-H Kim,J Nuyts,A Kyme,Z Kuncic and R Fulton A rigid motion correction method for helical computed tomography (CT) -Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047�2073 +Physics in Medicine & Biology, Phys. Med. Biol. 60 (2015) 2047 2073 DOI:10.1088/0031-9155/60/5/2047 @article{kim2015rigid, @@ -422,5 +422,17 @@ @article{lohvithee2017parameter publisher={IOP Publishing} } +TIGRE-DE + +Keeler A, Lehmann M, Luce J, Kaur M, Roeske J, Kang H. Technical note: TIGRE-DE for the creation of virtual monoenergetic images from dual-energy cone-beam CT. Med Phys. 2024; 1-8. https://doi.org/10.1002/mp.17002 + +@article{https://doi.org/10.1002/mp.17002, +author = {Keeler, Andrew and Lehmann, Mathias and Luce, Jason and Kaur, Mandeep and Roeske, John and Kang, Hyejoo}, +title = {Technical note: TIGRE-DE for the creation of virtual monoenergetic images from dual-energy cone-beam CT}, +journal = {Medical Physics}, +keywords = {cone beam CT, dual-energy, spectral imaging}, +doi = {https://doi.org/10.1002/mp.17002}, +url = {https://aapm.onlinelibrary.wiley.com/doi/abs/10.1002/mp.17002}, +eprint = {https://aapm.onlinelibrary.wiley.com/doi/pdf/10.1002/mp.17002}, End of Bib (leave this here) \ No newline at end of file