Skip to content

WIP: Pull request for merging LMU modifications #41

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SimonRit wants to merge 6 commits into
base: main
Choose a base branch
from
Open

WIP: Pull request for merging LMU modifications #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0d5f213
ENH: Implemented helium coefficients
Oct 8, 2025
c47f1be
ENH: Add LMU variance reconstruction code
Oct 8, 2025
037ab04
ENH: Limit the filling filter to 20 iterations
Oct 8, 2025
d045df1
ENH: Add intersection option to LMU fluence modulated PCT
Oct 8, 2025
ad9f0f0
COMP: Fix compilation for recent Geant4
Oct 8, 2025
2292bdf
ENH: Add time to list mode files produced from Loma Linda data
Oct 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions applications/pctbinning/pctbinning.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,8 @@ main(int argc, char * argv[])
projection->SetRobust(args_info.robust_flag);
projection->SetComputeScattering(args_info.scatwepl_given);
projection->SetComputeNoise(args_info.noise_given);
projection->SetComputeVariance(args_info.variance_given);
projection->SetParticle(args_info.particle_arg);

if (args_info.quadricIn_given)
{
Expand Down Expand Up @@ -129,6 +131,38 @@ main(int argc, char * argv[])
TRY_AND_EXIT_ON_ITK_EXCEPTION(cwriter->Update())
}

if (args_info.variance_given)
{

SmallHoleFiller<ProjectionFilter::VarianceImageType> fillerVar;
if (args_info.variance_given && args_info.fillvariance_flag)
{
fillerVar.SetImage(projection->GetVariance());
fillerVar.SetHolePixel(0.);
fillerVar.Fill();
}

typedef itk::ChangeInformationImageFilter<ProjectionFilter::VarianceImageType> CIIVarType;
CIIVarType::Pointer ciiVar = CIIVarType::New();
if (args_info.fillvariance_flag)
ciiVar->SetInput(fillerVar.GetOutput());
else
ciiVar->SetInput(projection->GetVariance());
ciiVar->ChangeOriginOn();
ciiVar->ChangeDirectionOn();
ciiVar->ChangeSpacingOn();
ciiVar->SetOutputDirection(projection->GetOutput()->GetDirection());
ciiVar->SetOutputOrigin(projection->GetOutput()->GetOrigin());
ciiVar->SetOutputSpacing(projection->GetOutput()->GetSpacing());

// Write
typedef itk::ImageFileWriter<ProjectionFilter::VarianceImageType> VarianceWriterType;
VarianceWriterType::Pointer vwriter = VarianceWriterType::New();
vwriter->SetFileName(args_info.variance_arg);
vwriter->SetInput(ciiVar->GetOutput());
TRY_AND_EXIT_ON_ITK_EXCEPTION(vwriter->Update())
}

