From e96e9e2b657ea956d3b4c6709f2021594ef1671d Mon Sep 17 00:00:00 2001 From: Yi Du Date: Fri, 20 May 2022 18:15:32 +0800 Subject: [PATCH 1/4] VarianCBCT version 2.0 --- .../Utilities/IO/VarianCBCT/BHCalibFromXML.m | 174 ++++++++++++++++++ MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m | 36 ++++ .../IO/VarianCBCT/BH_ObjectCalibLUT.m | 65 +++++++ .../IO/VarianCBCT/BH_ObjectRemapping.m | 118 ++++++++++++ .../IO/VarianCBCT/BH_RemappingFunc.m | 55 ++++++ .../IO/VarianCBCT/BH_SpectrumBowtieLUT.m | 115 ++++++++++++ .../BH_SpectrumBowtieLUT_deprecated.m | 49 +++++ .../IO/VarianCBCT/BH_SpectrumFilter.m | 24 +++ MATLAB/Utilities/IO/VarianCBCT/BlkLoader.m | 68 +++++++ MATLAB/Utilities/IO/VarianCBCT/CTDI.m | 40 ++++ .../DetectorPointScatterCorrection.m | 111 +++++++++++ .../Utilities/IO/VarianCBCT/EnforcePositive.m | 10 + .../Utilities/IO/VarianCBCT/GeometryFromXML.m | 12 +- MATLAB/Utilities/IO/VarianCBCT/HUMapping.m | 27 +++ MATLAB/Utilities/IO/VarianCBCT/LogNormal.m | 75 ++++++++ .../Utilities/IO/VarianCBCT/PolarFiltering.m | 44 +++++ MATLAB/Utilities/IO/VarianCBCT/ProjLoader.m | 75 ++++++++ MATLAB/Utilities/IO/VarianCBCT/RingRemoval.m | 70 +++++++ MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m | 53 ++++++ .../IO/VarianCBCT/SC_AmplitudeFactor.m | 35 ++++ .../Utilities/IO/VarianCBCT/SC_EdgeResponse.m | 23 +++ MATLAB/Utilities/IO/VarianCBCT/SC_FormFunc.m | 32 ++++ MATLAB/Utilities/IO/VarianCBCT/SC_GroupMask.m | 30 +++ .../IO/VarianCBCT/SC_SmoothThickness.m | 28 +++ .../IO/VarianCBCT/SC_ThicknessEstimator.m | 32 ++++ .../Utilities/IO/VarianCBCT/ScCalibFromXML.m | 44 +++++ .../IO/VarianCBCT/ScatterCorrection.m | 153 +++++++++++++++ .../IO/VarianCBCT/VarianDataLoader.m | 114 ++++++++---- MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp | 9 +- .../Utilities/IO/VarianCBCT/interp_weight.m | 35 ++++ MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp | 33 +++- .../IO/VarianCBCT/polar2cart/ImToPolar.m | 61 ++++++ .../IO/VarianCBCT/polar2cart/PolarToIm.m | 71 +++++++ .../IO/VarianCBCT/polar2cart/license.txt | 24 +++ MATLAB/Utilities/im2DDenoise.m | 24 +++ 35 files changed, 1923 insertions(+), 46 deletions(-) create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BHCalibFromXML.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_RemappingFunc.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT_deprecated.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumFilter.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/BlkLoader.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/CTDI.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/EnforcePositive.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/HUMapping.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/LogNormal.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/PolarFiltering.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/ProjLoader.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/RingRemoval.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_EdgeResponse.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_FormFunc.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_GroupMask.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_SmoothThickness.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/SC_ThicknessEstimator.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/ScCalibFromXML.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/interp_weight.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/polar2cart/ImToPolar.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/polar2cart/PolarToIm.m create mode 100644 MATLAB/Utilities/IO/VarianCBCT/polar2cart/license.txt create mode 100644 MATLAB/Utilities/im2DDenoise.m diff --git a/MATLAB/Utilities/IO/VarianCBCT/BHCalibFromXML.m b/MATLAB/Utilities/IO/VarianCBCT/BHCalibFromXML.m new file mode 100644 index 00000000..b60257c1 --- /dev/null +++ b/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 (yi.du@hotmail.com) + +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 diff --git a/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m new file mode 100644 index 00000000..a9155039 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m @@ -0,0 +1,36 @@ +function [proj_lg, BHCalib] = BHCorrection(datafolder, geo, ScanXML, proj_lg) +% 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); + +disp('BH correction is done.') + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m new file mode 100644 index 00000000..bcd8d83e --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m @@ -0,0 +1,65 @@ +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)); + +%% Photo flux based calibration +%{ +% bowtie thickness - > spectrum +for ii = 1: length(BHCalib.bowtie.sl) + % object thickness -> [min_bowtie_thickness: max_bowtie_thickness] + spec = specLUT(ii,:); + for jj = 1:length(object_sl) + % non-ideal attenuated signal + tmp = spec .* exp(-object_sl(jj).* miu); + % nonlinear BH projection + calibLUT(ii,jj) = -log(sum(tmp)/sum(spec)); + end +end +%} + +%% 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 + diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m new file mode 100644 index 00000000..3fc5c691 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m @@ -0,0 +1,118 @@ +function proj_BH = BH_ObjectRemapping(BHCalib, projlg) +%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; + +%% Remap non-ideal BH signal to linear object thickness +%{ +proj_BH = zeros(size(projlg)); + +% for 1D bowtie only +for jj = 1: nCol + tl_v = ulgd(1, jj); + % model the LUT for specific bowtie thickness + proj_signal_LUT = interp1(BHCalib.bowtie.sl, calibLUT, tl_v, 'spline'); + proj_BH(:,jj,:) = interp1(proj_signal_LUT, object_sl, projlg(:,jj,:), 'spline'); +end +%} + +%% 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(1); +reset(g); + +% gpuArray preparation +% proj_BH = single(zeros(size(projlg))); +% gproj_BH = zeros(size(proj_BH), 'gpuArray'); +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 +g = gpuDevice(1); +reset(g); + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_RemappingFunc.m b/MATLAB/Utilities/IO/VarianCBCT/BH_RemappingFunc.m new file mode 100644 index 00000000..c747be8e --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_RemappingFunc.m @@ -0,0 +1,55 @@ +function BHCalib = BH_RemappingFunc(BHCalib) +% remaps the non-linear projection signals via linear fitted models +% +% SYNOPSIS: BHCalib = BH_RemappingFunc(BHCalib) +% +% INPUT BH Calibration Structure +% +% OUTPUT BHCalib linear mapping functions +% +% REMARKS +% +% created with MATLAB ver.: 9.8.0.1873465 (R2020a) Update 8 on Microsoft Windows 10 Pro Version 10.0 (Build 18363) +% +% created by: Yi Du +% DATE: 16-May-2022 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +bowtie_sl = length(BHCalib.bowtie.sl); +object_sl = length(BHCalib.object.sl); + +% control percentage for linear fitting: 5% (default) +cp = 0.2; + +fitted_sl = round(object_sl *cp); +object_sl = BHCalib.object.sl(1:fitted_sl); + +% Set up fittype and options. +ft = fittype( 'a*x', 'independent', 'x', 'dependent', 'y'); +opts = fitoptions( 'Method', 'NonlinearLeastSquares' ); +opts.StartPoint = 0.01; + +cof = zeros(bowtie_sl,1); +rsquare = cof; + +%% Calculate Fitting Parameters +disp('Calculating Linear Fitting Model: ') +for ii = 1:bowtie_sl + [xData, yData] = prepareCurveData(object_sl, BHCalib.object.calibLUT(ii, 1:fitted_sl)); + % Fit model to data. + [func_ft, gof] = fit( xData, yData, ft, opts ); + cof(ii) = func_ft.a; + rsquare(ii) = gof.rsquare; +end + +%% Fitting Model to Return +BHCalib.object.linear_ft_cof = cof; + +%{ +plot(BHCalib.object.sl, BHCalib.object.calibLUT(ii,:), 'r'), grid on; hold on; +plot(BHCalib.object.sl,cof(ii).*BHCalib.object.sl, 'b'), grid on; +plot(BHCalib.object.sl,cof(ii).*BHCalib.object.sl - BHCalib.object.calibLUT(ii,:), 'k'), grid on; +%} + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m new file mode 100644 index 00000000..62c0f5e3 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m @@ -0,0 +1,115 @@ +function BHCalib = BH_SpectrumBowtieLUT(geo, BHCalib) +%SPECTRAPROCESS Summary of this function goes here +% Detailed explanation goes here + +%% Thin window filter has been applied to the spectrum already +% Source-to-detector distance: mm +SDD = geo.DSD; +% Source-to-bowtie distance: mm +SBD = BHCalib.bowtie.distance; + +%% [u,v] vector: mm +offset=geo.offDetector; +us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + offset(1); +vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2) + offset(2); + +%% Map Detector Plane to Bowtier Surface +ub = us * SBD/SDD; +% vb = vs * SBD/SDD; + +u_s = fix(min(ub)) + sign(min(us))*10; +u_l = fix(max(ub)) + sign(max(us))*10; + +sampling_rate = 4002; + +xxgd = linspace(u_s, u_l, sampling_rate); +y_l = ceil(max(BHCalib.bowtie.thickness)); +yygd = linspace(-y_l, y_l, sampling_rate); +% bowtie height = 55 mm, ref to TrueBeam Technical Reference +zzgd = linspace(-27.5, 27.5, 102); + +bt_img = zeros(length(xxgd), length(yygd)); + + +%% Fitting to Model Bowtie Profile: +[xData, yData] = prepareCurveData(BHCalib.bowtie.uu, BHCalib.bowtie.thickness); +%{ +%% Fitting Model: first fitting(linear) +% Set up fittype and options. +ft = fittype( 'poly1' ); +% Fit model to data. +[~, gof] = fit( xData, yData, ft ); +if( gof.rsquare < 0.995) + disp('It seems that BH correction is not required at all'); + typein = input('Continue BH correction? Y or N (recommended): ', 's'); + if(contains(typein, 'n', 'IgnoreCase', true)) + BHCalib = NaN; + return; + end +end +%} +%% Fitting Model: Smoothing Spline +% Set up fittype and options. +ft = fittype( 'smoothingspline' ); + +% Fit model to data. +[fitresult, gof] = fit(xData, yData, ft); +disp(['Bowtie Fitting Goodness (R^2): ' num2str(gof.rsquare)]) + +%% 3D bowtie img matrix +for ii = 1:length(xxgd) + tmp_thick = fitresult(xxgd(ii)); + tmp = ((yygd>0)&(yygd0) + RtnDirection = 'CC'; +else + RtnDirection = 'CW'; +end + +%% TrueBeam Version +if( strfind(ScanXML.Attributes.Version, '2.7')) + Version = 2.7; +elseif( strfind(ScanXML.Attributes.Version, '2.0')) + Version = 2.0; +else + error('Error in Version'); +end + +%% Bowtie +tag_bowtie = strfind(ScanXML.Acquisitions.Bowtie.Text, 'None'); + +%% Blk File Struct +% wi Bowtie1 +if(isempty(tag_bowtie)) + % TB 2.0 + if(Version == 2.0) + blkfilestr = dir([datafolder filesep 'Calibrations' filesep 'AIR-*' filesep '**' filesep 'FilterBowtie.xim']); + elseif(Version == 2.7) + blkfilestr = dir([datafolder filesep 'Calibrations' filesep 'AIR-*' filesep '**' filesep 'FilterBowtie_' RtnDirection '*.xim']); + end +else + blkfilestr = dir([datafolder filesep 'Calibrations' filesep 'AIR-*' filesep '**' filesep 'Filter.xim']); +end + +%% Load Blk and AirNorm for TB2.5 +if(Version == 2.0) + filename = fullfile(blkfilestr.folder, blkfilestr.name); + tmp = mexReadXim(filename); + Blk = rot90(single(tmp), -1); + tmp = ReadXim(filename, 0); + BlkAirNorm = single(tmp.properties.KVNormChamber); + Sec = []; + return; +end + +%% Cresent Artifacts correction +if(Version == 2.7) + for ii = 1:length(blkfilestr) + filename = fullfile(blkfilestr(ii).folder, blkfilestr(ii).name); + tmp = mexReadXim(filename); + Blk(:,:,ii) = rot90(single(tmp), -1); + tmp = ReadXim(filename, 0); + BlkAirNorm(ii) = single(tmp.properties.KVNormChamber); + % GantryRtn + Sec(ii) = single(tmp.properties.GantryRtn); + end +end + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/CTDI.m b/MATLAB/Utilities/IO/VarianCBCT/CTDI.m new file mode 100644 index 00000000..811ffdfe --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/CTDI.m @@ -0,0 +1,40 @@ +function CTDI = CTDI(datafolder) +% Calculate CT Dosimetry Index (cGy) of a CBCT scan +% CTDI = DoseFactor * (mAs/100) * Num_Projection +% Reference: TrueBeam Technical Reference Guide II: Imaging +% Date: 2021-11-23 +% Author: Yi Du +% Email: yi.du@hotmail.com + +% Load ScanXML +if(~exist('ScanXML', 'var')) + % Read ScanXML + filestr = dir([datafolder filesep 'Scan.xml']); + ScanXML = getfield(xml2struct(fullfile(filestr.folder, filestr.name)), 'Scan'); +end + +% Dose Factor: cGy /100 mAs +DoseFc = str2double(ScanXML.Acquisitions.DoseFactor.Text); +mA = str2double(ScanXML.Acquisitions.Current.Text); +mS = str2double(ScanXML.Acquisitions.PulseLength.Text); +% mAs +mAs = mA*mS/1000; + +filedir = dir([datafolder filesep 'Acquisitions' filesep '*' filesep 'Proj_*']); + +% FileSize Threshold 10 kB to Exclude Invalid Acquisition +threshold = 10*1024; +filesize = zeros(length(filedir),1); + +for ii=1:length(filedir) + filesize(ii) = (filedir(ii).bytes >threshold); +end + +% Projection Number +nProj = sum(filesize); + +% CTDIw, cGy +CTDI = DoseFc * (mAs/100) *nProj; + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m new file mode 100644 index 00000000..d1e61704 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m @@ -0,0 +1,111 @@ +function proj = DetectorPointScatterCorrection(proj, geo, ScCalib) +%% Detector Point Scatter Correction +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Date: 2021-03-26 +% Author: Yi Du (yi.du@hotmail.com) + +%% Empirical values from reference paper +% unit: cm-2 +% a0 = 3.43; (in refrence paper, but incorrect) +a0 = 1; + +a1 = 0.000309703536035; +a1 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit0.Text); + +% unit: cm-1 +a2 = 0.546566915157327; +a2 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit1.Text); + +% unit: cm-1 +a3 = 0.311272841141691; +a3 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit2.Text); + +% unit: cm-1 +a4 = 0.002472148007134; +a4 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit3.Text); + +a5 = -12.6606856375944; +a5 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit4.Text); + +% for amplitude normalization +CoverSPR = 0.04; + +% unit: mm +offset=geo.offDetector; + +% grid unit: mm +us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1); % + offset(1); +vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2); % + offset(2); +% unit mm - > cm +% unit converter: 1/10 for mm-> cm +mm2cm = 1/10; +us = us * mm2cm; +vs = vs * mm2cm; + +%% Downsampling +% about 10 mm in axial direction +dus = downsample(us, 26); +% about 4 mm in axial direction +dvs = downsample(vs, 10); + +ds_rate = 8; +dus = decimate(us, ds_rate); +dvs = decimate(vs, ds_rate); + + +%% Grid mesh +[uu,vv] = meshgrid(us,vs); %detector +[duu, dvv] = meshgrid(dus,dvs); %detector + +%% Scatter convolution kernel +grid = sqrt(duu.^2 + dvv.^2); +% a0 is the normalization factor +hd = a0*( a1* exp(-a2 * grid) + a3 * (exp( -a4 * ( grid - a5).^3 ))); + +% a0 = CoverSPR/sum(hd(:)); + +% normalized to 4% SPR coverage +hd = CoverSPR/sum(hd(:)) .* hd; + +%% GPU based +% reset(gpuDevice(1)); +gproj = gpuArray(single(proj)); + +%% 2D Convolution with downsampling and upsampling +for ii = 1:size(proj,3) + % CPU version + %{ + page = interp2(uu, vv, proj(:,:,ii), duu, dvv); + % gross scatter distribution + sc = conv2(page, hd, 'same'); + % upsample the scatter distribution to the same grid level as the + % measured intensity + scpage = interp2(duu, dvv, sc, uu, vv, 'spline'); + % primary = measure - scatter + proj(:,:,ii) = proj(:,:,ii) - scpage; + %} + + % GPU version + page = interp2(uu, vv, gproj(:,:,ii), duu, dvv); + % gross scatter distribution + sc = gather(conv2(page, hd, 'same')); + % upsample the scatter distribution to the same grid level as the + % measured intensity + scpage = interp2(duu, dvv, sc, uu, vv, 'spline'); + % primary = measure - scatter + gproj(:,:,ii) = gproj(:,:,ii) - scpage; +end + +proj = single(gather(gproj)); +% Reset GPU +reset(gpuDevice(1)); + +%% Cutoff for over-correction +proj(proj<0) = NaN; +for ii = 1:size(proj,3) + % default fast extrapolation: robust for noise and holes at boundaries + proj(:,:,ii) = single(inpaint_nans(double(proj(:,:,ii)), 2)); +end +proj(proj<0) = eps; + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/EnforcePositive.m b/MATLAB/Utilities/IO/VarianCBCT/EnforcePositive.m new file mode 100644 index 00000000..b16fee78 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/EnforcePositive.m @@ -0,0 +1,10 @@ +function proj = EnforcePositive(proj) +%% Remove anomalies +% in case of NaN or Inf +proj(isnan(proj)) = 0; +proj(isinf(proj)) = 0; + +%% Set all negative to zeros +proj(proj<0) = 0; + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/GeometryFromXML.m b/MATLAB/Utilities/IO/VarianCBCT/GeometryFromXML.m index dcac45fb..0dddc487 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/GeometryFromXML.m +++ b/MATLAB/Utilities/IO/VarianCBCT/GeometryFromXML.m @@ -1,6 +1,11 @@ function [geo, ScanXML, ReconXML] = GeometryFromXML(datafolder,varargin) +% Purpose: +% 1) load scan and reconstruction geometry from .xml files +% 2) generate geo for TIGRE % Dependence: xml2struct.m % Date: 2020-03-28 +% Author: Yi Du, yi.du@hotmail.com + if nargin>1 load_geo=varargin{1}~=0; else @@ -24,7 +29,7 @@ % Cone-beam CT scenario geo.mode = 'cone'; -% in geo, only 14 predefined fields can allowed +% in geo, only 14 predefined fields are allowed, hard-coded in TIGRE % determine arc/fan in FDK %{ % Circular Trajectory: Full/Half @@ -57,9 +62,10 @@ % Offset of Detector (mm) offset = str2double(ScanXML.Acquisitions.ImagerLat.Text); -geo.offDetector = [offset; 0 ]; +% offset orientation in VarianCBCT is opposite to TIGRE +geo.offDetector = [offset * (-1); 0 ]; % Offset of image from origin (mm) -geo.offOrigin =[0;0;0]; +geo.offOrigin =[0;0;0]; % Auxiliary % Variable to define accuracy of 'interpolated' projection diff --git a/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m b/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m new file mode 100644 index 00000000..f3eaee7b --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m @@ -0,0 +1,27 @@ +function HU3D = HUMapping(img3D, Protocol) +% Pixel value to HU value Mapping +% Method: using CATPhan 604 calibration data +% Input: +% img3D: reconstructed image matrix in pixel value +% Protocol: CBCT scan protocol +% Output: +% HU3D: image matrix in CT HU value +% Date: 2021-09-03 +% Author: Yi Du, yi.du@hotmail.com + +HU3D = zeros(numel(img3D), 1); +load('D:\MATLABWorkplace\Ongoing_Projects\Fronc_VarianCBCT\demo_Pixel2HU.mat', 'Pixel2HU'); + +tmp = getfield(Pixel2HU, Protocol); + +[xData, yData] = prepareCurveData(tmp(:,1), tmp(:,2)); + +% 2nd-order Polynominal Fitting :Set up fittype and options. +ft = fittype( 'poly2' ); + +% Fit model to data. +[fitresult, gof] = fit( xData, yData, ft); + +HU3D = reshape(fitresult(img3D), size(img3D)); + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m b/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m new file mode 100644 index 00000000..33203747 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m @@ -0,0 +1,75 @@ +function proj = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm) +% Log Normalization: +% Calculate Logrithmic Projections with AirNorm +% Tested on TrueBeam 2.5 and 2.7 +% Date: 2021-03-23 +% Author: Yi Du (yi.du@hotmail.com) + +%% TB version: CPU +%{ +% +tic +% Version = 2.5 +if(isempty(Sec)) + for ii = 1:length(angles) + CF = airnorm(ii)/BlkAirNorm; + proj(:,:,ii) = log(CF*Blk./(proj(:,:,ii) + eps) + eps); + end +% Version = 2.7 +else + % GantryRtn = KvSourceRtn - 90; + angles = angles - 90; + + % Flip for interpolation + if(Sec(2) - Sec(1) <0) + Sec = flip(Sec); + BlkAirNorm = flip(BlkAirNorm); + Blk = flip(Blk, 3); + end + + % interpolation weights + for ii = 1:length(angles) + [left, weights] = interp_weight(angles(ii), Sec); + interp_blk = weights(1) * Blk(:,:,left) + weights(2) * Blk(:,:,left+1); + % Correction factor + CF = airnorm(ii)/(0.5*BlkAirNorm(left) + 0.5*BlkAirNorm(left+1)); + proj(:,:,ii) = log(CF*interp_blk./(proj(:,:,ii)+eps) + eps); + end +end +toc +%} +%% TB version: GPU +gproj = gpuArray(single(proj)); +gBlk = gpuArray(Blk); +% Version = 2.5 +if(isempty(Sec)) + for ii = 1:length(angles) + CF = airnorm(ii)/BlkAirNorm; + gproj(:,:,ii) = log(CF*gBlk./(gproj(:,:,ii) + eps) + eps); + end +% Version = 2.7 +else + % GantryRtn = KvSourceRtn - 90; + angles = angles - 90; + + % Flip for interpolation + if(Sec(2) - Sec(1) <0) + Sec = flip(Sec); + BlkAirNorm = flip(BlkAirNorm); + gBlk = flip(gBlk, 3); + end + + % interpolation weights + for ii = 1:length(angles) + [left, weights] = interp_weight(angles(ii), Sec); + interp_blk = weights(1) * gBlk(:,:,left) + weights(2) * gBlk(:,:,left+1); + % Correction factor + CF = airnorm(ii)/(0.5*BlkAirNorm(left) + 0.5*BlkAirNorm(left+1)); + gproj(:,:,ii) = log(CF*interp_blk./(gproj(:,:,ii)+eps) + eps); + end +end + +proj = gather(gproj); +% GPU clear +reset(gpuDevice(1)); +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/PolarFiltering.m b/MATLAB/Utilities/IO/VarianCBCT/PolarFiltering.m new file mode 100644 index 00000000..c8d58cdb --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/PolarFiltering.m @@ -0,0 +1,44 @@ +function img3D_f = PolarFiltering(img3D) +% Polar Filtering in image domain for ring correction +% Method: dependent on third-party tool, polar2cart +% ring artifacts is first to be identified, and only one ring can be +% corrected in each run. +% Input: +% img3D: image matrix +% geo: geometry +% Output: +% img: filtered projection matrix +% Date: 2022-02-04 +% Author: Yi Du, yi.du@hotmail.com + + +imgSize = size(img3D, 1); +img3D_f = img3D; + +%% slice by slice +for ii = 1:size(img3D,3) + % Cartisan to Polar + imP = ImToPolar(img3D(:,:,ii), 0, 1, imgSize, imgSize); + sum_over_col = sum(imP,2); + + % find the indices of max and min values + [~, max_idx] = max(diff(sum_over_col), [], 'linear'); + [~, min_idx] = min(diff(sum_over_col), [], 'linear'); + + % median filtering + if(abs(max_idx - min_idx)>10) + continue; + else + if(max_idx > min_idx) + imP(min_idx:max_idx,:) = ordfilt2(imP(min_idx:max_idx,:), 5, ones(5,1)); + else + imP(max_idx:min_idx,:) = ordfilt2(imP(max_idx:min_idx,:), 5, ones(5,1)); + end + end + + img3D_f(:,:,ii) = PolarToIm(imP, 0, 1, imgSize, imgSize); + +end + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/ProjLoader.m b/MATLAB/Utilities/IO/VarianCBCT/ProjLoader.m new file mode 100644 index 00000000..983ba7f1 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/ProjLoader.m @@ -0,0 +1,75 @@ +function [proj, angles, airnorm] = ProjLoader(datafolder,varargin) +% [proj, angle, blk, projinfo] = BatchReadXim(foldername, varargin) +% Varian .xim files exported from TrueBeam on-board imager/EPID, Halcyon imager +% .xim files are exported under service mode, as a single .zip file +% Once unzipped, all relevant files are organized under the same root folder +% Input: +% datafolder : root folder of the upzipped files, +% or the subfolder Acquisitions where all xim files are stored +% thd : 0(default), no motion lag correcion +% thd, motion correcion threadhold +% Output: +% proj : .xim pixel_images frame by frame +% angles : KVSourceRtn at each frame +% projinfo : proj_*.xim info cell arrays +% Date: 2021-03-23 +% Author: Yi Du, yi.du@hotmail.com + +thd = 0; +if(nargin == 2) + thd = varargin{1}; +end + +fprintf("Dataset in: %s \n", datafolder); + +ximfilelist = dir([datafolder filesep '**' filesep 'Proj_*.xim']); + +proj = []; +angles = []; +airnorm = []; + +% +proj_no = length(ximfilelist); +% +count = 1; +for ii = 1:proj_no + ximfilename = fullfile(ximfilelist(ii).folder, ximfilelist(ii).name); + if(~mod(ii,50)) + fprintf("Loading: %d / %d \n", ii, proj_no); + end + [page, rtn] = mexReadXim(ximfilename); + if(~isempty(page)) + % load first page + if(count==1) + projinfo{count} = ReadXim(ximfilename, 0); + angles(count) = rtn; + proj(:,:,count) = rot90(single(page), -1); + + airnorm(count) = single(projinfo{count}.properties.KVNormChamber); + count = count + 1; + else + % proj info + projinfo{count} = ReadXim(ximfilename, 0); + + % whether threshold is applied + if(thd) + if(abs(rtn-angles(end))>thd) + angles(count) = rtn; + proj(:,:,count) = rot90(single(page), -1); + airnorm(count) = single(projinfo{count}.properties.KVNormChamber); + count = count + 1; + else + continue; + end + else + angles(count) = rtn; + proj(:,:,count) = rot90(single(page), -1); + airnorm(count) = single(projinfo{count}.properties.KVNormChamber); + count = count + 1; + end + end + end +end +fprintf("Loading complete \n"); + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/RingRemoval.m b/MATLAB/Utilities/IO/VarianCBCT/RingRemoval.m new file mode 100644 index 00000000..3ec346d4 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/RingRemoval.m @@ -0,0 +1,70 @@ +function proj = RingRemoval(proj, varargin) +% column-order based filtering +% Method: median or trimmean +% Input: +% proj: projection matrix +% varargin: filter type, filter kernel size +% Output: +% proj: filtered projection matrix +% Date: 2021-09-03 +% Author: Yi Du, yi.du@hotmail.com + +if(isempty(varargin)) + tag_filter = 'median'; + nkernel = 9; +elseif( length(varargin)==1 ) + if(ischar(varargin{1})) + if(contains(varargin{1},'median')) + tag_filter = 'median'; + elseif((contains(varargin{1},'trimmean'))) + tag_filter = 'trimmean'; + else + error('filter type error'); + end + else + error('filter type is missing'); + end + nkernel = 9; +elseif( length(varargin)== 2) + if(ischar(varargin{1})) + if(contains(varargin{1},'median')) + tag_filter = 'median'; + elseif((contains(varargin{1},'trimmean'))) + tag_filter = 'trimmean'; + else + error('filter type error'); + end + else + error('filter type is missing'); + end + + if(isnumeric(varargin{2})) + nkernel = varargin{2}; + end + + if(nkernel ~= round(varargin{2})) + error('filter kernel size error'); + end +else + error('input paramter error'); +end + + +% 3D projection matrix +if(contains(tag_filter,'median')) + idx = round( (nkernel + 1)/2); + for ii = 1:size(proj,3) + proj(:,:,ii) = ordfilt2(proj(:,:,ii), idx, ones(1,nkernel)); + end +elseif(contains(tag_filter,'trimmean')) + nb = round( (nkernel+1)/2) -1; + left = nb+1; + right = size(proj, 2) - left; + for ii = left: right + tmp_blk = proj(:, ii-nb : ii+nb, :); + proj(:,ii,:) = trimmean(tmp_blk, 40, 2, 'round'); + end +end + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m b/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m new file mode 100644 index 00000000..ad4b2350 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m @@ -0,0 +1,53 @@ +function kernel = SC_ASGkernel(sccalib, geo, dus, dvs) +%% Anti-sctter grid response function +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Input: +% geo: geometry structure +% dus: downsampled u vector +% dvs: downsampled v vector +% Output: +% kernel: anti-scatter grid +% +% Date: 2021-05-04 +% Author: Yi Du (yi.du@hotmail.com) + +%{ +%% Pixel coordinates +% center detector +offset = [0, 0]; + +% grid unit: mm +us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + offset(1); +vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2) + offset(2); +% unit mm - > cm +us = us/10; +vs = vs/10; + +%% Downsampling +% about 10 mm in axial direction (intra-plane) +dus = downsample(us, 26); +% about 4 mm in transaxial direction (inter-plane) +dvs = downsample(vs, 10); +%} + +%% Geometry +% unit mm +DSD = geo.DSD; + +%% Anti-Scatter Grid along X direction +% vs dimension +gamma = abs(rad2deg(atan(dvs'./DSD))); + +%% Transmission modelling +k = -0.15; +b = 1; +t_ratio = k.*gamma + b; + +t_ratio = k.*abs(dvs'/10) + b; + +%% Kernel: [nv, nu] +kernel = repmat(t_ratio, [1,length(dus)]); +efficiency = str2double(sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{1}.GridEfficiency.LamellaTransmission.Text); +kernel(kernel < efficiency) = efficiency; + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m b/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m new file mode 100644 index 00000000..329f83bb --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m @@ -0,0 +1,35 @@ +function cfactor = SC_AmplitudeFactor(blk, page, edgewt, sccalib) +% Amplitude Factor: ce_i with Edge Response Function (edgewt) +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Author: Yi Du (yi.du@hotmail.com) +% Date: 2021-05-24 + +%% group number +ngroup = length(sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel); + +% Amplitude factor groups +cfactor = []; + +for ii=1:ngroup + tmp = sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{ii}.ObjectScatterFit; + % Amplitude Factor + % unit: cm^2 - > mm^2 +% A = str2double(tmp.A.Text) / 100; + % unit: mm - > cm + A = str2double(tmp.A.Text) / 10; + % unitless + alpha = str2double(tmp.alpha.Text); + beta = str2double(tmp.beta.Text); + + % fill holes + term = (page + eps)./(blk + eps); + logterm = -log(term); + logterm(logterm<0) = NaN; + logterm = single(inpaint_nans(double(logterm), 2)); + + % amplitude factor wi edge response function as well + cfactor(:,:,ii) = A .*edgewt .* (term).^(alpha) .* ( logterm ).^(beta); +end + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_EdgeResponse.m b/MATLAB/Utilities/IO/VarianCBCT/SC_EdgeResponse.m new file mode 100644 index 00000000..5ba52360 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_EdgeResponse.m @@ -0,0 +1,23 @@ +function edgewt = SC_EdgeResponse(thickness) +% Ideal model: k = EdgeCoef0, (0.4), EdgeRange = EdgeRangeMM, (80) +% i.e., y = k*x + 0.78 +% Herein, we used an approximation method to model the edge response +% function +% Date: 2021-05-24 + + +% Emperical Thickness Threshold for Edge Segementation +edgewt = double(imbinarize(thickness, 50)); +tmpmask = edgewt; +h = fspecial('average', [25, 25]); +for ii = 1:5 + edgewt = imfilter(edgewt, h); +end +tmp = tmpmask.*edgewt; +tmp(tmp==0) = NaN; +tmp = rescale(tmp, 0.6, 1); +tmp(isnan(tmp)) = 0; +edgewt = tmp; + +end + diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_FormFunc.m b/MATLAB/Utilities/IO/VarianCBCT/SC_FormFunc.m new file mode 100644 index 00000000..ec23971c --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_FormFunc.m @@ -0,0 +1,32 @@ +function gform = SC_FormFunc(sccalib, dugd, dvgd) +%% Thickness-based Multiple Group Form Function kernels +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Date: 2021-05-18 +% Author: Yi Du + +%% group number +ngroup = length(sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel); + +% mm -> cm +unit_cvt = 1/10; +dugd = dugd *unit_cvt; +dvgd = dvgd * unit_cvt; + +%% Kernel form function groups +% Form function groups +gform = []; + +for ii=1:ngroup + tmp = sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{ii}.ObjectScatterFit; + % unit: cm^(-1) + sigma1 = str2double(tmp.sigma1.Text); + sigma2 = str2double(tmp.sigma2.Text); + % unitless + B = str2double(tmp.B.Text); + grid2 = dugd.^2 +dvgd.^2; + % Form Function + gform(:,:,ii) = exp( -0.5 * grid2 /(sigma1^2) ) + B * exp( -0.5 * grid2 /(sigma2^2) ); +end + + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_GroupMask.m b/MATLAB/Utilities/IO/VarianCBCT/SC_GroupMask.m new file mode 100644 index 00000000..7b5bcfca --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_GroupMask.m @@ -0,0 +1,30 @@ +function nmask = SC_GroupMask(thickness, ngroup, nbounds) +%% Generate n-group masks for estimated thickness (2D) +% +% SYNOPSIS: nmask = GroupMask(thickness, nbounds) +% +% INPUT nbounds: thickness based n-group bounds +% thickness: estimated thickness 2D matrix +% +% OUTPUT nmask(u, v, ngroup) +% +% REMARKS +% +% created with MATLAB ver.: 9.8.0.1538580 (R2020a) Update 6 on Microsoft Windows 10 Pro Version 10.0 (Build 18363) +% +% created by: Yi Du +% DATE: 18-May-2021 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% thickness based n-group +nmask = zeros(size(thickness, 1), size(thickness, 2), ngroup); + +%% n-group mask +for ii = 1:ngroup-1 + nmask(:, :, ii) = (thickness>nbounds(ii)).*(thicknessnbounds(ngroup); + +end \ No newline at end of file diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_SmoothThickness.m b/MATLAB/Utilities/IO/VarianCBCT/SC_SmoothThickness.m new file mode 100644 index 00000000..85a6af63 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_SmoothThickness.m @@ -0,0 +1,28 @@ +function thickness = SC_SmoothThickness(thickness, sccalib, step_du, step_dv) +%% Gaussian Filter to Smooth EstimatedThickness +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Input: +% thickness: Estimated object thickness, i.e., tau(x,y) in Reference +% Output: +% sccalib: Scatter Calibration Structure +% Date: 2021-05-05 +% Author: Yi Du (yi.du@hotmail.com) + +% unit: mm +sigma_u = str2double(sccalib.CalibrationResults.Globals.AsymPertSigmaMMu.Text); +sigma_v = str2double(sccalib.CalibrationResults.Globals.AsymPertSigmaMMv.Text); + +%% Gaussian filtering +% only one page +if(size(thickness, 3) == 1) + %% -------------- to check later the coordinate orientation + thickness = imgaussfilt(thickness, [sigma_v/step_dv sigma_u/step_du]); +% 3D thickness matrix +elseif(size(thickness, 3) >1) + for ii = 1: size(thickness, 3) + %% -------------- to check later the coordinate orientation + thickness(:,:,ii) = imgaussfilt(thickness(:,:,ii), [sigma_v/step_dv sigma_u/step_du]); + end +end + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_ThicknessEstimator.m b/MATLAB/Utilities/IO/VarianCBCT/SC_ThicknessEstimator.m new file mode 100644 index 00000000..413b51f2 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_ThicknessEstimator.m @@ -0,0 +1,32 @@ +function thickness = SC_ThicknessEstimator(blk, page, sccalib, step_du, step_dv) +%% Estimate Water-Equivalent Thickness (2D) +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Input: +% Blk(u,v): Total intensity, i.e., I_0 +% BlkAirNorm: +% Prm(u,v): Primary intensity, i.e., I_p +% AirNorm: Primary intensity airnorm chamber value +% sccalib: Scatter Calibration Structure +% Output: +% thickness(u,v): Estimated object thickness, i.e., tau(x,y) in Reference +% Date: 2021-05-05 +% Author: Yi Du (yi.du@hotmail.com) + +% mu H2O = 0.02 /mm +muH2O = str2double(sccalib.CalibrationResults.Globals.muH2O.Text); + +% unit mm +tmp = blk./page; +tmp(tmp<0)=0.0001; +thickness = log(tmp) /muH2O; + +% fill holes by interpolation +thickness(thickness<0) = NaN; + +thickness = single(inpaint_nans(double(thickness), 2)); + +%% Smooth the estimated thickness +% thickness(vv, uu, ntheta) +thickness = SC_SmoothThickness(thickness, sccalib, step_du, step_dv); + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/ScCalibFromXML.m b/MATLAB/Utilities/IO/VarianCBCT/ScCalibFromXML.m new file mode 100644 index 00000000..2f86e4a9 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/ScCalibFromXML.m @@ -0,0 +1,44 @@ +function ScCalibXML = ScCalibFromXML(datafolder) +%% Load Calibration.xml for Scatter Correction +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% Input: +% datafolder: Varian CBCT scan data folder +% Output: +% ScCalibXML: Scatter Calibration Structure +% Date: 2021-05-05 +% Author: Yi Du (yi.du@hotmail.com) + +foldername = [datafolder filesep 'Calibrations' filesep 'SC-*' filesep 'Factory']; +srcdir = dir([foldername filesep 'Calibration.xml']); +srcfilename = [srcdir.folder filesep 'Calibration.xml']; +tmpfilename = [srcdir.folder filesep 'tmp.xml']; +copyfile(srcfilename,tmpfilename); + +%% Delete redundent comment line that causes parsing error otherwise +fidsrc = fopen(srcfilename); +fidtmp = fopen(tmpfilename, 'wt'); +count = 0; +while(true) + tmp = fgetl(fidsrc); + % End-of-File + if(tmp == -1) + break; + end + + if(contains(tmp, '!--')) + continue; + end + + fprintf(fidtmp, '%s\n', tmp); + count = count+1; +end +fclose(fidsrc); +fclose(fidtmp); + +%% Export as struct +tmp = xml2struct(tmpfilename); +ScCalibXML = tmp.Calibration; +% delete temperary file +delete(tmpfilename); + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m new file mode 100644 index 00000000..e0afb49e --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m @@ -0,0 +1,153 @@ +function prim = ScatterCorrection(ScCalib, Blk, BlkAirNorm, proj, airnorm, geo) +% Scatter Correction Module +% Reference: Improved scatter correction using adaptive scatter kernel superposition +% No downsampling or upsampling is involved in current version +% Author: Yi Du (yi.du@hotmail.com) +% Date: 2021-05-24 + +%% Center Coordinates +%unit: mm +offset = geo.offDetector; +offset = [0, 0]; + +% grid unit: mm +us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) - offset(1); +vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2) - offset(2); + +%% Downsampling (not used in current version) +% downsampling rate +ds_rate = 12; +% unit: mm +% Downsampled to about 10 mm in axial direction in Ref +dus = decimate(us, ds_rate); +% Downsampled to about 4 mm in transaxial direction in Ref +dvs = decimate(vs, ds_rate); + +step_du = mean(diff(dus)); +step_dv = mean(diff(dvs)); + +%% X, Y grid +% unit: mm +[ugd,vgd] = meshgrid(us,vs); %detector +% downsampled grid +[dugd, dvgd] = meshgrid(dus,dvs); %detector + +%% Load Scatter Calibration +% ScCalib = ScCalibFromXML(datafolder); + +%% Blk scan +sBlk = sum(Blk, 3); +sAirNorm = sum(BlkAirNorm); + +%% n-thickness group number and boundaries +ngroup = length(ScCalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel); +nbounds = []; + +for ii=1:ngroup + % unit: mm + tmp = str2double(ScCalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{ii}.Thickness.Text); + nbounds = [nbounds, tmp]; +end + +%% k(y): anti-scatter grid kernel +ASG = SC_ASGkernel(ScCalib, geo, dus, dvs); + +%% Component Weights: gamma (gamma = 0 for SKS) +gamma = str2double(ScCalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{1}.ObjectScatterFit.gamma.Text); +% unit: cm-> mm +mm2cm = 1/10; + +% gamma = gamma * unit_cvt; + +%% iteration number +niter = 8; +% relaxation factor in iteration +lambda = 0.005; + +%% Primary signal matrix +prim = zeros(size(proj)); + +zeropage = zeros(size(dugd)); % preallocation + +for ii = 1: size(proj, 3) + if(~mod(ii,50)) + display([num2str(ii) '/' num2str(size(proj, 3))]); + end + + % blk: I_0 unattenuated blk signal + CF = sAirNorm/airnorm(ii); + blk = interp2(ugd, vgd, sBlk/CF, dugd, dvgd, 'linear', 0); + + % page: I_p primary signal + % note: detector point spread deconvolution should be done first + page = interp2(ugd, vgd, proj(:,:,ii), dugd, dvgd, 'linear', 0); + + % Is initialization + Is = zeropage; + comp1 = zeropage; + comp2 = zeropage; + + %% Iterative Correction + for jj = 1: niter + % Is previous scatter map + Is_prv = Is; + + % estimate thickness: mm + thickness = SC_ThicknessEstimator(blk, page, ScCalib, step_du, step_dv); + % smooth thickness + thickness = SC_SmoothThickness(thickness, ScCalib, step_du, step_dv); + + % Ri(x,y): group-based masks + nmask = SC_GroupMask(thickness, ngroup, nbounds); + + % gi(x,y): group-based form function + gform = SC_FormFunc(ScCalib, dugd, dvgd); + + % edge response function: about 7.5 cm inward + edgewt = SC_EdgeResponse(thickness); + + % cei(x,y): group-based amplitude factors + cfactor = SC_AmplitudeFactor(blk, page, edgewt, ScCalib); + + % mm -> cm + thickness = thickness * mm2cm; + %% n-group summation + term1 = repmat(page,[1,1, ngroup]).*nmask.*cfactor; + term2 = fft2(gform.* repmat(ASG, [1,1, ngroup])); + + tmp1 = sum(fft2(term1).*term2,3); + tmp2 = sum(fft2(repmat(thickness, [1,1, ngroup]).* term1).*term2,3); + + comp1 = real(ifft2(tmp1)); + comp2 = real(ifft2(tmp2)); + %% fASKS scatter correction + Is = (1 - gamma .*thickness).*comp1 + gamma.*comp2; + %sum(sum(Is_prv - Is)) + page = page + lambda * (Is_prv - Is); + page(page<0) = eps; + end + + %% Upsampling and cutoff for over-correction + % measured intensity + scmap = interp2(dugd, dvgd, Is, ugd, vgd, 'spline'); + % extrolation errors + scmap(isnan(scmap)) = eps; + % SF = Is; + scmap(scmap<0) = eps; + + % scatter fraction + SF = scmap./proj(:,:,ii); + % in case of zeros in proj + SF(SF==Inf) = NaN; + SF(SF>1000) = NaN; + SF = single(inpaint_nans(double(SF),0)); + % median filtering + SF = medfilt2(SF, [3 3]); + + % Scatter fraction cutoff + SF = min(SF, 0.95); + % primary + prim(:,:,ii) = proj(:,:,ii).*(1 - SF); +end + +end diff --git a/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m b/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m index f85d3fe5..05aef7c1 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m +++ b/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m @@ -1,55 +1,105 @@ -function [proj,geo, angles ] = VarianDataLoader(datafolder, varargin) -% VarianDataLoader Loads Varian CBCT projection, geomtry and angles data -% +function [proj_lg, geo, angles] = VarianDataLoader(datafolder, varargin) +% VarianDataLoader: Loads Varian CBCT projection, geomtry and angles data % Optional parameter: Motion lag correction. Default True. % % Load all dataset that are needed for reconstruction -% Date: 2020-04-16 +% Tested on TrueBeam 2.0 and 2.7 +% Date: 2021-04-02 % Author: Yi Du (yi.du@hotmail.com) -% datafolder = 'E:\BigData\Edge\CBCT_Export\2020-01-09_144244'; +% datafolder = '~/your_data_path/varian/2020-01-01_123456/'; + +%% GPU initialization +reset(gpuDevice(1)); + +%% Input Parser +% ACDC: acceleration & deceleration correction (default: true) +% DPS: Detector Point Spread correction (default: true) +% SC: Scatter Correction (default: true) +% BH: Beam Hardening correction (default: false, due to Bowtie, BH correction is not required at all.) +[tag_ACDC, tag_DPS, tag_SC, tag_BH] = parse_inputs(varargin{:}); %% Load geometry [geo, ScanXML] = GeometryFromXML(datafolder); -%% Motion lag correcion +%% Load Scatter Correction Calibration Parameter Structure +if(tag_DPS || tag_SC) + ScCalib = ScCalibFromXML(datafolder); +end + +%% Remove over-sampled projections due to acceleration and deceleration thd = 0; -if(~isempty(varargin)&&(varargin{1})) - thd = str2double(ScanXML.Acquisitions.Velocity.Text)... +if(tag_ACDC) + angular_interval = str2double(ScanXML.Acquisitions.Velocity.Text)... ./str2double(ScanXML.Acquisitions.FrameRate.Text); - thd = thd *0.95; + thd = angular_interval *0.9; +end + +%% Load proj and angles +disp('Loading Proj: ') +[proj, angles, airnorm] = ProjLoader(datafolder, thd); +% Detector point scatter correction +if(tag_DPS) + disp('Proj DPS: ') + proj = DetectorPointScatterCorrection(proj, geo, ScCalib); end -%% Load proj and angle -[proj, angles, blk] = BatchReadXim(datafolder, thd); +%% Load blank scan +disp('Loading Blk: ') +[Blk, Sec, BlkAirNorm] = BlkLoader(datafolder); +% Detector point scatter correction +if(tag_DPS) + disp('Blk DPS: ') + Blk = DetectorPointScatterCorrection(Blk, geo, ScCalib); +end +%% Scatter Correction +if(tag_SC) + disp('Scatter correction onging: ') + proj = ScatterCorrection(ScCalib, Blk, BlkAirNorm, proj, airnorm, geo); + disp('Scatter correction is completed.') +end -%% Logarithmic calculation -proj = log(repmat(blk, [1 1 size(proj,3)])./proj); +%% Airnorm and Logarithmic Normalization +proj_lg = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm); +disp('Log Normalization is completed.') +% remove anomolies +proj_lg = EnforcePositive(proj_lg); -%% Mediate filtering along colume-orth -for ii = 1:size(proj,3) - proj(:,:,ii) = ordfilt2(proj(:,:,ii), 5, ones(1,9)); +%% Beam Hardening correction is applied (kind of slow) +if(tag_BH) + [proj_lg, ~] = BHCorrection(datafolder, geo, ScanXML, proj_lg); end -% in case of abnormlies -proj(isnan(proj)) = 0; -proj(isinf(proj)) = 0; +%% Remove anomalies +proj_lg = EnforcePositive(proj_lg); -% all negative to zeros -proj(proj<0) = 0; +%% mediant filtering along colume-order +proj_lg = RingRemoval(proj_lg); -% double to single -proj = single(proj); +%% double to single +proj_lg = single(proj_lg); +angles = deg2rad(angles); -% degree to rad -angles = angles/180*pi; +% +disp('Data processing is complete! Ready for reconstruction: ') -%% Gantry Rotation correction -% Clockwise -if(angles(end) - angles(1)>0) - proj = flip(proj, 3); -% Counter-clockwise -> Clockwise -else - angles = flip(angles); end + +function [tag_ACDC, tag_DPS, tag_SC, tag_BH] = parse_inputs(varargin) +% create input parser +p = inputParser; +% add optional parameters +addParameter(p,'acdc', true); +addParameter(p,'dps', true); +addParameter(p,'sc', true); +% BH performance is very lame for unclear reason +addParameter(p,'bh', false); + +%execute +parse(p,varargin{:}); +%extract +tag_ACDC=p.Results.acdc; +tag_DPS=p.Results.dps; +tag_SC=p.Results.sc; +tag_BH=p.Results.bh; end diff --git a/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp b/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp index 138a9738..88bddf96 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp +++ b/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp @@ -3,7 +3,10 @@ #include #include -//typedef struct XimPara +// Purpose: To fast read .xim files +// Method: based on ReadXim.m by Fredrik Nordström 2015 +// Date: 2017.07 +// Author: Yi Du, yi.du@hotmail.com #ifndef STR_XIM #define STR_XIM @@ -22,7 +25,3 @@ typedef struct XimPara }XimPara; #endif -//#ifndef cReadXim_FUN -//#define cReadXim_FUN -// int cReadXim(char *XimFullFile, XimPara *XimStr, int *XimImg); -//#endif diff --git a/MATLAB/Utilities/IO/VarianCBCT/interp_weight.m b/MATLAB/Utilities/IO/VarianCBCT/interp_weight.m new file mode 100644 index 00000000..13ccd4be --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/interp_weight.m @@ -0,0 +1,35 @@ +function [ii, weights] = interp_weight(x, X) + % LOCATE Locate first node on grid below a given value. + % + % [ii, weights] = locate(x, X) returns the first node in X that is below + % each element in x and the relative proximities to the two closest nodes. + % + % X must be a monotonically increasing vector. x is a matrix (of any + % order). + + % Preallocate + ii = ones(size(x)); % Indices of first node below (or 1 if no nodes below) + weights = zeros([2, size(x)]); % Relative proximity of the two closest nodes + + % Find indices + for iX = 1:length(X) - 1 + ii(X(iX) <= x) = iX; + end + + % Find weights + below = x <= X(1); + weights(1, below) = 1; % All mass on the first node + weights(2, below) = 0; + + above = x >= X(end); + weights(1, above) = 0; + weights(2, above) = 1; % All mass on the last node + + interior = ~below & ~above; + xInterior = x(interior)'; + iiInterior = ii(interior); + XBelow = X(iiInterior)'; + XAbove = X(iiInterior + 1)'; + weights(:, interior) = ... + [XAbove - xInterior; xInterior - XBelow] ./ (XAbove - XBelow); +end \ No newline at end of file diff --git a/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp b/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp index 81d1e9e2..35fca59e 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp +++ b/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp @@ -14,6 +14,13 @@ #include "matrix.h" #include "XimPara.hpp" +#define GET_BIT(x,bit) ((x & (1 << bit)) >>bit) + +// Purpose: To fast read .xim files +// Method: based on ReadXim.m by Fredrik Nordström 2015 +// Date: 2017.07 +// Author: Yi Du, yi.du@hotmail.com + int cReadXim(char *XimFullFile, XimPara *XimStr, int *XimImg); @@ -131,11 +138,22 @@ int cReadXim(char *XimFullFile, for (int ii = 0; ii < LookUpTableSize; ii++) { // Load in the 8-bit date - //unsigned __int8 tmp = 0; - //fread(&tmp, sizeof(unsigned __int8), 1, fid); - int tmp = 0; - fread(&tmp, sizeof(uint8_T), 1, fid); + // Updated: 2021-11-05, Yi Du + uint8_T tmp =0; + fread(&tmp, 1, 1, fid); + int Bit2[4] = { 0 }; + Bit2[0] = GET_BIT(tmp,0) + GET_BIT(tmp,1) *2; + Bit2[1] = GET_BIT(tmp,2) + GET_BIT(tmp,3) *2; + Bit2[2] = GET_BIT(tmp,4) + GET_BIT(tmp,5) *2; + Bit2[3] = GET_BIT(tmp,6) + GET_BIT(tmp,7) *2; + // extract the lookup_table data + for (int jj = 0; jj < 4; jj++) + { + LookUpTable[ii * 4 + jj] = Bit2[jj]; + } + + /** Old Code with bug int Bit2[4] = { 0 }; // extract the lookup_table data @@ -147,6 +165,7 @@ int cReadXim(char *XimFullFile, //printf("Index = %d, LookUpTable = %d\n", ii * 4 + jj / 2, LookUpTable[ii * 4 + jj / 2]); } + **/ } // Skip compressed_pixel_buffer_size: passed @@ -159,9 +178,10 @@ int cReadXim(char *XimFullFile, int delta = 0; int LUT_Pos = 0; - // You must be very careful with all data types!!! + // Be very careful with all data types!!! int8_T tmp8 = 0; int16_T tmp16 = 0; + int32_T tmp32 = 0; for (int ImgTag = XimStr->ImgWidth + 1; ImgTag < (XimStr->ImgHeight) * (XimStr->ImgWidth); @@ -179,7 +199,8 @@ int cReadXim(char *XimFullFile, } else { - fread(&delta, sizeof(int32_T), 1, fid); + fread(&tmp32, sizeof(int32_T), 1, fid); + delta = int(tmp32); } XimImg[ImgTag] = delta + XimImg[ImgTag - 1] diff --git a/MATLAB/Utilities/IO/VarianCBCT/polar2cart/ImToPolar.m b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/ImToPolar.m new file mode 100644 index 00000000..6c676822 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/ImToPolar.m @@ -0,0 +1,61 @@ +function imP = ImToPolar (imR, rMin, rMax, M, N) +% IMTOPOLAR converts rectangular image to polar form. The output image is +% an MxN image with M points along the r axis and N points along the theta +% axis. The origin of the image is assumed to be at the center of the given +% image. The image is assumed to be grayscale. +% Bilinear interpolation is used to interpolate between points not exactly +% in the image. +% +% rMin and rMax should be between 0 and 1 and rMin < rMax. r = 0 is the +% center of the image and r = 1 is half the width or height of the image. +% +% V0.1 7 Dec 2007 (Created), Prakash Manandhar pmanandhar@umassd.edu + +[Mr Nr] = size(imR); % size of rectangular image +Om = (Mr+1)/2; % co-ordinates of the center of the image +On = (Nr+1)/2; +sx = (Mr-1)/2; % scale factors +sy = (Nr-1)/2; + +imP = zeros(M, N); + +delR = (rMax - rMin)/(M-1); +delT = 2*pi/N; + +% loop in radius and +for ri = 1:M +for ti = 1:N + r = rMin + (ri - 1)*delR; + t = (ti - 1)*delT; + x = r*cos(t); + y = r*sin(t); + xR = x*sx + Om; + yR = y*sy + On; + imP (ri, ti) = interpolate (imR, xR, yR); +end +end + +function v = interpolate (imR, xR, yR) + xf = floor(xR); + xc = ceil(xR); + yf = floor(yR); + yc = ceil(yR); + if xf == xc & yc == yf + v = imR (xc, yc); + elseif xf == xc + v = imR (xf, yf) + (yR - yf)*(imR (xf, yc) - imR (xf, yf)); + elseif yf == yc + v = imR (xf, yf) + (xR - xf)*(imR (xc, yf) - imR (xf, yf)); + else + A = [ xf yf xf*yf 1 + xf yc xf*yc 1 + xc yf xc*yf 1 + xc yc xc*yc 1 ]; + r = [ imR(xf, yf) + imR(xf, yc) + imR(xc, yf) + imR(xc, yc) ]; + a = A\double(r); + w = [xR yR xR*yR 1]; + v = w*a; + end diff --git a/MATLAB/Utilities/IO/VarianCBCT/polar2cart/PolarToIm.m b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/PolarToIm.m new file mode 100644 index 00000000..2c4051ea --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/PolarToIm.m @@ -0,0 +1,71 @@ +function imR = PolarToIm (imP, rMin, rMax, Mr, Nr) +% POLARTOIM converts polar image to rectangular image. +% +% V0.1 16 Dec, 2007 (Created) Prakash Manandhar, pmanandhar@umassd.edu +% +% This is the inverse of ImToPolar. imP is the polar image with M rows and +% N columns of data (double data between 0 and 1). M is the number of +% samples along the radius from rMin to rMax (which are between 0 and 1 and +% rMax > rMin). Mr and Nr are the number of pixels in the rectangular +% domain. The center of the image is assumed to be the origin for the polar +% co-ordinates, and half the width of the image corresponds to r = 1. +% Bilinear interpolation is performed for points not in the imP image and +% points not between rMin and rMax are rendered as zero. The output is a Mr +% x Nr grayscale image (with double values between 0.0 and 1.0). + + +imR = zeros(Mr, Nr); +Om = (Mr+1)/2; % co-ordinates of the center of the image +On = (Nr+1)/2; +sx = (Mr-1)/2; % scale factors +sy = (Nr-1)/2; + +[M N] = size(imP); + +delR = (rMax - rMin)/(M-1); +delT = 2*pi/N; + +for xi = 1:Mr +for yi = 1:Nr + x = (xi - Om)/sx; + y = (yi - On)/sx; + r = sqrt(x*x + y*y); + if r >= rMin & r <= rMax + t = atan2(y, x); + if t < 0 + t = t + 2*pi; + end + imR (xi, yi) = interpolate (imP, r, t, rMin, rMax, M, N, delR, delT); + end +end +end + +function v = interpolate (imP, r, t, rMin, rMax, M, N, delR, delT) + ri = 1 + (r - rMin)/delR; + ti = 1 + t/delT; + rf = floor(ri); + rc = ceil(ri); + tf = floor(ti); + tc = ceil(ti); + if tc > N + tc = tf; + end + if rf == rc & tc == tf + v = imP (rc, tc); + elseif rf == rc + v = imP (rf, tf) + (ti - tf)*(imP (rf, tc) - imP (rf, tf)); + elseif tf == tc + v = imP (rf, tf) + (ri - rf)*(imP (rc, tf) - imP (rf, tf)); + else + A = [ rf tf rf*tf 1 + rf tc rf*tc 1 + rc tf rc*tf 1 + rc tc rc*tc 1 ]; + z = [ imP(rf, tf) + imP(rf, tc) + imP(rc, tf) + imP(rc, tc) ]; + a = A\double(z); + w = [ri ti ri*ti 1]; + v = w*a; + end diff --git a/MATLAB/Utilities/IO/VarianCBCT/polar2cart/license.txt b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/license.txt new file mode 100644 index 00000000..6607e051 --- /dev/null +++ b/MATLAB/Utilities/IO/VarianCBCT/polar2cart/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2009, Prakash Manandhar +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/MATLAB/Utilities/im2DDenoise.m b/MATLAB/Utilities/im2DDenoise.m new file mode 100644 index 00000000..31780ccb --- /dev/null +++ b/MATLAB/Utilities/im2DDenoise.m @@ -0,0 +1,24 @@ +function img = im2DDenoise(img, type, kernelsize) +%img2DDenoise removes noise of image with 2D kernel-based methods +% Currentyl only median and average filter is supported. +% type: 'median' or 'average' +% kernelsize: filter kernel size +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% Date: 2021-09-21 +% Author: Yi Du, +% Contact: yi.du@hotmail.com + + +if strcmp(type,'median') + SliceNum = size(img,3); + for kk =1:SliceNum + img(:,:,kk) = medfilt2(img(:,:,kk), [kernelsize, kernelsize]); + end +elseif strcmp(type,'average') + h = fspecial('average', kernelsize); + img = imfilter(img, h, 'same'); +end + +end + From 2798d93cecd291956269ac6d50f2c879c86847b7 Mon Sep 17 00:00:00 2001 From: Yi Du Date: Mon, 23 May 2022 17:23:49 +0800 Subject: [PATCH 2/4] Add one more output: KVNormChamber --- MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp | 1 + MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp | 21 ++++++++++++++----- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp b/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp index 88bddf96..670c2d3e 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp +++ b/MATLAB/Utilities/IO/VarianCBCT/XimPara.hpp @@ -22,6 +22,7 @@ typedef struct XimPara int Compression_Indicator; // Data number in Rec Image Matrix double GantryRtn; // Gantry rotation angle + int KVNormChamber; // KV norm chamber reading, date: 2022-05-23 }XimPara; #endif diff --git a/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp b/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp index 35fca59e..453c4278 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp +++ b/MATLAB/Utilities/IO/VarianCBCT/mexReadXim.cpp @@ -84,6 +84,9 @@ void mexFunction( // KVSourceRtn = GantryRtn + 90 deg; double KVSourceRtn = para->GantryRtn + 90; plhs[1] = mxCreateDoubleScalar(KVSourceRtn); + + double NormChamberReading = para->KVNormChamber * 1.0; + plhs[2] = mxCreateDoubleScalar(NormChamberReading); } @@ -271,7 +274,7 @@ int cReadXim(char *XimFullFile, { int pName_len = 0; // Only load the property name rather than the content - char pName[64] = { 0 }; + char pName[128] = { 0 }; int pType = 0; for (int ii = 0; ii < nProperties; ii++) { @@ -281,13 +284,21 @@ int cReadXim(char *XimFullFile, fread(pName, sizeof(char)* pName_len, 1, fid); // load property data type fread(&pType, sizeof(int), 1, fid); + + //printf("%s\n", pName); // extract the Gantry Rotation Angle if (!strcmp(pName, "GantryRtn")) { fread(&(XimStr->GantryRtn), sizeof(double), 1, fid); - break; +// continue; } + else if(!strcmp(pName, "KVNormChamber")) + { + //printf("KVNormChamber"); + fread(&(XimStr->KVNormChamber), sizeof(int), 1, fid); + break; + } else { switch (pType) @@ -313,14 +324,14 @@ int cReadXim(char *XimFullFile, { int skiplen = 0; fread(&skiplen, sizeof(int), 1, fid); - fseek(fid, sizeof(double) * skiplen, SEEK_CUR); + fseek(fid, sizeof(double) * skiplen /8, SEEK_CUR); break; } case 5: { int skiplen = 0; fread(&skiplen, sizeof(int), 1, fid); - fseek(fid, sizeof(int) * skiplen, SEEK_CUR); + fseek(fid, sizeof(int) * skiplen /4, SEEK_CUR); break; } break; @@ -328,7 +339,7 @@ int cReadXim(char *XimFullFile, } // reset all the temporary variables pName_len = 0; - memset(pName, 0, 64*sizeof(char)); + memset(pName, 0, 128*sizeof(char)); pType = 0; } From 2a2b428fbb7f2c994d4b885edc8c1805fb1b24c9 Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Mon, 23 May 2022 10:39:38 +0100 Subject: [PATCH 3/4] Add pixel2HU and warning --- Common/data/demo_Pixel2HU.mat | Bin 0 -> 11496 bytes MATLAB/Utilities/IO/VarianCBCT/HUMapping.m | 7 ++++++- 2 files changed, 6 insertions(+), 1 deletion(-) create mode 100644 Common/data/demo_Pixel2HU.mat diff --git a/Common/data/demo_Pixel2HU.mat b/Common/data/demo_Pixel2HU.mat new file mode 100644 index 0000000000000000000000000000000000000000..5fe2081c77bd4438e6a28b348ac6cbf39442b4a4 GIT binary patch literal 11496 zcmeHMYfMx}6rQ^)vMLH%>#H_iA8jZB3d=)FtOpe9f}kspf&``4HTis;vu@~@09w-V3YGUO1Prs2jN}LBvR^GS{x$P6CL*JCxop!Lvbj9j zA)B3(iDyxo#v2}x3n-wQ|Auh?W&Wrv3Xj#B|Ds`!Hva=IQ&szqBf$N)m95jl`B+-u zel=zsl~1ttAQYvP1T+-s@CP5#_`EJK*(9fvG(3U-2I4>N-%gjK+~guHlrrbdPZp5I z8-GIsX+76BGI{ABTF>LFsO_$2V6uhzJ#Fnj-}gGleDUio4(Nt-$GQI?sQyxd0FMLw z91CHvxY_>WC`qiq_o;z5JR%ox>>y)>umJ@Y5d?u&@MDy7Ek5TIecXq}ghXZL05*Y@I1)3T-u4X`(FarNT% zW>5xBeDjdPk(Kh6jfqrVNaZziDCc^I@@c!Ll=2f5l%K9vc+dS1g$O<9`Ce!*Ut*KV zKH&WVF-|&Uh%+W-8R&h0xGvu}&--48y&kSii>TiRB9V&7KsL{Pi0=<$(0!C}6SX=* zH~eomZ|6i-${p#g@>{?^O(ESVzs=rC9#=mZ?XB`#piu;&mJxoNAG)$H{r;o8b4N7h z5_mSDJ;(ls8@6_QfkKNtaHQ`0jd{~*U{7G!zPDG_LbTLY{a$)86dx51WnQlXZQZw_ z+jiAMMB9<`FAr>l)<37#XJ6|de&3h&a_l@9Iq_Nv%s-i-H^*;-^6$n}+AcOjp!4`^ z`;XT^L-}$`_NiJ}H{N{XLS_lf73aPfa%(%J{r1!7XR;gN&D{@M;-ChbU__(nbCMW-et zS9$tHJmeQjgDe#m7_3E>1?7^>O5ePsjnosKQYsb6bL29YO?eKZEQ{SCEl0>G+bXO$ z<41X#OZWj<2xYcH-()?O F{sA*T%+mk> literal 0 HcmV?d00001 diff --git a/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m b/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m index f3eaee7b..08b74b33 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m +++ b/MATLAB/Utilities/IO/VarianCBCT/HUMapping.m @@ -10,7 +10,12 @@ % Author: Yi Du, yi.du@hotmail.com HU3D = zeros(numel(img3D), 1); -load('D:\MATLABWorkplace\Ongoing_Projects\Fronc_VarianCBCT\demo_Pixel2HU.mat', 'Pixel2HU'); +if ~ isfile('../../../../Common/data/pixel2HU.mat ') + warning('Using DEMO pixel2HU transform. We recommend creating your own and storing it in TIGRE/Common/data/pixel2HU.mat') + load('../../../../Common/data/demo_Pixel2HU.mat', 'Pixel2HU'); +else + load('../../../../Common/data/Pixel2HU.mat', 'Pixel2HU'); +end tmp = getfield(Pixel2HU, Protocol); From 42cc0ee2246212bc4c0739977a9b2b6fed73690d Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Mon, 23 May 2022 11:09:01 +0100 Subject: [PATCH 4/4] Add gpu selection (and defaults), remove comments --- MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m | 4 ++-- .../IO/VarianCBCT/BH_ObjectCalibLUT.m | 15 ------------ .../IO/VarianCBCT/BH_ObjectRemapping.m | 22 +++--------------- .../IO/VarianCBCT/BH_SpectrumBowtieLUT.m | 23 +------------------ .../BH_SpectrumBowtieLUT_deprecated.m | 7 +++--- .../DetectorPointScatterCorrection.m | 7 ++++-- MATLAB/Utilities/IO/VarianCBCT/LogNormal.m | 6 +++-- MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m | 20 ---------------- .../IO/VarianCBCT/SC_AmplitudeFactor.m | 2 -- .../IO/VarianCBCT/ScatterCorrection.m | 3 --- .../IO/VarianCBCT/VarianDataLoader.m | 20 +++++++++------- 11 files changed, 31 insertions(+), 98 deletions(-) diff --git a/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m index a9155039..4bc43a2b 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m +++ b/MATLAB/Utilities/IO/VarianCBCT/BHCorrection.m @@ -1,4 +1,4 @@ -function [proj_lg, BHCalib] = BHCorrection(datafolder, geo, ScanXML, proj_lg) +function [proj_lg, BHCalib] = BHCorrection(datafolder, geo, ScanXML, proj_lg,gpuids) % Entry Function For BH Correction % Detailed explanation goes here @@ -28,7 +28,7 @@ % BH correction via reference object (water) BHCalib = BH_RemappingFunc(BHCalib); -proj_lg = BH_ObjectRemapping(BHCalib, proj_lg); +proj_lg = BH_ObjectRemapping(BHCalib, proj_lg, gpuids); disp('BH correction is done.') diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m index bcd8d83e..c15bf013 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectCalibLUT.m @@ -15,21 +15,6 @@ % object thickness look-up table: [bowtie.sl, object_sl] calibLUT = zeros(length(BHCalib.bowtie.sl), length(object_sl)); -%% Photo flux based calibration -%{ -% bowtie thickness - > spectrum -for ii = 1: length(BHCalib.bowtie.sl) - % object thickness -> [min_bowtie_thickness: max_bowtie_thickness] - spec = specLUT(ii,:); - for jj = 1:length(object_sl) - % non-ideal attenuated signal - tmp = spec .* exp(-object_sl(jj).* miu); - % nonlinear BH projection - calibLUT(ii,jj) = -log(sum(tmp)/sum(spec)); - end -end -%} - %% Energy Flux based calibration % scintillator energy absorption coefficient scintillator_ac = BHCalib.scintillator.ac(1: length(miu)); diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m index 3fc5c691..13596e99 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_ObjectRemapping.m @@ -1,8 +1,9 @@ -function proj_BH = BH_ObjectRemapping(BHCalib, projlg) +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); @@ -13,19 +14,6 @@ % Object Sampling Length object_sl = BHCalib.object.sl; -%% Remap non-ideal BH signal to linear object thickness -%{ -proj_BH = zeros(size(projlg)); - -% for 1D bowtie only -for jj = 1: nCol - tl_v = ulgd(1, jj); - % model the LUT for specific bowtie thickness - proj_signal_LUT = interp1(BHCalib.bowtie.sl, calibLUT, tl_v, 'spline'); - proj_BH(:,jj,:) = interp1(proj_signal_LUT, object_sl, projlg(:,jj,:), 'spline'); -end -%} - %% CPU Version: very very slow %{ proj_BH = zeros(size(projlg)); @@ -66,12 +54,9 @@ %% GPU version: acceptable % GPU Reset -g = gpuDevice(1); +g = gpuDevice(gpuids.devices(0)+1); reset(g); -% gpuArray preparation -% proj_BH = single(zeros(size(projlg))); -% gproj_BH = zeros(size(proj_BH), 'gpuArray'); gprojlg = gpuArray(single(projlg)); gsl = gpuArray(BHCalib.bowtie.sl); @@ -111,7 +96,6 @@ proj_BH(proj_BH<0) = eps; % GPU reset -g = gpuDevice(1); reset(g); end diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m index 62c0f5e3..32e4093d 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT.m @@ -33,21 +33,7 @@ %% Fitting to Model Bowtie Profile: [xData, yData] = prepareCurveData(BHCalib.bowtie.uu, BHCalib.bowtie.thickness); -%{ -%% Fitting Model: first fitting(linear) -% Set up fittype and options. -ft = fittype( 'poly1' ); -% Fit model to data. -[~, gof] = fit( xData, yData, ft ); -if( gof.rsquare < 0.995) - disp('It seems that BH correction is not required at all'); - typein = input('Continue BH correction? Y or N (recommended): ', 's'); - if(contains(typein, 'n', 'IgnoreCase', true)) - BHCalib = NaN; - return; - end -end -%} + %% Fitting Model: Smoothing Spline % Set up fittype and options. ft = fittype( 'smoothingspline' ); @@ -78,13 +64,6 @@ ulgd = Ax(bt_img3D, geo_bt, 0, 'interpolated'); clearvars bt_img bt_img3D xxgd yygd -%% -%{ -% tranverse length vector: mm -ul = interp1(BHCalib.bowtie.uu, BHCalib.bowtie.thickness, ub); -% tranverse length grid: mm -ulgd = repmat(ul, [length(vs), 1]); -%} %% bowtie material attenuation look-up table [Min, Max] = bounds(ulgd, 'all'); diff --git a/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT_deprecated.m b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT_deprecated.m index ccaaefe3..71303515 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT_deprecated.m +++ b/MATLAB/Utilities/IO/VarianCBCT/BH_SpectrumBowtieLUT_deprecated.m @@ -24,12 +24,13 @@ %% bowtie material attenuation look-up table [Min, Max] = bounds(ulgd(:)); % sampling length -sl = linspace(Min, Max, 100); +n_samples=100; +sl = linspace(Min, Max, n_samples); % LUT: [sl, miu_bowtie] -atten_tab = zeros(length(sl), length(BHCalib.bowtie.ac)); +atten_tab = zeros(n_samples, length(BHCalib.bowtie.ac)); % I/I_0 = exp(- length * miu): dependent on incident energy -for ii = 1:length(sl) +for ii = 1:n_samples atten_tab(ii,:) = exp(-sl(ii).* BHCalib.bowtie.ac'); end diff --git a/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m index d1e61704..fc3f0183 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m +++ b/MATLAB/Utilities/IO/VarianCBCT/DetectorPointScatterCorrection.m @@ -1,9 +1,10 @@ -function proj = DetectorPointScatterCorrection(proj, geo, ScCalib) +function proj = DetectorPointScatterCorrection(proj, geo, ScCalib, gpuids) %% Detector Point Scatter Correction % Reference: Improved scatter correction using adaptive scatter kernel superposition % Date: 2021-03-26 % Author: Yi Du (yi.du@hotmail.com) + %% Empirical values from reference paper % unit: cm-2 % a0 = 3.43; (in refrence paper, but incorrect) @@ -69,6 +70,7 @@ %% GPU based % reset(gpuDevice(1)); +reset(gpuDevice(gpuids.devices(0)+1)); gproj = gpuArray(single(proj)); %% 2D Convolution with downsampling and upsampling @@ -98,7 +100,8 @@ proj = single(gather(gproj)); % Reset GPU -reset(gpuDevice(1)); +reset(gpuDevice(gpuids.devices(0)+1)); + %% Cutoff for over-correction proj(proj<0) = NaN; diff --git a/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m b/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m index 33203747..7f82abe0 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m +++ b/MATLAB/Utilities/IO/VarianCBCT/LogNormal.m @@ -1,4 +1,4 @@ -function proj = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm) +function proj = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm,gpuids) % Log Normalization: % Calculate Logrithmic Projections with AirNorm % Tested on TrueBeam 2.5 and 2.7 @@ -39,6 +39,7 @@ toc %} %% TB version: GPU +gpuDevice(gpuids.devices(0)+1); gproj = gpuArray(single(proj)); gBlk = gpuArray(Blk); % Version = 2.5 @@ -71,5 +72,6 @@ proj = gather(gproj); % GPU clear -reset(gpuDevice(1)); +reset(gpuDevice(gpuids.devices(0)+1)); + end diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m b/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m index ad4b2350..a4befdc8 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_ASGkernel.m @@ -11,25 +11,6 @@ % Date: 2021-05-04 % Author: Yi Du (yi.du@hotmail.com) -%{ -%% Pixel coordinates -% center detector -offset = [0, 0]; - -% grid unit: mm -us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1) + offset(1); -vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2) + offset(2); -% unit mm - > cm -us = us/10; -vs = vs/10; - -%% Downsampling -% about 10 mm in axial direction (intra-plane) -dus = downsample(us, 26); -% about 4 mm in transaxial direction (inter-plane) -dvs = downsample(vs, 10); -%} - %% Geometry % unit mm DSD = geo.DSD; @@ -41,7 +22,6 @@ %% Transmission modelling k = -0.15; b = 1; -t_ratio = k.*gamma + b; t_ratio = k.*abs(dvs'/10) + b; diff --git a/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m b/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m index 329f83bb..091cef8a 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m +++ b/MATLAB/Utilities/IO/VarianCBCT/SC_AmplitudeFactor.m @@ -13,8 +13,6 @@ for ii=1:ngroup tmp = sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{ii}.ObjectScatterFit; % Amplitude Factor - % unit: cm^2 - > mm^2 -% A = str2double(tmp.A.Text) / 100; % unit: mm - > cm A = str2double(tmp.A.Text) / 10; % unitless diff --git a/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m b/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m index e0afb49e..2f04eb83 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m +++ b/MATLAB/Utilities/IO/VarianCBCT/ScatterCorrection.m @@ -32,9 +32,6 @@ % downsampled grid [dugd, dvgd] = meshgrid(dus,dvs); %detector -%% Load Scatter Calibration -% ScCalib = ScCalibFromXML(datafolder); - %% Blk scan sBlk = sum(Blk, 3); sAirNorm = sum(BlkAirNorm); diff --git a/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m b/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m index 05aef7c1..cd514c25 100644 --- a/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m +++ b/MATLAB/Utilities/IO/VarianCBCT/VarianDataLoader.m @@ -1,22 +1,23 @@ function [proj_lg, geo, angles] = VarianDataLoader(datafolder, varargin) % VarianDataLoader: Loads Varian CBCT projection, geomtry and angles data % Optional parameter: Motion lag correction. Default True. -% +% gpuids: Usable GPUs. Default, the first one. % Load all dataset that are needed for reconstruction % Tested on TrueBeam 2.0 and 2.7 % Date: 2021-04-02 % Author: Yi Du (yi.du@hotmail.com) % datafolder = '~/your_data_path/varian/2020-01-01_123456/'; -%% GPU initialization -reset(gpuDevice(1)); %% Input Parser % ACDC: acceleration & deceleration correction (default: true) % DPS: Detector Point Spread correction (default: true) % SC: Scatter Correction (default: true) % BH: Beam Hardening correction (default: false, due to Bowtie, BH correction is not required at all.) -[tag_ACDC, tag_DPS, tag_SC, tag_BH] = parse_inputs(varargin{:}); +[tag_ACDC, tag_DPS, tag_SC, tag_BH, gpuids] = parse_inputs(varargin{:}); + +%% GPU initialization +reset(gpuDevice(gpuids.devices(0)+1)); %% Load geometry [geo, ScanXML] = GeometryFromXML(datafolder); @@ -40,7 +41,7 @@ % Detector point scatter correction if(tag_DPS) disp('Proj DPS: ') - proj = DetectorPointScatterCorrection(proj, geo, ScCalib); + proj = DetectorPointScatterCorrection(proj, geo, ScCalib, gpuids); end %% Load blank scan @@ -59,14 +60,14 @@ end %% Airnorm and Logarithmic Normalization -proj_lg = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm); +proj_lg = LogNormal(proj, angles, airnorm, Blk, Sec, BlkAirNorm, gpuids); disp('Log Normalization is completed.') % remove anomolies proj_lg = EnforcePositive(proj_lg); %% Beam Hardening correction is applied (kind of slow) if(tag_BH) - [proj_lg, ~] = BHCorrection(datafolder, geo, ScanXML, proj_lg); + [proj_lg, ~] = BHCorrection(datafolder, geo, ScanXML, proj_lg, gpuids); end %% Remove anomalies @@ -85,7 +86,7 @@ end -function [tag_ACDC, tag_DPS, tag_SC, tag_BH] = parse_inputs(varargin) +function [tag_ACDC, tag_DPS, tag_SC, tag_BH, gpuids] = parse_inputs(varargin) % create input parser p = inputParser; % add optional parameters @@ -95,6 +96,8 @@ % BH performance is very lame for unclear reason addParameter(p,'bh', false); +addParameters(p,'gpuids',GpuIds()) + %execute parse(p,varargin{:}); %extract @@ -102,4 +105,5 @@ tag_DPS=p.Results.dps; tag_SC=p.Results.sc; tag_BH=p.Results.bh; +gpuids=p.Results.gpuids; end