WIP: Pull request for merging LMU modifications

applications/pctbinning/pctbinning.cxx

@@ -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)
   {
@@ -128,6 +130,37 @@ main(int argc, char * argv[])
     TRY_AND_EXIT_ON_ITK_EXCEPTION(cwriter->Update())
   }
 
+  if (args_info.variance_given)
+  {
+
+    auto varianceHoleFilter = pct::HoleFillingImageFilter<OutputImageType>::New();
+    if (args_info.variance_given && args_info.fillvariance_flag)
+    {
+      varianceHoleFilter->SetInput(projection->GetVariance());
+      TRY_AND_EXIT_ON_ITK_EXCEPTION(varianceHoleFilter->Update());
+    }
+
+    typedef itk::ChangeInformationImageFilter<ProjectionFilter::VarianceImageType> CIIVarType;
+    CIIVarType::Pointer                                                            ciiVar = CIIVarType::New();
+    if (args_info.fillvariance_flag)
+      ciiVar->SetInput(varianceHoleFilter->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

applications/pctbinning/pctbinning.ggo

@@ -6,9 +6,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
@@ -18,6 +20,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

applications/pctfdk/pctfdk.cxx

@@ -5,6 +5,7 @@
 #include "rtkProjectionsReader.h"
 #include "pctDDParkerShortScanImageFilter.h"
 #include "pctFDKDDConeBeamReconstructionFilter.h"
+#include "pctFDKDDConeBeamVarianceReconstructionFilter.h"
 
 #include <itkRegularExpressionSeriesFileNames.h>
 #include <itkImageFileWriter.h>
@@ -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());

applications/pctfdk/pctfdk.ggo

@@ -13,6 +13,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

applications/pctlomalinda/pctlomalinda.py

@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+import argparse
+import uproot
+import itk
+from itk import PCT as pct
+
+import numpy as np
+
+
+def build_parser():
+
+    parser = pct.PCTArgumentParser(description="Convert Loma Linda data to PCT data")
+    parser.add_argument(
+        "-i", "--input", help="Root phase space file of particles", required=True
+    )
+    parser.add_argument("-o", "--output", help="Output file name", required=True)
+    parser.add_argument(
+        "--plane-in",
+        help="Plane position of incoming protons",
+        required=True,
+        type=float,
+    )
+    parser.add_argument(
+        "--plane-out",
+        help="Plane position of outgoing protons",
+        required=True,
+        type=float,
+    )
+    parser.add_argument(
+        "--min-run", help="Minimum run (inclusive)", default=0, type=int
+    )
+    parser.add_argument(
+        "--max-run", help="Maximum run (exclusive)", default=1e6, type=int
+    )
+    parser.add_argument(
+        "--verbose", "-v", help="Verbose execution", default=False, action="store_true"
+    )
+    parser.add_argument(
+        "--ps", help="Name of tree in input phase space", default="PhaseSpace"
+    )
+
+    return parser
+
+
+def process(args_info: argparse.Namespace):
+
+    if args_info.verbose:
+
+        def verbose(message):
+            print(message)
+
+    else:
+
+        def verbose(message):
+            pass
+
+    tree = uproot.open(args_info.input)[args_info.ps]
+    pairs = tree.arrays(library="np")
+    verbose("Read input phase space:\n" + str(pairs))
+
+    # should match what it in the file, no checks is performed
+    pairs["u_hit1"] = np.full_like(pairs["u_hit1"], args_info.plane_in)
+    pairs["u_hit2"] = np.full_like(pairs["u_hit2"], args_info.plane_out)
+
+    verbose("Filter RunIDs…")
+    # "projection_angle" acts as the run ID here
+    interception = np.logical_and(
+        pairs["projection_angle"] < args_info.max_run,
+        pairs["projection_angle"] >= args_info.min_run,
+    )
+    for k, v in pairs.items():
+        pairs[k] = v[interception]
+
+    verbose("Calculating directions…")
+    dir_in_t = pairs["t_hit1"] - pairs["t_hit0"]
+    dir_in_v = pairs["v_hit1"] - pairs["v_hit0"]
+    dir_in_u = pairs["u_hit1"] - pairs["u_hit0"]
+    norm_in = np.sqrt(dir_in_t**2 + dir_in_v**2 + dir_in_u**2)
+    dir_out_t = pairs["t_hit3"] - pairs["t_hit2"]
+    dir_out_v = pairs["v_hit3"] - pairs["v_hit2"]
+    dir_out_u = pairs["u_hit3"] - pairs["u_hit2"]
+    norm_out = np.sqrt(dir_out_t**2 + dir_out_v**2 + dir_out_u**2)
+
+    ComponentType = itk.ctype("float")
+    PixelType = itk.Vector[ComponentType, 3]
+    ImageType = itk.Image[PixelType, 2]
+
+    runs = np.unique(pairs["projection_angle"])
+    number_of_runs = len(runs)
+    verbose("Identified number of runs: " + str(number_of_runs))
+    run_range = range(args_info.min_run, min(number_of_runs, args_info.max_run))
+    for r in run_range:
+        ps_run = pairs["projection_angle"] == runs[r]
+        len_ps_run = ps_run.sum()
+        if len_ps_run == 0:
+            continue
+
+        ps_np = np.empty(shape=(len_ps_run, 5, 3), dtype=np.float32)
+        ps_np[:, 0, 0] = pairs["t_hit1"][ps_run]
+        ps_np[:, 0, 1] = pairs["v_hit1"][ps_run]
+        ps_np[:, 0, 2] = pairs["u_hit1"][ps_run]
+        ps_np[:, 1, 0] = pairs["t_hit2"][ps_run]
+        ps_np[:, 1, 1] = pairs["v_hit2"][ps_run]
+        ps_np[:, 1, 2] = pairs["u_hit2"][ps_run]
+        ps_np[:, 2, 0] = dir_in_t[ps_run] / norm_in[ps_run]
+        ps_np[:, 2, 1] = dir_in_v[ps_run] / norm_in[ps_run]
+        ps_np[:, 2, 2] = dir_in_u[ps_run] / norm_in[ps_run]
+        ps_np[:, 3, 0] = dir_in_t[ps_run] / norm_out[ps_run]
+        ps_np[:, 3, 1] = dir_in_v[ps_run] / norm_out[ps_run]
+        ps_np[:, 3, 2] = dir_in_u[ps_run] / norm_out[ps_run]
+        ps_np[:, 4, 0] = 0.0  # WEPL is already present in the ROOT file
+        ps_np[:, 4, 1] = pairs["calculated_WEPL"][ps_run]
+        ps_np[:, 4, 2] = pairs["timestamp"][ps_run]
+
+        df_itk = itk.GetImageFromArray(ps_np, ttype=ImageType)
+
+        output_file = args_info.output.replace(".", f"{r:04d}.")
+        itk.imwrite(df_itk, output_file)
+        verbose(f"Wrote file {output_file}.")
+
+
+def main(argv=None):
+    parser = build_parser()
+    args_info = parser.parse_args(argv)
+    process(args_info)
+
+
+if __name__ == "__main__":
+    main()

applications/pctpaircuts/pctpaircuts.cxx

@@ -12,6 +12,7 @@
 
 #include <itkImageFileWriter.h>
 #include <itkRegularExpressionSeriesFileNames.h>
+#include <itkMersenneTwisterRandomVariateGenerator.h>
 #include <itkTimeProbe.h>
 
 #define PAIRS_IN_RAM 1000000
@@ -55,6 +56,14 @@ main(int argc, char * argv[])
   ConvFunc = new pct::Functor::IntegratedBetheBlochProtonStoppingPowerInverse<float, double>(
     68.9984 * CLHEP::eV, 600. * CLHEP::MeV, 0.1 * CLHEP::keV);
 
+  // RNG for fluence
+  using GeneratorType = itk::Statistics::MersenneTwisterRandomVariateGenerator;
+  auto generator = GeneratorType::New();
+  if (args_info.seed_given)
+    generator->Initialize(args_info.seed_arg);
+  else
+    generator->Initialize();
+
   // Read pairs
   using VectorType = itk::Vector<float, 3>;
   using PairsImageType = itk::Image<VectorType, 2>;
@@ -283,6 +292,41 @@ main(int argc, char * argv[])
         ++it;
       }
 