if (args_info.scatwepl_given)
{
// Write
Expand Down
3 changes: 3 additions & 0 deletions applications/pctbinning/pctbinning.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,11 @@ option "input" i "Input file name containing the proton pairs"
option "output" o "Output file name" string no
option "elosswepl" - "Output file name (alias for --output)" string no
option "count" c "Image of count of proton pairs per pixel" string no
option "variance" - "Image of variance of WEPL values per pixel" string no
option "scatwepl" - "Image of scattering WEPL of proton pairs per pixel" string no
option "noise" - "Image of WEPL variance per pixel" string no
option "robust" r "Use robust estimation of scattering using 19.1 %ile." flag off
option "particle" p "Particle used for imaging (proton or alpha)" string no default="proton"

option "source" s "Source position" double no default="0."
option "quadricIn" - "Parameters of the entrance quadric support function, see http://education.siggraph.org/static/HyperGraph/raytrace/rtinter4.htm" double multiple no
Expand All @@ -19,6 +21,7 @@ option "mlptrackeruncert" - "Consider tracker uncertainties in MLP [Krah 2018
option "mlppolydeg" - "Degree of the polynomial to approximate 1/beta^2p^2" int no default="5"
option "ionpot" - "Ionization potential used in the reconstruction in eV" double no default="68.9984"
option "fill" - "Fill holes, i.e. pixels that were not hit by protons" flag off
option "fillvariance" - "Fill holes, i.e. pixels that were not hit by protons (for variance projections)" flag off
option "trackerresolution" - "Tracker resolution in mm" double no
option "trackerspacing" - "Tracker pair spacing in mm" double no
option "materialbudget" - "Material budget x/X0 of tracker" double no
Expand Down
13 changes: 12 additions & 1 deletion applications/pctfdk/pctfdk.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "rtkProjectionsReader.h"
#include "pctDDParkerShortScanImageFilter.h"
#include "pctFDKDDConeBeamReconstructionFilter.h"
#include "pctFDKDDConeBeamVarianceReconstructionFilter.h"

#include <itkRegularExpressionSeriesFileNames.h>
#include <itkImageFileWriter.h>
Expand Down Expand Up @@ -62,7 +63,17 @@ main(int argc, char * argv[])
rtk::SetConstantImageSourceFromGgo<ConstantImageSourceType, args_info_pctfdk>(constantImageSource, args_info);

// FDK reconstruction filtering
auto feldkamp = pct::FDKDDConeBeamReconstructionFilter<OutputImageType>::New();
using FDKCPUType = pct::FDKDDConeBeamReconstructionFilter<OutputImageType>;
FDKCPUType::Pointer feldkamp = nullptr;
if (args_info.variance_flag)
{
feldkamp = pct::FDKDDConeBeamVarianceReconstructionFilter<OutputImageType>::New();
}
else
{
feldkamp = FDKCPUType::New();
}

feldkamp->SetInput(0, constantImageSource->GetOutput());
feldkamp->SetProjectionStack(pssf->GetOutput());
feldkamp->SetGeometry(geometryReader->GetOutputObject());
Expand Down
1 change: 1 addition & 0 deletions applications/pctfdk/pctfdk.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ section "Ramp filter"
option "pad" - "Data padding parameter to correct for truncation" double no default="0.0"
option "hann" - "Cut frequency for hann window in ]0,1] (0.0 disables it)" double no default="0.0"
option "hannY" - "Cut frequency for hann window in ]0,1] (0.0 disables it)" double no default="0.0"
option "variance" - "Use variance reconstruction kernel" flag off

section "Volume properties"
option "origin" - "Origin (default=centered)" double multiple no
Expand Down
4 changes: 2 additions & 2 deletions applications/pctschulte/pctschulte.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,8 @@ main(int argc, char * argv[])
{
os << E / CLHEP::MeV;
os << ',';
os << bb.GetValue(E * CLHEP::MeV, args_info.ionpot_arg * CLHEP::eV) / (1000. * CLHEP::kg / CLHEP::m3) * CLHEP::g /
(CLHEP::MeV * CLHEP::cm2);
os << bb.GetValue(E * CLHEP::MeV, args_info.ionpot_arg * CLHEP::eV, "proton") / (1000. * CLHEP::kg / CLHEP::m3) *
CLHEP::g / (CLHEP::MeV * CLHEP::cm2);
os << std::endl;
}
os.close();
Expand Down
5 changes: 3 additions & 2 deletions geant4/DetectorConstruction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,8 +145,9 @@ DetectorConstruction::Construct()
<< " Total Length(mm)= " << 2.0 * targetZ / mm << " ###" << G4endl;

// colors
G4VisAttributes zero = G4VisAttributes::Invisible;
fLogicWorld->SetVisAttributes(&zero);
// GD: The two lines below are not compatible with Geant4.11
// G4VisAttributes zero = G4VisAttributes::Invisible;
// fLogicWorld->SetVisAttributes(&zero);

G4VisAttributes regWcolor(G4Colour(0.3, 0.3, 0.3));
fLogicCheck->SetVisAttributes(&regWcolor);
Expand Down
6 changes: 3 additions & 3 deletions geant4/G4EmUserPhysics.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@
#include "G4LossTableManager.hh"
#include "G4VEnergyLossProcess.hh"
#include "G4VProcess.hh"
#include "G4EmProcessOptions.hh"
// #include "G4EmProcessOptions.hh" // GD: obsolete in 11.3
#include "G4AntiProton.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
Expand Down Expand Up @@ -95,8 +95,8 @@ G4EmUserPhysics::ConstructProcess()
}
}

G4EmProcessOptions opt;
opt.SetVerbose(fVerbose);
// G4EmProcessOptions opt;
// opt.SetVerbose(fVerbose);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4 changes: 2 additions & 2 deletions geant4/PhysicsList.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@
#include "G4IonBinaryCascadePhysics.hh"
#include "G4IonPhysics.hh"
#include "G4EmExtraPhysics.hh"
#include "G4EmProcessOptions.hh"
// #include "G4EmProcessOptions.hh" // GD: obsolete in 11.3

#include "G4HadronPhysicsFTFP_BERT.hh"
#include "G4HadronPhysicsFTFP_BERT_HP.hh"
Expand All @@ -88,7 +88,7 @@
#include "G4Proton.hh"

#include "G4SystemOfUnits.hh"

#include "G4UnitsTable.hh" // GD: needed in 11.3
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

