diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index 5f51925..aa20cc8 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/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) { @@ -129,6 +131,38 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(cwriter->Update()) } + if (args_info.variance_given) + { + + SmallHoleFiller fillerVar; + if (args_info.variance_given && args_info.fillvariance_flag) + { + fillerVar.SetImage(projection->GetVariance()); + fillerVar.SetHolePixel(0.); + fillerVar.Fill(); + } + + typedef itk::ChangeInformationImageFilter 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 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 diff --git a/applications/pctbinning/pctbinning.ggo b/applications/pctbinning/pctbinning.ggo index 86a4a06..bcc54fe 100644 --- a/applications/pctbinning/pctbinning.ggo +++ b/applications/pctbinning/pctbinning.ggo @@ -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 @@ -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 diff --git a/applications/pctfdk/pctfdk.cxx b/applications/pctfdk/pctfdk.cxx index 6ce02bb..e095c32 100644 --- a/applications/pctfdk/pctfdk.cxx +++ b/applications/pctfdk/pctfdk.cxx @@ -5,6 +5,7 @@ #include "rtkProjectionsReader.h" #include "pctDDParkerShortScanImageFilter.h" #include "pctFDKDDConeBeamReconstructionFilter.h" +#include "pctFDKDDConeBeamVarianceReconstructionFilter.h" #include #include @@ -62,7 +63,17 @@ main(int argc, char * argv[]) rtk::SetConstantImageSourceFromGgo(constantImageSource, args_info); // FDK reconstruction filtering - auto feldkamp = pct::FDKDDConeBeamReconstructionFilter::New(); + using FDKCPUType = pct::FDKDDConeBeamReconstructionFilter; + FDKCPUType::Pointer feldkamp = nullptr; + if (args_info.variance_flag) + { + feldkamp = pct::FDKDDConeBeamVarianceReconstructionFilter::New(); + } + else + { + feldkamp = FDKCPUType::New(); + } + feldkamp->SetInput(0, constantImageSource->GetOutput()); feldkamp->SetProjectionStack(pssf->GetOutput()); feldkamp->SetGeometry(geometryReader->GetOutputObject()); diff --git a/applications/pctfdk/pctfdk.ggo b/applications/pctfdk/pctfdk.ggo index 7c51685..77e1b14 100644 --- a/applications/pctfdk/pctfdk.ggo +++ b/applications/pctfdk/pctfdk.ggo @@ -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 diff --git a/applications/pctschulte/pctschulte.cxx b/applications/pctschulte/pctschulte.cxx index b71d0e8..3848bd9 100644 --- a/applications/pctschulte/pctschulte.cxx +++ b/applications/pctschulte/pctschulte.cxx @@ -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(); diff --git a/geant4/DetectorConstruction.cc b/geant4/DetectorConstruction.cc index ccaf11a..a5b94ee 100644 --- a/geant4/DetectorConstruction.cc +++ b/geant4/DetectorConstruction.cc @@ -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(®Wcolor); diff --git a/geant4/G4EmUserPhysics.cc b/geant4/G4EmUserPhysics.cc index 278d923..d676f51 100644 --- a/geant4/G4EmUserPhysics.cc +++ b/geant4/G4EmUserPhysics.cc @@ -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" @@ -95,8 +95,8 @@ G4EmUserPhysics::ConstructProcess() } } - G4EmProcessOptions opt; - opt.SetVerbose(fVerbose); + // G4EmProcessOptions opt; + // opt.SetVerbose(fVerbose); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/PhysicsList.cc b/geant4/PhysicsList.cc index 248fec8..cccbaed 100644 --- a/geant4/PhysicsList.cc +++ b/geant4/PhysicsList.cc @@ -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" @@ -88,7 +88,7 @@ #include "G4Proton.hh" #include "G4SystemOfUnits.hh" - +#include "G4UnitsTable.hh" // GD: needed in 11.3 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... PhysicsList::PhysicsList() diff --git a/include/SmallHoleFiller.hxx b/include/SmallHoleFiller.hxx index 5de6d7e..6a7067b 100644 --- a/include/SmallHoleFiller.hxx +++ b/include/SmallHoleFiller.hxx @@ -51,7 +51,7 @@ SmallHoleFiller::Fill() // Initialize by setting the output image to the input image. DeepCopy(this->Image, this->Output); unsigned int numberOfIterations = 0; - while (HasEmptyPixels()) + while (HasEmptyPixels() && numberOfIterations < 20) { std::cout << "Iteration " << numberOfIterations << "..." << std::endl; Iterate(); diff --git a/include/pctBetheBlochFunctor.h b/include/pctBetheBlochFunctor.h index e1ea5bb..b6c449c 100644 --- a/include/pctBetheBlochFunctor.h +++ b/include/pctBetheBlochFunctor.h @@ -7,6 +7,7 @@ #ifdef PCT_GEANT4 # include "geant4/pctGeant4.h" # include +# include # include # include # include @@ -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; @@ -89,10 +109,12 @@ template 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(m_S.GetLowEnergyLimit() / m_BinSize); @@ -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++) { @@ -152,6 +174,7 @@ class IntegratedBetheBlochProtonStoppingPowerInverse std::vector m_Length; double m_BinSize; + std::string m_Particle; std::vector m_LUT; }; } // end namespace Functor diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.h b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h new file mode 100644 index 0000000..5b8d14a --- /dev/null +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.h @@ -0,0 +1,68 @@ +#ifndef __pctFDKDDConeBeamVarianceReconstructionFilter_h +#define __pctFDKDDConeBeamVarianceReconstructionFilter_h + +#include "pctFFTVarianceImageFilter.h" +#include "pctFDKDDConeBeamReconstructionFilter.h" +#include "pctFFTVarianceImageFilter.h" + +#include +#include + +/** \class FDKDDConeBeamReconstructionFilter + * TODO + * + * \author Jannis Dickmann + */ +namespace pct +{ + +template +class ITK_EXPORT FDKDDConeBeamVarianceReconstructionFilter + : public pct::FDKDDConeBeamReconstructionFilter +{ +public: + typedef pct::FDKDDConeBeamReconstructionFilter Baseclass; + + /** Standard class typedefs. */ + typedef FDKDDConeBeamVarianceReconstructionFilter Self; + typedef itk::ImageToImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer 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 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 diff --git a/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx b/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx new file mode 100644 index 0000000..0daa2cc --- /dev/null +++ b/include/pctFDKDDConeBeamVarianceReconstructionFilter.hxx @@ -0,0 +1,33 @@ +namespace pct +{ + +template +FDKDDConeBeamVarianceReconstructionFilter:: + 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 diff --git a/include/pctFDKDDWeightProjectionFilter.h b/include/pctFDKDDWeightProjectionFilter.h index cc0ffb5..6cd8d88 100644 --- a/include/pctFDKDDWeightProjectionFilter.h +++ b/include/pctFDKDDWeightProjectionFilter.h @@ -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() {} @@ -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 diff --git a/include/pctFDKDDWeightProjectionFilter.hxx b/include/pctFDKDDWeightProjectionFilter.hxx index 58a7da3..2c2b470 100644 --- a/include/pctFDKDDWeightProjectionFilter.hxx +++ b/include/pctFDKDDWeightProjectionFilter.hxx @@ -27,6 +27,11 @@ FDKDDWeightProjectionFilter::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]; + } } } diff --git a/include/pctFFTVarianceImageFilter.h b/include/pctFFTVarianceImageFilter.h new file mode 100644 index 0000000..dc6918b --- /dev/null +++ b/include/pctFFTVarianceImageFilter.h @@ -0,0 +1,96 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef pctFFTVarianceImageFilter_h +#define pctFFTVarianceImageFilter_h + +#include "rtkFFTRampImageFilter.h" +#include "rtkFFTProjectionsConvolutionImageFilter.h" +#include "itkHalfHermitianToRealInverseFFTImageFilter.h" +#include "itkImageFileWriter.h" + +#include +#include + +namespace pct +{ + +/** \class FFTVarianceImageFilter + * \brief Implements the variance image filter of the filtered backprojection algorithm. + * + * The filter code is based on FFTConvolutionImageFilter by Gaetan Lehmann + * (see http://hdl.handle.net/10380/3154) + * + * \test rtkvariancefiltertest.cxx + * + * \author Simon Rit + * + * \ingroup ImageToImageFilter + */ + +template +class ITK_EXPORT FFTVarianceImageFilter : public rtk::FFTRampImageFilter +{ +public: + typedef rtk::FFTRampImageFilter Baseclass; + + /** Standard class typedefs. */ + typedef FFTVarianceImageFilter Self; + typedef rtk::FFTRampImageFilter Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Some convenient typedefs. */ + typedef typename Baseclass::InputImageType InputImageType; + typedef typename Baseclass::OutputImageType OutputImageType; + typedef typename Baseclass::FFTPrecisionType FFTPrecisionType; + typedef typename Baseclass::IndexType IndexType; + typedef typename Baseclass::SizeType SizeType; + + typedef typename Baseclass::FFTInputImageType FFTInputImageType; + typedef typename Baseclass::FFTInputImagePointer FFTInputImagePointer; + typedef typename Baseclass::FFTOutputImageType FFTOutputImageType; + typedef typename Baseclass::FFTOutputImagePointer FFTOutputImagePointer; + + /** Standard New method. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkOverrideGetNameOfClassMacro(FFTVarianceImageFilter); + +protected: + FFTVarianceImageFilter(); + ~FFTVarianceImageFilter() {} + + /** Creates and return a pointer to one line of the variance kernel in Fourier space. + * Used in generate data functions. */ + void + UpdateFFTProjectionsConvolutionKernel(const SizeType size) override; + + FFTVarianceImageFilter(const Self &); // purposely not implemented + void + operator=(const Self &); // purposely not implemented +}; // end of class + +} // namespace pct + +#ifndef ITK_MANUAL_INSTANTIATION +# include "pctFFTVarianceImageFilter.hxx" +#endif + +#endif diff --git a/include/pctFFTVarianceImageFilter.hxx b/include/pctFFTVarianceImageFilter.hxx new file mode 100644 index 0000000..19be1b1 --- /dev/null +++ b/include/pctFFTVarianceImageFilter.hxx @@ -0,0 +1,311 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef pctFFTVarianceImageFilter_hxx +#define pctFFTVarianceImageFilter_hxx + +namespace pct +{ + +template +FFTVarianceImageFilter::FFTVarianceImageFilter() +{ + // this->m_HannCutFrequency = 0.; + // this->m_CosineCutFrequency = 0.; + // this->m_HammingFrequency = 0.; + // this->m_HannCutFrequency = 0.; + // this->m_HannCutFrequencyY = 0.; + // this->m_RamLakCutFrequency = 0.; + // this->m_SheppLoganCutFrequency = 0.; +} + +template +void +FFTVarianceImageFilter::UpdateFFTProjectionsConvolutionKernel( + const SizeType s) +{ + // if(this->m_KernelFFT.GetPointer() != ITK_NULLPTR && s == this->m_PreviousKernelUpdateSize) + // { + // return; + // } + // this->m_PreviousKernelUpdateSize = s; + + const int width = s[0]; + const int height = s[1]; + + // Allocate kernel + SizeType size; + size.Fill(1); + size[0] = width; + FFTInputImagePointer kernel = FFTInputImageType::New(); + kernel->SetRegions(size); + kernel->Allocate(); + kernel->FillBuffer(0.); + + // Compute kernel in space domain (see Kak & Slaney, chapter 3 equation 61 + // page 72) although spacing is not squared according to equation 69 page 75 + double spacing = this->GetInput()->GetSpacing()[0]; + IndexType ix, jx; + ix.Fill(0); + jx.Fill(0); + kernel->SetPixel(ix, 1. / (4. * spacing)); + for (ix[0] = 1, jx[0] = size[0] - 1; ix[0] < typename IndexType::IndexValueType(size[0] / 2); ix[0] += 2, jx[0] -= 2) + { + double v = ix[0] * vnl_math::pi; + v = -1. / (v * v * spacing); + kernel->SetPixel(ix, v); + kernel->SetPixel(jx, v); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // typedef itk::ImageRegionIteratorWithIndex InverseIteratorType2; + // InverseIteratorType2 itiK2(kernel, kernel->GetLargestPossibleRegion()); + // std::ofstream out0("kernel_real.csv"); + // itiK2.GoToBegin(); + // for(; !itiK2.IsAtEnd(); ++itiK2) + // { + // out0 << itiK2.Get() << std::endl; + // } + // out0.close(); + // // DEBUG + + // FFT kernel + typedef itk::RealToHalfHermitianForwardFFTImageFilter FFTType; + typename FFTType::Pointer fftK = FFTType::New(); + fftK->SetInput(kernel); + fftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK->Update(); + this->m_KernelFFT = fftK->GetOutput(); + + // Windowing (if enabled) + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType itK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion()); + + unsigned int n = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); + + itK.GoToBegin(); + if (this->GetHannCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetHannCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = 0.5 * (1 + std::cos(itk::Math::pi * i / ncut)); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetCosineCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetCosineCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = std::cos(0.5 * itk::Math::pi * i / ncut); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetHammingFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetHammingFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + double k = 0.54 + 0.46 * (std::cos(itk::Math::pi * i / ncut)); + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else if (this->GetRamLakCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetRamLakCutFrequency())); + for (unsigned int i = 0; i < ncut; i++, ++itK) + { + } + } + else if (this->GetSheppLoganCutFrequency() > 0.) + { + const unsigned int ncut = itk::Math::Round(n * std::min(1.0, this->GetSheppLoganCutFrequency())); + // sinc(0) --> is 1 + ++itK; + for (unsigned int i = 1; i < ncut; i++, ++itK) + { + double x = 0.5 * vnl_math::pi * i / ncut; + double k = std::sin(x) / x; + itK.Set(itK.Get() * TFFTPrecision(k)); + } + } + else + { + itK.GoToReverseBegin(); + ++itK; + } + + for (; !itK.IsAtEnd(); ++itK) + { + itK.Set(itK.Get() * TFFTPrecision(0.)); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out1("kernel_fourier.csv"); + // itK.GoToBegin(); + // for(; !itK.IsAtEnd(); ++itK) + // { + // out1 << itK.Get().real() << " " << itK.Get().imag() << std::endl; + // } + // out1.close(); + // // DEBUG + + // iFFT kernel + typedef itk::HalfHermitianToRealInverseFFTImageFilter IFFTType; + typename IFFTType::Pointer ifftK = IFFTType::New(); + ifftK->SetInput(this->m_KernelFFT); + ifftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + ifftK->Update(); + typename FFTType::InputImageType::Pointer KernelIFFT = ifftK->GetOutput(); + + // calculate ratio g_C/g^2 + typedef itk::ImageRegionIteratorWithIndex InverseIteratorType; + InverseIteratorType itiK(KernelIFFT, KernelIFFT->GetLargestPossibleRegion()); + typename FFTType::InputImageType::PixelType sumgc = 0.; + typename FFTType::InputImageType::PixelType sumg2 = 0.; + typename FFTType::InputImageType::IndexType idxshifted; + const unsigned int maxidx = KernelIFFT->GetLargestPossibleRegion().GetSize()[0]; + for (; !itiK.IsAtEnd(); ++itiK) + { + idxshifted = itiK.GetIndex(); + if (idxshifted[0] == 0) + idxshifted[0] = maxidx - 1; + else + idxshifted[0] -= 1; + + sumgc += itiK.Get() * KernelIFFT->GetPixel(idxshifted); + sumg2 += itiK.Get() * itiK.Get(); + } + const typename FFTType::InputImageType::PixelType ratiogcg2 = sumgc / sumg2; + + // numerical integration to calculate f_interp + const double aprecision = 0.00001; + double finterp = 0.; + for (unsigned int i = 0; i < int(1. / aprecision); i++) + { + const double a = double(i) * aprecision; + finterp += (1 - a) * (1 - a) + 2 * ratiogcg2 * (1 - a) * a + a * a; + } + finterp *= aprecision; + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out9("kernel_window_real.csv"); + // itiK.GoToBegin(); + // for(; !itiK.IsAtEnd(); ++itiK) + // { + // out9 << itiK.Get() << std::endl; + // } + // out9.close(); + // // DEBUG + + // square kernel and multiply with finterp + itiK.GoToBegin(); + for (; !itiK.IsAtEnd(); ++itiK) + { + itiK.Set(itiK.Get() * itiK.Get() * finterp); + } + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out2("kernel_sq_real.csv"); + // itiK.GoToBegin(); + // for(; !itiK.IsAtEnd(); ++itiK) + // { + // out2 << itiK.Get() << std::endl; + // } + // out2.close(); + // // DEBUG + + // FFT kernel + typename FFTType::Pointer fftK2 = FFTType::New(); + fftK2->SetInput(ifftK->GetOutput()); + fftK2->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + fftK2->Update(); + this->m_KernelFFT = fftK2->GetOutput(); + + // ///////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out3("kernel_sq_fourier.csv"); + // IteratorType itK2(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion() ); + // itK2.GoToBegin(); + // for(; !itK2.IsAtEnd(); ++itK2) + // { + // out3 << itK2.Get().real() << " " << itK2.Get().imag() << std::endl; + // } + // out3.close(); + // // DEBUG + + // Replicate and window if required + if (this->GetHannCutFrequencyY() > 0.) + { + std::cerr << "HannY not implemented for variance image filter." << std::endl; + + size.Fill(1); + size[0] = this->m_KernelFFT->GetLargestPossibleRegion().GetSize(0); + size[1] = height; + + const unsigned int ncut = itk::Math::Round((height / 2 + 1) * std::min(1.0, this->GetHannCutFrequencyY())); + + this->m_KernelFFT = FFTOutputImageType::New(); + this->m_KernelFFT->SetRegions(size); + this->m_KernelFFT->Allocate(); + this->m_KernelFFT->FillBuffer(0.); + + IteratorType itTwoDK(this->m_KernelFFT, this->m_KernelFFT->GetLargestPossibleRegion()); + for (unsigned int j = 0; j < ncut; j++) + { + itK.GoToBegin(); + const TFFTPrecision win(0.5 * (1 + std::cos(itk::Math::pi * j / ncut))); + for (unsigned int i = 0; i < size[0]; ++itK, ++itTwoDK, i++) + { + itTwoDK.Set(win * itK.Get()); + } + } + itTwoDK.GoToReverseBegin(); + for (unsigned int j = 1; j < ncut; j++) + { + itK.GoToReverseBegin(); + const TFFTPrecision win(0.5 * (1 + std::cos(itk::Math::pi * j / ncut))); + for (unsigned int i = 0; i < size[0]; --itK, --itTwoDK, i++) + { + itTwoDK.Set(win * itK.Get()); + } + } + } + + this->m_KernelFFT->DisconnectPipeline(); + + // ////////////////////////////////////////////////////////////////////////////////////////////// + // // DEBUG + // std::ofstream out4("finterp.csv"); + // out4 << finterp << std::endl; + // out4 << ratiogcg2 << std::endl; + // out4.close(); + // std::cerr << std::endl << "finterp = " << finterp << std::endl; + // std::cerr << "ratiogcg2 = " << ratiogcg2 << std::endl; + // std::cerr << "Aborting here since filters were alrady written out." << std::endl; + // exit(9); + // // DEBUG +} + +} // namespace pct +#endif diff --git a/include/pctProtonPairsToDistanceDrivenProjection.h b/include/pctProtonPairsToDistanceDrivenProjection.h index af7fe00..54f908c 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.h +++ b/include/pctProtonPairsToDistanceDrivenProjection.h @@ -29,6 +29,9 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection using CountImageType = itk::Image; using CountImagePointer = CountImageType::Pointer; + typedef itk::Image VarianceImageType; + typedef VarianceImageType::Pointer VarianceImagePointer; + using AngleImageType = itk::Image; using AngleImagePointer = AngleImageType::Pointer; @@ -78,6 +81,9 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection /** Get/Set the count of proton pairs per pixel. */ itkGetMacro(Count, CountImagePointer); + /** Get/Set the variance of proton pairs per pixel. */ + itkGetMacro(Variance, VarianceImagePointer); + /** Get/Set the angle of proton pairs per pixel. */ itkGetMacro(Angle, AngleImagePointer); @@ -110,6 +116,16 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection itkGetConstMacro(ComputeNoise, bool); itkBooleanMacro(ComputeNoise); + /** Do the variance projections + ** Default is off. */ + itkSetMacro(ComputeVariance, bool); + itkGetConstMacro(ComputeVariance, bool); + itkBooleanMacro(ComputeVariance); + + /** Get/Set the particle type. */ + itkGetMacro(Particle, std::string); + itkSetMacro(Particle, std::string); + protected: ProtonPairsToDistanceDrivenProjection(); virtual ~ProtonPairsToDistanceDrivenProjection() {} @@ -150,6 +166,12 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection CountImagePointer m_Count; std::vector m_Counts; + /** Calculate variance in each thread */ + VarianceImagePointer m_Variance; + std::vector m_Variances; + VarianceImagePointer m_VarianceSq; + std::vector m_VariancesSq; + AngleImagePointer m_Angle; std::vector m_Angles; std::vector> m_AnglesVectors; @@ -181,7 +203,11 @@ class ITK_TEMPLATE_EXPORT ProtonPairsToDistanceDrivenProjection ProtonPairsImageType::Pointer m_ProtonPairs; bool m_Robust; bool m_ComputeScattering; + bool m_ComputeVariance; bool m_ComputeNoise; + + /** The particle type, enum comes from gengetopt header. */ + std::string m_Particle; }; } // end namespace pct diff --git a/include/pctProtonPairsToDistanceDrivenProjection.hxx b/include/pctProtonPairsToDistanceDrivenProjection.hxx index 37db9f6..c9c3f41 100644 --- a/include/pctProtonPairsToDistanceDrivenProjection.hxx +++ b/include/pctProtonPairsToDistanceDrivenProjection.hxx @@ -15,6 +15,7 @@ ProtonPairsToDistanceDrivenProjection::ProtonPairsToD : m_Robust(false) , m_ComputeScattering(false) , m_ComputeNoise(false) + , m_ComputeVariance(false) { this->DynamicMultiThreadingOff(); this->SetNumberOfWorkUnits(itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads()); @@ -26,6 +27,13 @@ ProtonPairsToDistanceDrivenProjection::BeforeThreaded { m_Outputs.resize(this->GetNumberOfWorkUnits()); m_Counts.resize(this->GetNumberOfWorkUnits()); + + if (m_ComputeVariance) + { + m_Variances.resize(this->GetNumberOfWorkUnits()); + m_VariancesSq.resize(this->GetNumberOfWorkUnits()); + } + if (m_ComputeScattering) { m_Angles.resize(this->GetNumberOfWorkUnits()); @@ -40,8 +48,9 @@ ProtonPairsToDistanceDrivenProjection::BeforeThreaded if (m_QuadricOut.GetPointer() == NULL) m_QuadricOut = m_QuadricIn; + const int zsq = (m_Particle == "alpha") ? 4 : 1; // helium energy is 800 MeV instead of 200 MeV for protons m_ConvFunc = new Functor::IntegratedBetheBlochProtonStoppingPowerInverse( - m_IonizationPotential, 600. * CLHEP::MeV, 0.1 * CLHEP::keV); + m_IonizationPotential, zsq * 600. * CLHEP::MeV, 0.1 * CLHEP::keV, m_Particle); // Read pairs using ReaderType = itk::ImageFileReader; @@ -57,12 +66,16 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera const OutputImageRegionType & itkNotUsed(outputRegionForThread), rtk::ThreadIdType threadId) { + if (m_Particle != "proton" && m_MostLikelyPathType != "schulte") + itkGenericExceptionMacro("Only Schulte's MLP can be used with particles which are not proton"); + // Create MLP depending on type pct::MostLikelyPathFunction::Pointer mlp; if (m_MostLikelyPathType == "polynomial") mlp = pct::ThirdOrderPolynomialMLPFunction::New(); else if (m_MostLikelyPathType == "krah") { + mlp = pct::ThirdOrderPolynomialMLPFunction::New(); pct::PolynomialMLPFunction::Pointer mlp_poly; mlp_poly = pct::PolynomialMLPFunction::New(); mlp_poly->SetPolynomialDegree(m_MostLikelyPathPolynomialDegree); @@ -74,7 +87,9 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera } else if (m_MostLikelyPathType == "schulte") { - mlp = pct::SchulteMLPFunction::New(); + pct::SchulteMLPFunction::Pointer mlpSchulte = pct::SchulteMLPFunction::New(); + mlpSchulte->SetParticle(m_Particle); + mlp = mlpSchulte; } else { @@ -92,6 +107,19 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera m_Counts[threadId]->Allocate(); m_Counts[threadId]->FillBuffer(0); + if (m_ComputeVariance) + { + m_Variances[threadId] = VarianceImageType::New(); + m_Variances[threadId]->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_Variances[threadId]->Allocate(); + m_Variances[threadId]->FillBuffer(0); + + m_VariancesSq[threadId] = VarianceImageType::New(); + m_VariancesSq[threadId]->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_VariancesSq[threadId]->Allocate(); + m_VariancesSq[threadId]->FillBuffer(0); + } + if (m_ComputeScattering && (!m_Robust || threadId == 0)) // Note NK: is this condition correct? Should it not be !(m_Robust || threadId==0) ? { @@ -118,6 +146,11 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera { m_Outputs[0] = this->GetOutput(); m_Count = m_Counts[0]; + if (m_ComputeVariance) + { + m_Variance = m_Variances[0]; + m_VarianceSq = m_VariancesSq[0]; + } if (m_ComputeScattering) { m_Angle = m_Angles[0]; @@ -152,7 +185,15 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera typename OutputImageType::PixelType * imgData = m_Outputs[threadId]->GetBufferPointer(); typename OutputImageType::PixelType * imgSquaredData = NULL; unsigned int * imgCountData = m_Counts[threadId]->GetBufferPointer(); - float * imgAngleData = NULL, *imgAngleSqData = NULL; + + float *imgVarianceData = NULL, *imgVarianceSqData = NULL; + if (m_ComputeVariance) + { + imgVarianceData = m_Variances[threadId]->GetBufferPointer(); + imgVarianceSqData = m_VariancesSq[threadId]->GetBufferPointer(); + } + + float *imgAngleData = NULL, *imgAngleSqData = NULL; if (m_ComputeScattering && !m_Robust) { imgAngleData = m_Angles[threadId]->GetBufferPointer(); @@ -431,6 +472,11 @@ ProtonPairsToDistanceDrivenProjection::ThreadedGenera { imgSquaredData[idx] += value * value; } + if (m_ComputeVariance) + { + imgVarianceData[idx] += value; + imgVarianceSqData[idx] += value * value; + } imgCountData[idx]++; if (m_ComputeScattering) { @@ -545,6 +591,61 @@ ProtonPairsToDistanceDrivenProjection::AfterThreadedG } } + if (m_ComputeVariance) + { + typedef itk::ImageRegionIterator ImageVarianceIteratorType; + ImageVarianceIteratorType itVarianceOut(m_Variances[0], m_Outputs[0]->GetLargestPossibleRegion()); + + typedef itk::ImageRegionIterator ImageVarianceSqIteratorType; + ImageVarianceSqIteratorType itVarianceSqOut(m_VariancesSq[0], m_Outputs[0]->GetLargestPossibleRegion()); + + for (unsigned int i = 1; i < this->GetNumberOfWorkUnits(); i++) + { + if (m_Outputs[i].GetPointer() == NULL) + continue; + ImageVarianceIteratorType itVarianceOutThread(m_Variances[i], m_Outputs[i]->GetLargestPossibleRegion()); + ImageVarianceSqIteratorType itVarianceSqOutThread(m_VariancesSq[i], m_Outputs[i]->GetLargestPossibleRegion()); + + while (!itVarianceOut.IsAtEnd()) + { + itVarianceOut.Set(itVarianceOut.Get() + itVarianceOutThread.Get()); + ++itVarianceOutThread; + ++itVarianceOut; + + itVarianceSqOut.Set(itVarianceSqOut.Get() + itVarianceSqOutThread.Get()); + ++itVarianceSqOutThread; + ++itVarianceSqOut; + } + itVarianceOut.GoToBegin(); + itVarianceSqOut.GoToBegin(); + } + + // Set variance image information + m_Variance->SetSpacing(this->GetOutput()->GetSpacing()); + m_Variance->SetOrigin(this->GetOutput()->GetOrigin()); + + itCOut.GoToBegin(); + while (!itCOut.IsAtEnd()) + { + if (itCOut.Get() >= 2) + { + const float vsum = itVarianceOut.Get(); + const float vsqsum = itVarianceSqOut.Get(); + const float n = itCOut.Get(); + const float variance_wepl = (vsqsum - vsum * vsum / n) / (n - 1); + const float variance_proj = variance_wepl / n; + itVarianceOut.Set(variance_proj); + } + else + { + itVarianceOut.Set(0.); + } + ++itCOut; + ++itVarianceOut; + ++itVarianceSqOut; + } + } + if (m_ComputeScattering) { using ImageAngleIteratorType = itk::ImageRegionIterator; @@ -625,6 +726,8 @@ ProtonPairsToDistanceDrivenProjection::AfterThreadedG m_Outputs.resize(0); m_SquaredOutputs.resize(0); m_Counts.resize(0); + m_Variances.resize(0); + m_VariancesSq.resize(0); m_Angles.resize(0); m_AnglesSq.resize(0); m_AnglesVectors.resize(0); diff --git a/include/pctSchulteMLPFunction.h b/include/pctSchulteMLPFunction.h index f46a873..c7acdce 100644 --- a/include/pctSchulteMLPFunction.h +++ b/include/pctSchulteMLPFunction.h @@ -37,6 +37,12 @@ static const double a2 = -9.986645e-08 * aunit / (CLHEP::cm2); static const double a3 = 2.026409e-08 * aunit / (CLHEP::cm3); static const double a4 = -1.420501e-09 * aunit / (CLHEP::cm3 * CLHEP::cm); static const double a5 = 3.899100e-11 * aunit / (CLHEP::cm3 * CLHEP::cm2); +static const double a0He = 4.665553448715765888e-07 * aunit; +static const double a1He = 2.817219469097221029e-08 * aunit / (CLHEP::cm); +static const double a2He = -3.534701631484195519e-09 * aunit / (CLHEP::cm2); +static const double a3He = 7.963860554062791474e-10 * aunit / (CLHEP::cm3); +static const double a4He = -5.635546522420130709e-11 * aunit / (CLHEP::cm3 * CLHEP::cm); +static const double a5He = 1.635484142530505741e-12 * aunit / (CLHEP::cm3 * CLHEP::cm2); #elif defined(NICO_COEFF) static const double a0 = 3.599e-06 * aunit; static const double a1 = 7.303e-08 * aunit / (CLHEP::cm); @@ -69,6 +75,11 @@ class ConstantPartOfIntegrals double v = 1. + 0.038 * std::log((diffU)*invX0); return c * v * v; } + static double + GetValueHe(const double ux, const double uy) + { + return GetValue(ux, uy); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 8 @@ -88,6 +99,20 @@ class IntegralForSigmaSqTheta static const double a5Theta = a5 / 6.; return u * (a0Theta + u * (a1Theta + u * (a2Theta + u * (a3Theta + u * (a4Theta + u * a5Theta))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0Theta = a0He; + static const double a1Theta = a1He / 2.; + static const double a2Theta = a2He / 3.; + static const double a3Theta = a3He / 4.; + static const double a4Theta = a4He / 5.; + static const double a5Theta = a5He / 6.; + return u * (a0Theta + u * (a1Theta + u * (a2Theta + u * (a3Theta + u * (a4Theta + u * a5Theta))))); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 9 @@ -107,6 +132,20 @@ class IntegralForSigmaSqTTheta static const double a5TTheta = a5 / 7.; return u * u * (a0TTheta + u * (a1TTheta + u * (a2TTheta + u * (a3TTheta + u * (a4TTheta + u * a5TTheta))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0TTheta = a0He / 2.; + static const double a1TTheta = a1He / 3.; + static const double a2TTheta = a2He / 4.; + static const double a3TTheta = a3He / 5.; + static const double a4TTheta = a4He / 6.; + static const double a5TTheta = a5He / 7.; + return u * u * (a0TTheta + u * (a1TTheta + u * (a2TTheta + u * (a3TTheta + u * (a4TTheta + u * a5TTheta))))); + } }; // [Schulte, Med Phys, 2008], integral for sigma in equation 7 @@ -126,6 +165,20 @@ class IntegralForSigmaSqT static const double a5T = a5 / 8.; return u * u * u * (a0T + u * (a1T + u * (a2T + u * (a3T + u * (a4T + u * a5T))))); } + + static double + GetValueHe(const double u) + { + // Multiplied with polynomial integral + // Table 2, (a) in [Williams, 2004] as well as numbers in [Li, 2006] + static const double a0T = a0He / 3.; + static const double a1T = a1He / 4.; + static const double a2T = a2He / 5.; + static const double a3T = a3He / 6.; + static const double a4T = a4He / 7.; + static const double a5T = a5He / 8.; + return u * u * u * (a0T + u * (a1T + u * (a2T + u * (a3T + u * (a4T + u * a5T))))); + } }; } // end namespace SchulteMLP @@ -176,6 +229,18 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction void EvaluateError(const double u1, itk::Matrix & error); + /** Get/Set the particle type. */ + void + SetParticle(const std::string particle) + { + m_Particle = particle; + } + std::string + GetParticle() + { + return m_Particle; + } + #ifdef MLP_TIMING /** Print timing information */ virtual void @@ -238,6 +303,8 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction double m_IntForSigmaSqTTheta2; double m_IntForSigmaSqT2; + std::string m_Particle{ "proton" }; + #ifdef MLP_TIMING itk::TimeProbe m_EvaluateProbe1; itk::TimeProbe m_EvaluateProbe2; @@ -246,5 +313,4 @@ class PCT_EXPORT SchulteMLPFunction : public MostLikelyPathFunction } // end namespace pct - #endif diff --git a/pctpairprotonsLMU_IMPCT.cxx b/pctpairprotonsLMU_IMPCT.cxx index 0a21c6e..995632d 100644 --- a/pctpairprotonsLMU_IMPCT.cxx +++ b/pctpairprotonsLMU_IMPCT.cxx @@ -160,6 +160,33 @@ main(int argc, char * argv[]) treeIn->GetEntry(iIn); treeOut->GetEntry(iOut); + // Fluence modulation for ideal data, ROI always 0 + // Part for proton track segment - corcular ROI intersection + // http://mathworld.wolfram.com/Circle-LineIntersection.html + bool interceptionFlag = true; + if (args_info.fmpct_flag) + { + double radius = args_info.roiR_arg; // ROI radius in mm + if (radius != 0) + { + std::cout << "Warning: ROI=0 for ideal data" << std::endl; + } + // double dx = pd.position2[2]-pd.position1[2]; + // double dy = pd.position2[0]-pd.position1[0]; + // double dr_sq = (dx*dx) + (dy*dy); + // double D = (pd.position1[2]*pd.position2[0]) - (pd.position2[2]*pd.position1[0]); + // double Delta = (radius*radius*dr_sq) - (D*D); + double randomNumber = ((double)rand() / (RAND_MAX)); + // std::cout << "RAND_MAX = " << RAND_MAX << std::endl; + // std::cout << "randomNumber = " << randomNumber << std::endl; + // no interception, then discard proton with randomNumber > modF_arg + // if (Delta<0) { + if (randomNumber > args_info.modF_arg) + interceptionFlag = false; + //} + // std::cout << "Delta, randomNumber, interceptionFlag = " << Delta << " " << randomNumber << " " << + // interceptionFlag <=args_info.minRun_arg && args_info.runID_arg=args_info.minRun_arg && args_info.runID_arg position0; itk::Vector position1; itk::Vector position2; @@ -66,7 +66,7 @@ BranchParticleToPhaseSpace(struct ParticleInfo & piInput, struct ParticleData & SetTreeBranch(tree, "projection_angle", &piInput.runID); ; SetTreeBranch(tree, "calculated_WEPL", &pdInput.wepl); - + SetTreeBranch(tree, "timestamp", &pdInput.time); // WARNING: X and Z are purposely swap... SetTreeBranch(tree, "v_hit0", pdInput.position0.GetDataPointer() + 1); @@ -207,7 +207,7 @@ main(int argc, char * argv[]) pdIn.wepl = 0.; pdOut.wepl = pd.wepl; pdIn.time = 0.; - pdOut.time = 0.; + pdOut.time = pd.time; #if 0 @@ -244,9 +244,9 @@ main(int argc, char * argv[]) if (pairs[i].size()) { std::ostringstream os; - os << itksys::SystemTools::GetFilenameWithoutLastExtension(args_info.output_arg) << std::setw(4) - << std::setfill('0') << i << itksys::SystemTools::GetFilenameLastExtension(args_info.output_arg); - std::cout << "Writing into file:" << os.str() << std::endl; + os << itksys::SystemTools::GetFilenamePath(args_info.output_arg) << "/" + << itksys::SystemTools::GetFilenameWithoutLastExtension(args_info.output_arg) std::cout + << "Writing into file:" << os.str() << std::endl; WritePairs(pairs[i], os.str()); } } diff --git a/src/pctSchulteMLPFunction.cxx b/src/pctSchulteMLPFunction.cxx index a9d1146..8e58c67 100644 --- a/src/pctSchulteMLPFunction.cxx +++ b/src/pctSchulteMLPFunction.cxx @@ -105,9 +105,22 @@ SchulteMLPFunction ::Init(const VectorType posIn, m_uOrigin = posIn[2]; m_u2 = posOut[2] - m_uOrigin; - m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(m_u2); - m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(m_u2); - m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(m_u2); + if (m_Particle == "proton") + { + m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(m_u2); + m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(m_u2); + m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(m_u2); + } + else if (m_Particle == "alpha") + { + m_IntForSigmaSqTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValueHe(m_u2); + m_IntForSigmaSqTTheta2 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValueHe(m_u2); + m_IntForSigmaSqT2 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValueHe(m_u2); + } + else + { + throw std::invalid_argument("Invalid particle in Schulte MLP function."); + } // Parameters vectors m_x0[0] = posIn[0]; @@ -141,27 +154,51 @@ SchulteMLPFunction ::Evaluate(const double u, double & x, double & y, double & d InverseMatrix(R1_Inv); // Constants used in both integrals - const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(u1); - const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(u1); - const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(u1); - - // Construct Sigma1 (equations 6-9) - m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; - m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; - m_Sigma1(1, 0) = m_Sigma1(0, 1); - m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; - m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(m_u0, u1); - - double sigma1 = std::sqrt(m_Sigma2(1, 1)); - - // Construct Sigma2 (equations 15-18) - m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; - m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; - m_Sigma2(1, 0) = m_Sigma2(0, 1); - m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; - m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1, m_u2); + if (m_Particle == "proton") + { + const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValue(u1); + const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValue(u1); + const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValue(u1); + + // Construct Sigma1 (equations 6-9) + m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; + m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; + m_Sigma1(1, 0) = m_Sigma1(0, 1); + m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; + m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(m_u0, u1); + + double sigma1 = std::sqrt(m_Sigma2(1, 1)); + + // Construct Sigma2 (equations 15-18) + m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; + m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; + m_Sigma2(1, 0) = m_Sigma2(0, 1); + m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; + m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1, m_u2); + } + else if (m_Particle == "alpha") + { + const double intForSigmaSqTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTheta ::GetValueHe(u1); + const double intForSigmaSqTTheta1 = Functor::SchulteMLP::IntegralForSigmaSqTTheta::GetValueHe(u1); + const double intForSigmaSqT1 = Functor::SchulteMLP::IntegralForSigmaSqT ::GetValueHe(u1); + + // Construct Sigma1 (equations 6-9) + m_Sigma1(1, 1) = intForSigmaSqTheta1 /* - m_IntForSigmaSqTheta0*/; + m_Sigma1(0, 1) = u1 * m_Sigma1(1, 1) - intForSigmaSqTTheta1 /* + m_IntForSigmaSqTTheta0*/; + m_Sigma1(1, 0) = m_Sigma1(0, 1); + m_Sigma1(0, 0) = u1 * (2 * m_Sigma1(0, 1) - u1 * m_Sigma1(1, 1)) + intForSigmaSqT1 /* - m_IntForSigmaSqT0*/; + m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValueHe(m_u0, u1); + + // Construct Sigma2 (equations 15-18) + m_Sigma2(1, 1) = m_IntForSigmaSqTheta2 - intForSigmaSqTheta1; + m_Sigma2(0, 1) = m_u2 * m_Sigma2(1, 1) - m_IntForSigmaSqTTheta2 + intForSigmaSqTTheta1; + m_Sigma2(1, 0) = m_Sigma2(0, 1); + m_Sigma2(0, 0) = m_u2 * (2 * m_Sigma2(0, 1) - m_u2 * m_Sigma2(1, 1)) + m_IntForSigmaSqT2 - intForSigmaSqT1; + m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValueHe(u1, m_u2); + } #ifdef MLP_TIMING + m_EvaluateProbe1.Stop(); m_EvaluateProbe2.Start(); #endif