+      // Fluence filtering
+      if (args_info.fluence_given)
+      {
+        if (generator->GetVariate() >= args_info.fluence_arg)
+          continue;
+      }
+
+      // Check ROI intersection
+      if (args_info.roiradius_given)
+      {
+        float radius = args_info.roiradius_arg;
+        float dx = pOut[0] - pIn[0];
+        float dz = pOut[2] - pIn[2];
+        float dr_sq = (dx * dx) + (dz * dz);
+        float d = (pIn[0] * pOut[2]) - (pOut[0] * pIn[2]);
+        float delta = (radius * radius * dr_sq) - (d * d);
+        if (delta < 0.)
+          continue;
+      }
+
+      // Filter low or high WEPLs
+      if (data[0] == 0.)
+      {
+        if ((args_info.minwepl_given && (data[1] < args_info.minwepl_arg)) ||
+            (args_info.maxwepl_given && (args_info.maxwepl_arg < data[1])))
+          continue;
+      }
+      else
+      {
+        if (args_info.minwepl_given || args_info.maxwepl_given)
+          std::cerr << "Cannot specify WEPL bounds because the input file does not contain WEPL data. Aborting."
+                    << std::endl;
+        return EXIT_FAILURE;
+      }
+
       static double mag =
         (args_info.source_arg == 0.) ? 1. : (args_info.source_arg - pOut[2]) / (args_info.source_arg - pIn[2]);
 

applications/pctpaircuts/pctpaircuts.ggo

@@ -12,6 +12,11 @@ option "robustopt" - "Use newer options for robust cut."            int       no
 option "primaries" - "Consider only primary protons"                flag      off
 option "nonuclear" - "Consider only primary protons without nuclear interactions"                flag      off
 option "wet" - "Write WET instead of initial energy"                flag      off
+option "roiradius" - "Radius of the ROI in mm"                      double    no
+option "minwepl"   - "Minimum WEPL"                                 double    no
+option "maxwepl"   - "Maximum WEPL"                                 double    no
+option "fluence"   - "Fluence level as fraction of full fluence"    double    no
+option "seed"      - "Random seed (for reproducibility)"            int       no
 
 section "Projections parameters"
 option "origin"    - "Origin (default=centered)" double multiple no

include/pctBetheBlochFunctor.h

@@ -29,8 +29,13 @@ class BetheBlochProtonStoppingPower
   ~BetheBlochProtonStoppingPower() {}
 
   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;
@@ -55,10 +60,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);
@@ -71,7 +78,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++)
       {
@@ -118,6 +125,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse
   std::vector<TOutput>                            m_Length;
 
   double               m_BinSize;
+  std::string          m_Particle;
   std::vector<TOutput> m_LUT;
 };
 } // end namespace Functor