PhysicsList::PhysicsList()
Expand Down
2 changes: 1 addition & 1 deletion include/SmallHoleFiller.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ SmallHoleFiller<TImage>::Fill()
// Initialize by setting the output image to the input image.
DeepCopy<TImage>(this->Image, this->Output);
unsigned int numberOfIterations = 0;
while (HasEmptyPixels())
while (HasEmptyPixels() && numberOfIterations < 20)
{
std::cout << "Iteration " << numberOfIterations << "..." << std::endl;
Iterate();
Expand Down
45 changes: 34 additions & 11 deletions include/pctBetheBlochFunctor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#ifdef PCT_GEANT4
# include "geant4/pctGeant4.h"
# include <G4Material.hh>
# include <G4Alpha.hh>
# include <G4Proton.hh>
# include <G4BetheBlochModel.hh>
# include <G4NistManager.hh>
Expand Down Expand Up @@ -43,17 +44,36 @@ class BetheBlochProtonStoppingPower

#ifdef PCT_GEANT4
TOutput
GetValue(const TInput e, const double itkNotUsed(I)) const
GetValue(const TInput e, const double itkNotUsed(I), const std::string particle) const
{
return m_G4BetheBlochModel->ComputeDEDXPerVolume(
G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), // G4Material::GetMaterial("Water"),
G4Proton::Proton(),
G4double(e),
G4double(1 * CLHEP::km));
if (particle == "proton")
{
return m_G4BetheBlochModel->ComputeDEDXPerVolume(G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"),
G4Proton::Proton(),
G4double(e),
G4double(1 * CLHEP::km));
}
else if (particle == "alpha")
{
return m_G4BetheBlochModel->ComputeDEDXPerVolume(G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"),
G4Alpha::Alpha(),
G4double(e),
G4double(1 * CLHEP::km));
}
else
{
throw std::invalid_argument("Invalid particle type in Bethe Bloch Functor.");
}

#else
TOutput
GetValue(const TInput e, const double I) const
GetValue(const TInput e, const double I, const std::string particle) const
{
if (particle != "proton")
{
// only implemented for protons
throw std::invalid_argument("Invalid particle type in Bethe Block Functor.");
}
/** Physical constants */
static const double K = 4. * CLHEP::pi * CLHEP::classic_electr_radius * CLHEP::classic_electr_radius *
CLHEP::electron_mass_c2 * 3.343e+23 / CLHEP::cm3;
Expand Down Expand Up @@ -89,10 +109,12 @@ template <class TInput, class TOutput>
class IntegratedBetheBlochProtonStoppingPowerInverse
{
public:
IntegratedBetheBlochProtonStoppingPowerInverse(const double I,
const double maxEnergy,
const double binSize = 1. * CLHEP::keV)
IntegratedBetheBlochProtonStoppingPowerInverse(const double I,
const double maxEnergy,
const double binSize = 1. * CLHEP::keV,
const std::string particle = "proton")
: m_BinSize(binSize)
, m_Particle(particle)
{
unsigned int lowBinLimit, numberOfBins;
lowBinLimit = itk::Math::Ceil<unsigned int, double>(m_S.GetLowEnergyLimit() / m_BinSize);
Expand All @@ -105,7 +127,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse
}
for (unsigned int i = lowBinLimit; i < numberOfBins; i++)
{
m_LUT[i] = m_LUT[i - 1] + binSize / m_S.GetValue(TOutput(i) * binSize, I);
m_LUT[i] = m_LUT[i - 1] + binSize / m_S.GetValue(TOutput(i) * binSize, I, m_Particle);
// Create inverse lut, i.e., get energy from length in water
for (unsigned j = m_Length.size(); j < unsigned(m_LUT[i] * CLHEP::mm) + 1; j++)
{
Expand Down Expand Up @@ -152,6 +174,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse
std::vector<TOutput> m_Length;

double m_BinSize;
std::string m_Particle;
std::vector<TOutput> m_LUT;
};
} // end namespace Functor
Expand Down
68 changes: 68 additions & 0 deletions include/pctFDKDDConeBeamVarianceReconstructionFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
#ifndef __pctFDKDDConeBeamVarianceReconstructionFilter_h
#define __pctFDKDDConeBeamVarianceReconstructionFilter_h

#include "pctFFTVarianceImageFilter.h"
#include "pctFDKDDConeBeamReconstructionFilter.h"
#include "pctFFTVarianceImageFilter.h"

#include <itkExtractImageFilter.h>
#include <itkTimeProbe.h>

/** \class FDKDDConeBeamReconstructionFilter
* TODO
*
* \author Jannis Dickmann
*/
namespace pct
{

template <class TInputImage, class TOutputImage = TInputImage, class TFFTPrecision = double>
class ITK_EXPORT FDKDDConeBeamVarianceReconstructionFilter
: public pct::FDKDDConeBeamReconstructionFilter<TInputImage, TOutputImage>
{
public:
typedef pct::FDKDDConeBeamReconstructionFilter<TInputImage, TOutputImage> Baseclass;

/** Standard class typedefs. */
typedef FDKDDConeBeamVarianceReconstructionFilter Self;
typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;

/** Some convenient typedefs. */
typedef typename Baseclass::InputImageType InputImageType;
typedef typename Baseclass::OutputImageType OutputImageType;

typedef typename Baseclass::ProjectionStackType ProjectionStackType;
typedef typename Baseclass::ProjectionStackPointer ProjectionStackPointer;

/** Typedefs of each subfilter of this composite filter */
typedef typename Baseclass::ExtractFilterType ExtractFilterType;
typedef typename Baseclass::WeightFilterType WeightFilterType;
typedef pct::FFTVarianceImageFilter<ProjectionStackType, ProjectionStackType, TFFTPrecision> VarianceFilterType;
typedef typename Baseclass::BackProjectionFilterType BackProjectionFilterType;

/** Standard New method. */
itkNewMacro(Self);

/** Runtime information support. */
itkOverrideGetNameOfClassMacro(FDKDDConeBeamVarianceReconstructionFilter);

protected:
FDKDDConeBeamVarianceReconstructionFilter();
~FDKDDConeBeamVarianceReconstructionFilter() {}

private:
// purposely not implemented
FDKDDConeBeamVarianceReconstructionFilter(const Self &);
void
operator=(const Self &);
}; // end of class

} // end namespace pct

#ifndef ITK_MANUAL_INSTANTIATION
# include "pctFDKDDConeBeamVarianceReconstructionFilter.hxx"
#endif

#endif
33 changes: 33 additions & 0 deletions include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
namespace pct
{

template <class TInputImage, class TOutputImage, class TFFTPrecision>
FDKDDConeBeamVarianceReconstructionFilter<TInputImage, TOutputImage, TFFTPrecision>::
FDKDDConeBeamVarianceReconstructionFilter()
{
this->SetNumberOfRequiredInputs(1);

// Create each filter of the composite filter
this->m_ExtractFilter = ExtractFilterType::New();
this->m_WeightFilter = WeightFilterType::New();
this->m_RampFilter = VarianceFilterType::New();
this->m_BackProjectionFilter = BackProjectionFilterType::New();

// Permanent internal connections
this->m_WeightFilter->SetInput(this->m_ExtractFilter->GetOutput());
this->m_RampFilter->SetInput(this->m_WeightFilter->GetOutput());
this->m_BackProjectionFilter->SetProjectionStack(this->m_RampFilter->GetOutput());

// Default parameters
#if ITK_VERSION_MAJOR >= 4
this->m_ExtractFilter->SetDirectionCollapseToSubmatrix();
#endif
this->m_WeightFilter->InPlaceOn();
this->m_BackProjectionFilter->InPlaceOn();
this->m_BackProjectionFilter->SetTranspose(true);

this->m_WeightFilter->SetVarianceReconstruction(true);
}


} // end namespace pct
6 changes: 6 additions & 0 deletions include/pctFDKDDWeightProjectionFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,10 @@ class ITK_TEMPLATE_EXPORT FDKDDWeightProjectionFilter : public itk::InPlaceImage
itkGetMacro(Geometry, rtk::ThreeDCircularProjectionGeometry::Pointer);
itkSetMacro(Geometry, rtk::ThreeDCircularProjectionGeometry::Pointer);

/** Get/ Set variance mode */
itkGetMacro(VarianceReconstruction, bool);
itkSetMacro(VarianceReconstruction, bool);

protected:
FDKDDWeightProjectionFilter() {}
~FDKDDWeightProjectionFilter() {}
Expand All @@ -55,6 +59,8 @@ class ITK_TEMPLATE_EXPORT FDKDDWeightProjectionFilter : public itk::InPlaceImage
virtual void
DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;

bool m_VarianceReconstruction = false;

private:
FDKDDWeightProjectionFilter(const Self &); // purposely not implemented
void
Expand Down
5 changes: 5 additions & 0 deletions include/pctFDKDDWeightProjectionFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,11 @@ FDKDDWeightProjectionFilter<TInputImage, TOutputImage>::BeforeThreadedGenerateDa
const double rampFactor = sdd / (2. * sid);
m_AngularWeightsAndRampFactor[k] *= rampFactor;
}
if (this->m_VarianceReconstruction)
{
// square if in variance mode
m_AngularWeightsAndRampFactor[k] = m_AngularWeightsAndRampFactor[k] * m_AngularWeightsAndRampFactor[k];
}
}
}

Expand Down
Loading
Loading