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/pctfilterprotons/pctfilterprotons.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+import argparse
+import itk
+from itk import PCT as pct
+import numpy as np
+
+
+def build_parser():
+
+    parser = pct.PCTArgumentParser(
+        description="Filter protons according to various criteria"
+    )
+    parser.add_argument("-i", "--input", help="Input list-mode file", required=True)
+    parser.add_argument("-o", "--output", help="Output list-mode file", required=True)
+    parser.add_argument("--roi-radius", help="Radius of the ROI in mm", type=float)
+    parser.add_argument(
+        "--fluence", help="Fluence level as fraction of full fluence", type=float
+    )
+    parser.add_argument("--min-wepl", help="Minimum WEPL", type=float)
+    parser.add_argument("--max-wepl", help="Maximum WEPL", type=float)
+    parser.add_argument("--seed", help="Random seed (for reproducibility)", type=int)
+    parser.add_argument(
+        "--verbose", "-v", help="Verbose execution", default=False, action="store_true"
+    )
+
+    return parser
+
+
+def process(args_info: argparse.Namespace):
+
+    if args_info.verbose:
+
+        def verbose(message):
+            print(message)
+
+    else:
+
+        def verbose(message):
+            pass
+
+    pairs_itk = itk.imread(args_info.input)
+    pairs = itk.GetArrayFromImage(pairs_itk)
+
+    u_in = np.s_[:, 0, 0]
+    w_in = np.s_[:, 0, 2]
+    u_out = np.s_[:, 1, 0]
+    w_out = np.s_[:, 1, 2]
+    e_in = np.s_[:, 4, 0]
+    wepl = np.s_[:, 4, 1]
+
+    if np.any(pairs[e_in] != 0.0) and (
+        args_info.min_wepl is not None or args_info.max_wepl is not None
+    ):
+        raise ValueError(
+            "Cannot filter WEPLs because input file does not contain WEPL (e_in != 0)! Aborting."
+        )
+
+    interception = np.full_like(pairs[wepl], True)
+
+    if args_info.roi_radius is not None:
+        # Circular ROI intersection
+        verbose("Filtering ions outside of the ROI…")
+        radius = args_info.roi_radius
+        du = pairs[u_out] - pairs[u_in]
+        dt = pairs[w_out] - pairs[w_in]
+        dr_sq = (du * du) + (dt * dt)
+        d = (pairs[u_in] * pairs[w_out]) - (pairs[u_out] * pairs[w_in])
+        delta = (radius * radius * dr_sq) - (d * d)
+        interception = np.logical_and(interception, delta >= 0.0)
+
+    if args_info.fluence is not None:
+        # Fluence filtering
+        verbose("Applying fluence level…")
+        rng = np.random.default_rng(args_info.seed)
+        fluence_filter = rng.random(interception.shape)
+        interception = np.logical_and(interception, fluence_filter < args_info.fluence)
+
+    if args_info.min_wepl is not None:
+        verbose("Filtering low WEPLs…")
+        interception = np.logical_and(interception, pairs[wepl] >= args_info.min_wepl)
+
+    if args_info.max_wepl is not None:
+        verbose("Filtering high WEPLs…")
+        interception = np.logical_and(interception, pairs[wepl] <= args_info.max_wepl)
+
+    verbose("Applying interception filter…")
+
+    ComponentType = itk.ctype("float")
+    PixelType = itk.Vector[ComponentType, 3]
+    ImageType = itk.Image[PixelType, 2]
+
+    output_itk = itk.GetImageFromArray(pairs[interception], ttype=ImageType)
+    itk.imwrite(output_itk, args_info.output)
+    verbose(f"Wrote file {args_info.output}.")
+
+
+def main(argv=None):
+    parser = build_parser()
+    args_info = parser.parse_args(argv)
+    process(args_info)
+
+
+if __name__ == "__main__":
+    main()

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()

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