Merge pull request #367 from yidu-bjcancer/master

VarianCBCT version 2.0

MATLAB/Utilities/IO/VarianCBCT/BHCalibFromXML.m

@@ -0,0 +1,174 @@
+function [BHCalib, BHCalibXML] = BHCalibFromXML(datafolder, ScanXML)
+%% Load Calibration.xml for Beam Hardening Correction
+% Reference: Improved scatter correction using adaptive scatter kernel superposition
+% Input:
+%               datafolder: Varian CBCT scan data folder
+%               ScanXML: scan geometry structure
+% Output:
+%               BHCalib: key information for Beam Hardening Calibration
+%               BHCalibXML: Beam Hardening Calibration Structure
+% Date: 2021-06-01
+% Author: Yi Du ([email protected])
+
+srcfilename = [datafolder filesep 'Calibrations' filesep 'ASC' filesep 'Factory' filesep 'Calibration.xml'];
+
+%% Export as struct
+tmp = xml2struct(srcfilename);
+BHCalibXML = tmp.Calibration;
+
+%% Restructure BHCalibXML: spectrum data
+% Energy Bin Wid
+BHCalibXML.BHCalib.energybin = str2double(BHCalibXML.CalibrationResults.EnergyBinWidth.Text);
+
+for ii =1:length(BHCalibXML.CalibrationResults.Spectra.SpectrumProperties)
+    % scan voltage
+    tf = strcmpi(ScanXML.Acquisitions.Voltage.Text, BHCalibXML.CalibrationResults.Spectra.SpectrumProperties{ii}.Voltage.Text);
+    if(tf)
+        voltage = str2double(ScanXML.Acquisitions.Voltage.Text);
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Spectra.SpectrumProperties{ii}.Flux.float);
+        for jj = 1:length(tmp)
+            flux(jj) = str2double(tmp(jj).Text);
+        end
+    % voltage
+    BHCalibXML.BHCalib.voltage = voltage;
+    % flux from source
+    BHCalibXML.BHCalib.source.flux = flux;
+    % normalized spectrum from source
+    BHCalibXML.BHCalib.source.spec = flux./max(flux(:));
+    break;
+    end
+end
+
+BHCalibXML.BHCalib.source.kV = 1:BHCalibXML.BHCalib.energybin:BHCalibXML.BHCalib.voltage;
+
+%% Restructure BHCalibXML: Filter type and thickness
+for ii =1:length(BHCalibXML.CalibrationResults.Filters.FilterProperties)
+    % Filter type
+    tf = strcmpi(ScanXML.Acquisitions.KVFilter.Text, BHCalibXML.CalibrationResults.Filters.FilterProperties{ii}.Id.Text);
+    if(tf)
+        % filter name        
+        BHCalibXML.BHCalib.filter.name = BHCalibXML.CalibrationResults.Filters.FilterProperties{ii}.Id.Text;
+        % filer material
+        BHCalibXML.BHCalib.filter.material = BHCalibXML.CalibrationResults.Filters.FilterProperties{ii}.Material.Text;
+        % filer thickness (unit: mm)
+        BHCalibXML.BHCalib.filter.thickness = str2double(BHCalibXML.CalibrationResults.Filters.FilterProperties{ii}.Thickness.Text);
+        break;
+    end        
+end
+
+%% Restructure BHCalibXML: bowtie coordinates and profiles
+% Bowtie type in specific scan is defined in ScanXML
+for ii =1:length(BHCalibXML.CalibrationResults.Bowties.BowtieProperties)
+    % may be some items are comments rather than structured bowtie info
+    if(~isfield(BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}, 'Id'))
+        continue;
+    end
+    tf = strcmpi(ScanXML.Acquisitions.Bowtie.Text, BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.Id.Text);
+    if(tf)
+        % Bowtie Id
+        BHCalibXML.BHCalib.bowtie.name = BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.Id.Text;
+        % Bowtie distance in mm
+        BHCalibXML.BHCalib.bowtie.distance = str2double(BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.SourceDist.Text);
+        if(isfield(BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}, 'Coordinate'))
+            % material
+            BHCalibXML.BHCalib.bowtie.material = BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.Material.Text;
+            % uu, profile
+            tmpuu = cell2mat(BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.Coordinate.float);
+            tmpprofile = cell2mat(BHCalibXML.CalibrationResults.Bowties.BowtieProperties{ii}.Profile.float);
+            for jj = 1:length(tmpuu)
+                coordinates(jj) = str2double(tmpuu(jj).Text);
+                profiles(jj) = str2double(tmpprofile(jj).Text);
+            end
+            % Coordinates: u-direction, mm
+            BHCalibXML.BHCalib.bowtie.uu = coordinates;
+            % thickness: mm
+            BHCalibXML.BHCalib.bowtie.thickness = profiles;
+        end
+    end
+end
+
+%% Restructure BHCalibXML: object material
+% BH calibration material: water in default
+BHCalibXML.BHCalib.object.material = BHCalibXML.CalibrationResults.ObjectMaterialId.Text;
+% maximal thickness for object calibration: mm
+BHCalibXML.BHCalib.object.thicknessmax = str2double(BHCalibXML.CalibrationResults.ObjectThicknessMax.Text);
+% thickness sampling delta: mm
+BHCalibXML.BHCalib.object.delta = str2double(BHCalibXML.CalibrationResults.ObjectThicknessDelta.Text);
+
+%% Restructure BHCalibXML:  material coefficient database
+% attenuation coeeficient: /mm
+for ii =1:length(BHCalibXML.CalibrationResults.Materials.MaterialProperties)
+    % Filter material
+    filter_tf = strcmpi(BHCalibXML.BHCalib.filter.material, BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text);
+    if(filter_tf)
+        ac = [];
+        % attenuation coefficients: /mm
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.AttenuationCoefficient.float);
+        for jj = 1:length(tmp)
+            ac(jj) = str2double(tmp(jj).Text);
+        end
+        BHCalibXML.BHCalib.filter.ac = ac;
+    end
+
+    % Object material
+    obj_tf = strcmpi(BHCalibXML.BHCalib.object.material, BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text);
+    if(obj_tf)
+        ac = [];
+        % attenuation coefficients: /mm
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.AttenuationCoefficient.float);
+        for jj = 1:length(tmp)
+            ac(jj) = str2double(tmp(jj).Text);
+        end
+        BHCalibXML.BHCalib.object.ac = ac;
+    end
+
+    % Bowtie material
+    bowtie_tf = strcmpi(BHCalibXML.BHCalib.bowtie.material, BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text);
+    if(bowtie_tf)
+        ac = [];
+        % attenuation coefficients: /mm
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.AttenuationCoefficient.float);
+        for jj = 1:length(tmp)
+            ac(jj) = str2double(tmp(jj).Text);
+        end
+        BHCalibXML.BHCalib.bowtie.ac = ac;
+    end
+
+    % Bone material for further correction
+    bone_tf = contains(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text, 'bone', 'IgnoreCase',true);
+    if(bone_tf)
+        ac = [];
+        % bone material name
+        BHCalibXML.BHCalib.bone.material = BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text;   
+        % attenuation coefficients: /mm
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.AttenuationCoefficient.float);
+        for jj = 1:length(tmp)
+            ac(jj) = str2double(tmp(jj).Text);
+        end
+        BHCalibXML.BHCalib.bone.ac = ac;
+    end
+
+    % Scitillator material for detection response (PaxScan 4030CB: CsI: Ti)
+    scintillator_tf = contains(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text, 'cesium', 'IgnoreCase', true);
+    if(scintillator_tf)
+        ac = [];
+        % detector scintillator material
+        BHCalibXML.BHCalib.scintillator.material = BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.Id.Text;   
+
+        % CsI thickness: this parameter is not important at all.
+        BHCalibXML.BHCalib.scintillator.thickness = 0.6;   
+
+        % attenuation coefficients: /mm
+        tmp = cell2mat(BHCalibXML.CalibrationResults.Materials.MaterialProperties{ii}.EnergyAbsorptionCoefficient.float);
+        for jj = 1:length(tmp)
+            ac(jj) = str2double(tmp(jj).Text);
+        end
+        % energy absorption coefficient
+        BHCalibXML.BHCalib.scintillator.ac = ac;
+    end
+
+end
+
+%% Key Information for further beam hardening calibration
+BHCalib = BHCalibXML.BHCalib;
+end

MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m

@@ -0,0 +1,36 @@
+function [proj_lg, BHCalib] = BHCorrection(datafolder, geo, ScanXML, proj_lg,gpuids)
+% Entry Function For BH Correction
+%   Detailed explanation goes here
+
+
+disp('BH correction is very memory intensive. The effect is to be justified. ');
+typein = input('Continue BH correction? Y or N (recommended): ', 's');
+if(contains(typein, 'n', 'IgnoreCase', true))
+    BHCalib = NaN;
+    disp('BH correction is skipped.');
+    return;
+end
+
+disp('Beam Hardening Correction is on-going: be patient... ');
+
+% Key calibration information
+BHCalib = BHCalibFromXML(datafolder, ScanXML);
+
+% Precompute filter attenuated spectrum: debug pass
+BHCalib = BH_SpectrumFilter(BHCalib);
+
+% Precompute bowtie attenuated spectra LUT
+BHCalib = BH_SpectrumBowtieLUT(geo, BHCalib);
+
+% Build reference object (water) attanuation LUT
+BHCalib = BH_ObjectCalibLUT(BHCalib);
+
+% BH correction via reference object (water)
+BHCalib = BH_RemappingFunc(BHCalib);
+
+proj_lg = BH_ObjectRemapping(BHCalib, proj_lg, gpuids);
+
+disp('BH correction is done.')
+
+end
+

MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m

@@ -0,0 +1,50 @@
+function BHCalib = BH_ObjectCalibLUT(BHCalib)
+%BHWATERCALIB Summary of this function goes here
+%   Detailed explanation goes here
+
+% max object thickness: mm
+MaxThickness = BHCalib.object.thicknessmax;
+
+% object sampling length: mm
+object_sl = linspace(0, MaxThickness, ceil(MaxThickness)*10);
+
+%% spectra[sI, energy]
+specLUT = BHCalib.bowtie.specLUT;
+miu = BHCalib.object.ac;
+
+% object thickness look-up table: [bowtie.sl, object_sl]
+calibLUT = zeros(length(BHCalib.bowtie.sl), length(object_sl));
+
+%% Energy Flux based calibration
+% scintillator energy absorption coefficient
+scintillator_ac = BHCalib.scintillator.ac(1: length(miu));
+scintillator_thick = BHCalib.scintillator.thickness;
+
+disp('Calculating Ojbect Look-up Table: ')
+% bowtie thickness - > spectrum
+for ii = 1: length(BHCalib.bowtie.sl)    
+    % object thickness -> [min_bowtie_thickness: max_bowtie_thickness]
+    spec = specLUT(ii,:);
+
+    % incident energy fluence
+    fluence_0 = sum(spec .* BHCalib.source.kV .* scintillator_ac .* scintillator_thick);
+
+    for jj = 1:length(object_sl)
+        % attenuated spectram 
+        tmp = spec .* exp(-object_sl(jj).* miu);
+
+        % attenuated energy fluence
+        fluence_1 = sum(tmp .*BHCalib.source.kV .*scintillator_ac.*scintillator_thick);
+
+        % logarithmic projection signal: 
+        calibLUT(ii,jj) = -log(fluence_1/fluence_0);
+    end
+end
+%% non ideal projection signal LUT: [bowtie.sl, object_sl]
+BHCalib.object.calibLUT = calibLUT;
+
+% object sampling length
+BHCalib.object.sl = object_sl;
+
+end
+

MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m

@@ -0,0 +1,102 @@
+function proj_BH = BH_ObjectRemapping(BHCalib, projlg, gpuids)
+%BH_OBJECTREMAPPING Summary of this function goes here
+%   Detailed explanation goes here
+% lgproj: log normalized projection (non-ideal BH projection)
+
+
+% bowtie tranvers length grid
+ulgd = BHCalib.bowtie.ulgd;
+[nRow, nCol] = size(ulgd);
+
+% Object Calibration LUT: [bowtie_sl, object_sl]
+calibLUT = BHCalib.object.calibLUT;
+
+% Object Sampling Length
+object_sl = BHCalib.object.sl;
+
+%% CPU Version: very very slow
+%{
+proj_BH = zeros(size(projlg));
+
+for ii = 1: nRow
+    if(~mod(ii,50))
+        disp(ii);
+    end
+    for jj = 1: nCol
+        tl_v = ulgd(ii, jj);
+        % model the LUT for specific bowtie thickness
+        proj_signal_LUT = interp1(BHCalib.bowtie.sl, calibLUT, tl_v, 'spline');
+
+    end
+end
+%}
+
+%% Parfor CPU version: very slow
+%{
+proj_BH = zeros(size(projlg));
+
+sl = BHCalib.bowtie.sl;
+tic
+for ii = 1: nRow
+    if(~mod(ii,50))
+        disp(ii);
+        toc
+        tic
+    end
+    parfor jj = 1: nCol
+        tl_v = ulgd(ii, jj);
+        % model the LUT for specific bowtie thickness
+        proj_signal_LUT = interp1(sl, calibLUT, tl_v, 'spline');
+        proj_BH(ii,jj,:) = interp1(proj_signal_LUT, object_sl, projlg(ii,jj,:), 'spline');
+    end
+end
+%}
+
+%% GPU version: acceptable
+% GPU Reset
+g = gpuDevice(gpuids.devices(0)+1);
+reset(g);
+
+gprojlg = gpuArray(single(projlg));
+
+gsl = gpuArray(BHCalib.bowtie.sl);
+gLUT = gpuArray(calibLUT);
+
+gobj_sl = gpuArray(object_sl);
+proj_signal_LUT = gpuArray(zeros(1, size(gLUT,2)));
+
+% 1D interpolation pixel-by-pixel
+for ii = 1: nRow
+    if(~mod(ii,50))
+        disp([num2str(ii) '/' num2str(nRow)]);
+    end
+    for jj = 1: nCol
+        tl_v = ulgd(ii, jj);
+        % model the LUT for specific bowtie thickness: BUG! in interp1 for
+        % spline interpolation
+        % proj_signal_LUT = interp1(gsl, gLUT, tl_v, 'spline');
+        proj_signal_LUT = interp1(gsl, gLUT, tl_v);
+
+        eqv_thickness = interp1(proj_signal_LUT, gobj_sl, gprojlg(ii,jj,:));
+
+        % fitting coeficiency
+        ft_cof = interp1(gsl, BHCalib.object.linear_ft_cof, tl_v);
+        gprojlg(ii, jj, :) = ft_cof.* eqv_thickness;
+
+    end
+end
+
+proj_BH = double(gather(gprojlg));
+
+proj_BH(proj_BH<0) = NaN;
+for ii = 1:size(proj_BH,3)
+    % default fast extrapolation: robust for noise and holes at boundaries
+    proj_BH(:,:,ii) = inpaint_nans(double(proj_BH(:,:,ii)), 2);
+end
+proj_BH(proj_BH<0) = eps;
+
+% GPU reset
+reset(g);
+
+end
+