From 111db1c2131fcaed78b78d081878c4236452a193 Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Fri, 18 Oct 2019 16:51:56 +0200 Subject: [PATCH 1/8] Added option to use Collins-Fekete binning --- pctProtonPairsToBackProjection.h | 10 ++++++- pctProtonPairsToBackProjection.txx | 42 +++++++++++++++++++++++++++--- pctbackprojectionbinning.cxx | 1 + pctbackprojectionbinning.ggo | 2 +- 4 files changed, 49 insertions(+), 6 deletions(-) diff --git a/pctProtonPairsToBackProjection.h b/pctProtonPairsToBackProjection.h index d4d0f4b..457b294 100644 --- a/pctProtonPairsToBackProjection.h +++ b/pctProtonPairsToBackProjection.h @@ -29,7 +29,7 @@ class ITK_EXPORT ProtonPairsToBackProjection : typedef itk::Image ProtonPairsImageType; typedef ProtonPairsImageType::Pointer ProtonPairsImagePointer; - typedef typename itk::Image< unsigned int, + typedef typename itk::Image< double, TInputImage::ImageDimension> CountImageType; typedef typename CountImageType::Pointer CountImagePointer; @@ -90,6 +90,11 @@ class ITK_EXPORT ProtonPairsToBackProjection : itkSetMacro(DisableRotation, bool); itkBooleanMacro(DisableRotation); + /** Get / Set the boolean to use Collins-Fekete binning. Default is off. */ + itkGetMacro(WeightsCF, bool); + itkSetMacro(WeightsCF, bool); + itkBooleanMacro(WeightsCF); + protected: ProtonPairsToBackProjection(); virtual ~ProtonPairsToBackProjection() {} @@ -133,6 +138,9 @@ class ITK_EXPORT ProtonPairsToBackProjection : bool m_DisableRotation = false; std::mutex m_Mutex; + + /** Use binning as in Collins-Fekete */ + bool m_WeightsCF; }; } // end namespace pct diff --git a/pctProtonPairsToBackProjection.txx b/pctProtonPairsToBackProjection.txx index 0b15ad8..37626f9 100644 --- a/pctProtonPairsToBackProjection.txx +++ b/pctProtonPairsToBackProjection.txx @@ -96,7 +96,7 @@ ProtonPairsToBackProjection const typename OutputImageType::SpacingType imgSpacing = this->GetInput()->GetSpacing(); typename OutputImageType::PixelType *imgData = this->GetOutput()->GetBufferPointer(); - unsigned int *imgCountData = m_Counts->GetBufferPointer(); + double *imgCountData = m_Counts->GetBufferPointer(); itk::Vector imgSpacingInv; double minSpacing = imgSpacing[0]; for(unsigned int i=0; i<3; i++) @@ -179,6 +179,14 @@ ProtonPairsToBackProjection value = m_ConvFunc->GetValue(eOut, eIn); // convert to WEPL ++it; + VectorType nucInfo(0.); + + if(it.GetIndex()[0] != 0) + { + nucInfo = it.Get(); + ++it; + } + // Move straight to entrance and exit shapes VectorType pSIn = pIn; VectorType pSOut = pOut; @@ -211,6 +219,9 @@ ProtonPairsToBackProjection // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); + int totalLength = 0.; + std::vector offsets; + for(unsigned int k=0; k idx[1]>=0 && idx[1]<(int)imgSize[1] && idx[2]>=0 && idx[2]<(int)imgSize[2]) { - typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); + if(m_WeightsCF) + { + idx[2] = 0; + offsets.push_back(this->GetOutput()->ComputeOffset(idx)); + totalLength++; + } + else + { + typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); + m_Mutex.lock(); + imgData[ offset ] += value; + imgCountData[ offset ]++; + m_Mutex.unlock(); + } + } + } + if(m_WeightsCF) + { + std::vector channels(offsets); + std::sort( channels.begin(), channels.end() ); + channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() ); + for(unsigned int k=0; kSetMostLikelyPathType( args_info.mlptype_arg ); projection->SetIonizationPotential( args_info.ionpot_arg * CLHEP::eV ); projection->SetDisableRotation( args_info.norotation_flag ); + projection->SetWeightsCF( args_info.weightsCF_flag ); // Geometry if(args_info.verbose_flag) diff --git a/pctbackprojectionbinning.ggo b/pctbackprojectionbinning.ggo index 2597056..3eeea30 100644 --- a/pctbackprojectionbinning.ggo +++ b/pctbackprojectionbinning.ggo @@ -16,6 +16,7 @@ option "geometry" - "XML geometry file name" option "bpVal" - "Input backprojection image values" string no option "bpCount" - "Input backprojection image counts" string no option "norotation" - "Bin in parallel coordinate system" flag off +option "weightsCF" - "Use Collins-Fekete binning" flag off section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no @@ -23,4 +24,3 @@ option "dimension" - "Dimension" int multiple no default="256 option "spacing" - "Spacing" double multiple no default="1" option "direction" - "Direction" double multiple no option "like" - "Copy information from this image (origin, dimension, spacing, direction)" string no - From 81f901f53ec1d5605113160c89b2c0d39998c72b Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Tue, 22 Oct 2019 11:32:34 +0200 Subject: [PATCH 2/8] Added option to used Collins-Fekete weights in distance-driven binning --- pctProtonPairsToDistanceDrivenProjection.h | 10 +++++- pctProtonPairsToDistanceDrivenProjection.txx | 32 ++++++++++++++++++-- pctbinning.cxx | 1 + pctbinning.ggo | 2 +- 4 files changed, 40 insertions(+), 5 deletions(-) diff --git a/pctProtonPairsToDistanceDrivenProjection.h b/pctProtonPairsToDistanceDrivenProjection.h index 18cc725..5318617 100644 --- a/pctProtonPairsToDistanceDrivenProjection.h +++ b/pctProtonPairsToDistanceDrivenProjection.h @@ -26,7 +26,7 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : typedef itk::Image ProtonPairsImageType; typedef ProtonPairsImageType::Pointer ProtonPairsImagePointer; - typedef itk::Image CountImageType; + typedef itk::Image CountImageType; typedef CountImageType::Pointer CountImagePointer; typedef itk::Image AngleImageType; @@ -84,6 +84,11 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : itkGetConstMacro(ComputeScattering, bool); itkBooleanMacro(ComputeScattering); + /** Get / Set the boolean to use Collins-Fekete binning. Default is off. */ + itkGetMacro(WeightsCF, bool); + itkSetMacro(WeightsCF, bool); + itkBooleanMacro(WeightsCF); + protected: ProtonPairsToDistanceDrivenProjection(); virtual ~ProtonPairsToDistanceDrivenProjection() {} @@ -135,6 +140,9 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : ProtonPairsImageType::Pointer m_ProtonPairs; bool m_Robust; bool m_ComputeScattering; + + /** Use binning as in Collins-Fekete */ + bool m_WeightsCF; }; } // end namespace pct diff --git a/pctProtonPairsToDistanceDrivenProjection.txx b/pctProtonPairsToDistanceDrivenProjection.txx index 79e6989..2e1fb6d 100644 --- a/pctProtonPairsToDistanceDrivenProjection.txx +++ b/pctProtonPairsToDistanceDrivenProjection.txx @@ -107,7 +107,7 @@ ProtonPairsToDistanceDrivenProjection const unsigned long npixelsPerSlice = imgSize[0] * imgSize[1]; typename OutputImageType::PixelType *imgData = m_Outputs[threadId]->GetBufferPointer(); - unsigned int *imgCountData = m_Counts[threadId]->GetBufferPointer(); + double *imgCountData = m_Counts[threadId]->GetBufferPointer(); float *imgAngleData = NULL, *imgAngleSqData = NULL; if(m_ComputeScattering && !m_Robust) { @@ -230,6 +230,9 @@ ProtonPairsToDistanceDrivenProjection // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); + int totalLength = 0.; + std::vector offsets; + for(unsigned int k=0; k if(i>=0 && i<(int)imgSize[0] && j>=0 && j<(int)imgSize[1]) { + //const unsigned long idx; const unsigned long idx = i+j*imgSize[0]+k*npixelsPerSlice; - imgData[ idx ] += value; - imgCountData[ idx ]++; + if(m_WeightsCF) + { + offsets.push_back(i+j*imgSize[0]); + totalLength++; + } + else + { + imgData[ idx ] += value; + imgCountData[ idx ]++; + } + if(m_ComputeScattering) { if(m_Robust) @@ -282,6 +295,19 @@ ProtonPairsToDistanceDrivenProjection } } } + if(m_WeightsCF) + { + std::vector channels(offsets); + std::sort( channels.begin(), channels.end() ); + channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() ); + for(unsigned int k=0; kSetIonizationPotential( args_info.ionpot_arg * CLHEP::eV ); projection->SetRobust( args_info.robust_flag ); projection->SetComputeScattering( args_info.scatwepl_given ); + projection->SetWeightsCF( args_info.weightsCF_flag ); if(args_info.quadricIn_given) { diff --git a/pctbinning.ggo b/pctbinning.ggo index 25422e0..dfa31b8 100644 --- a/pctbinning.ggo +++ b/pctbinning.ggo @@ -16,6 +16,7 @@ option "quadricOut" - "Parameters of the exit quadric support function, see http option "mlptype" - "Type of most likely path (schulte or polynomial)" string no default="schulte" 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 "weightsCF" - "Use Collins-Fekete binning" flag off section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no @@ -23,4 +24,3 @@ option "dimension" - "Dimension" int multiple no default="256 option "spacing" - "Spacing" double multiple no default="1" option "direction" - "Direction" double multiple no option "like" - "Copy information from this image (origin, dimension, spacing, direction)" string no - From a95e491e953fa0ece0a2cea6c3d735a91e5a107b Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Fri, 8 Nov 2019 16:26:50 +0100 Subject: [PATCH 3/8] Fixed error with Collins-Fekete binning in distance-driven binning --- pctProtonPairsToDistanceDrivenProjection.txx | 25 ++++++++++---------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/pctProtonPairsToDistanceDrivenProjection.txx b/pctProtonPairsToDistanceDrivenProjection.txx index 2e1fb6d..478e2dd 100644 --- a/pctProtonPairsToDistanceDrivenProjection.txx +++ b/pctProtonPairsToDistanceDrivenProjection.txx @@ -295,19 +295,18 @@ ProtonPairsToDistanceDrivenProjection } } } - if(m_WeightsCF) - { - std::vector channels(offsets); - std::sort( channels.begin(), channels.end() ); - channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() ); - for(unsigned int k=0; k channels(offsets); + std::sort( channels.begin(), channels.end() ); + channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() ); + for(unsigned int k=0; k Date: Wed, 4 Dec 2019 11:13:33 +0100 Subject: [PATCH 4/8] Fixed issue with nonuclear option --- pctpairprotons.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pctpairprotons.cxx b/pctpairprotons.cxx index e115aba..0a49ff0 100644 --- a/pctpairprotons.cxx +++ b/pctpairprotons.cxx @@ -120,9 +120,9 @@ void WritePairs(const std::vector< std::pair > &pairs ++it; if(size[0] == 6) { - eet[0] = particlesInfo[i].creatorProcess; - eet[1] = particlesInfo[i].nuclearProcess; - eet[2] = particlesInfo[i].order; + nuclearinfo[0] = particlesInfo[i].creatorProcess; + nuclearinfo[1] = particlesInfo[i].nuclearProcess; + nuclearinfo[2] = particlesInfo[i].order; it.Set( nuclearinfo ); ++it; From 06266dc3f5073eed6fd1eab790b70dbd6c23b3f5 Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Wed, 18 Dec 2019 15:36:20 +0100 Subject: [PATCH 5/8] Use constant total length rather than length inside image grid for Collins-Fekete binning --- pctProtonPairsToDistanceDrivenProjection.txx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pctProtonPairsToDistanceDrivenProjection.txx b/pctProtonPairsToDistanceDrivenProjection.txx index 478e2dd..0f17e75 100644 --- a/pctProtonPairsToDistanceDrivenProjection.txx +++ b/pctProtonPairsToDistanceDrivenProjection.txx @@ -230,7 +230,6 @@ ProtonPairsToDistanceDrivenProjection // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); - int totalLength = 0.; std::vector offsets; for(unsigned int k=0; k if(m_WeightsCF) { offsets.push_back(i+j*imgSize[0]); - totalLength++; } else { @@ -303,7 +301,7 @@ ProtonPairsToDistanceDrivenProjection channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() ); for(unsigned int k=0; k Date: Tue, 19 May 2020 00:45:51 +0200 Subject: [PATCH 6/8] Added option to produce MLP uncertainty map including tracker information --- pctProtonPairsToBackProjection.h | 8 ++ pctProtonPairsToBackProjection.txx | 51 ++++++-- pctProtonPairsToDistanceDrivenProjection.h | 13 ++ pctProtonPairsToDistanceDrivenProjection.txx | 90 +++++++++++-- pctSchulteMLPFunction.h | 14 ++ pctSchulteMLPFunction.txx | 129 ++++++++++++++++--- pctbackprojectionbinning.cxx | 1 + pctbackprojectionbinning.ggo | 1 + pctbinning.cxx | 11 ++ pctbinning.ggo | 1 + 10 files changed, 284 insertions(+), 35 deletions(-) diff --git a/pctProtonPairsToBackProjection.h b/pctProtonPairsToBackProjection.h index 457b294..81fa3c9 100644 --- a/pctProtonPairsToBackProjection.h +++ b/pctProtonPairsToBackProjection.h @@ -95,6 +95,11 @@ class ITK_EXPORT ProtonPairsToBackProjection : itkSetMacro(WeightsCF, bool); itkBooleanMacro(WeightsCF); + /** Get / Set the boolean to output MLP uncertainty map. Default is off. */ + itkGetMacro(SigmaMap, bool); + itkSetMacro(SigmaMap, bool); + itkBooleanMacro(SigmaMap); + protected: ProtonPairsToBackProjection(); virtual ~ProtonPairsToBackProjection() {} @@ -141,6 +146,9 @@ class ITK_EXPORT ProtonPairsToBackProjection : /** Use binning as in Collins-Fekete */ bool m_WeightsCF; + + /** Output sigma map */ + bool m_SigmaMap; }; } // end namespace pct diff --git a/pctProtonPairsToBackProjection.txx b/pctProtonPairsToBackProjection.txx index 37626f9..b1c4f5b 100644 --- a/pctProtonPairsToBackProjection.txx +++ b/pctProtonPairsToBackProjection.txx @@ -62,15 +62,16 @@ ProtonPairsToBackProjection [this, iProj](const ProtonPairsImageType::RegionType & outputRegionForThread) { // Create MLP depending on type - pct::MostLikelyPathFunction::Pointer mlp; - if(m_MostLikelyPathType == "polynomial") - mlp = pct::ThirdOrderPolynomialMLPFunction::New(); - else if (m_MostLikelyPathType == "schulte") - mlp = pct::SchulteMLPFunction::New(); - else - { - itkGenericExceptionMacro("MLP must either be schulte or polynomial, not [" << m_MostLikelyPathType << ']'); - } + //pct::MostLikelyPathFunction::Pointer mlp; + pct::SchulteMLPFunction::Pointer mlp = pct::SchulteMLPFunction::New(); + //if(m_MostLikelyPathType == "polynomial") + // mlp = pct::ThirdOrderPolynomialMLPFunction::New(); + //else if (m_MostLikelyPathType == "schulte") + // mlp = pct::SchulteMLPFunction::New(); + //else + // { + // itkGenericExceptionMacro("MLP must either be schulte or polynomial, not [" << m_MostLikelyPathType << ']'); + // } // Get geometry information. We need the rotation matrix alone to rotate the direction // and the same matrix combined with mm (physical point) to voxel conversion for the volume. @@ -203,6 +204,13 @@ ProtonPairsToBackProjection if(pSOut[2]pOut[2]) pSOut = pOut + dOut * farDistOut; } + else + { + // pSIn = pOut; + // pSOut = pIn; + pSIn[2] = 0; + pSOut[2] = 0; + } } // Normalize direction with respect to z @@ -215,15 +223,18 @@ ProtonPairsToBackProjection VectorType dCurr = dIn; VectorType pCurr(0.); + itk::Matrix mlp_error; + mlp_error.Fill(0); // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); - + mlp->SetTrackerInfo(200,200,10,0.15,0.005); int totalLength = 0.; std::vector offsets; for(unsigned int k=0; k pCurr[0] = pIn[0]+z*dIn[0]; pCurr[1] = pIn[1]+z*dIn[1]; dCurr = dIn; + //mlp->EvaluateError(zmm[k],mlp_error); + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); + std_error=mlp_error(0,0); ///2 std::sqrt( + //std_error=0; } else if(pCurr[2]>=pSOut[2]) //after exit { @@ -238,6 +253,10 @@ ProtonPairsToBackProjection pCurr[0] = pSOut[0]+z*dOut[0]; pCurr[1] = pSOut[1]+z*dOut[1]; dCurr = dOut; + //mlp->EvaluateError(zmm[k],mlp_error); + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); + std_error=mlp_error(0,0); ///2 std::sqrt( + //std_error=0; } else //MLP { @@ -247,6 +266,10 @@ ProtonPairsToBackProjection mlp->Evaluate(zmm[k], pCurr[0], pCurr[1]); dCurr[0] = pCurr[0] - dCurr[0]; dCurr[1] = pCurr[1] - dCurr[1]; + //mlp->EvaluateError(zmm[k],mlp_error); + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); + std_error=mlp_error(0,0); ///2 std::sqrt( + //std_error=mlp_error(1,1); ///2 std::sqrt( } dCurr[2] *= -1.; pCurr[2] *= -1.; @@ -305,6 +328,14 @@ ProtonPairsToBackProjection offsets.push_back(this->GetOutput()->ComputeOffset(idx)); totalLength++; } + else if(m_SigmaMap && std_error==std_error) + { + typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); + m_Mutex.lock(); + imgData[ offset ] += std_error; + imgCountData[ offset ]++; + m_Mutex.unlock(); + } else { typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); diff --git a/pctProtonPairsToDistanceDrivenProjection.h b/pctProtonPairsToDistanceDrivenProjection.h index 5318617..1fe3385 100644 --- a/pctProtonPairsToDistanceDrivenProjection.h +++ b/pctProtonPairsToDistanceDrivenProjection.h @@ -68,6 +68,9 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : /** Get/Set the angle of proton pairs per pixel. */ itkGetMacro(Angle, AngleImagePointer); + /** Get/Set the MLP uncertainty of proton pairs per pixel. */ + itkGetMacro(Sigma, OutputImagePointer); + /** Get/Set the ionization potential used in the Bethe-Bloch equation. */ itkGetMacro(IonizationPotential, double); itkSetMacro(IonizationPotential, double); @@ -89,6 +92,11 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : itkSetMacro(WeightsCF, bool); itkBooleanMacro(WeightsCF); + /** Get / Set the boolean to output MLP uncertainty map. Default is off. */ + itkGetMacro(SigmaMap, bool); + itkSetMacro(SigmaMap, bool); + itkBooleanMacro(SigmaMap); + protected: ProtonPairsToDistanceDrivenProjection(); virtual ~ProtonPairsToDistanceDrivenProjection() {} @@ -126,6 +134,8 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : std::vector m_Outputs; std::vector m_AngleOutputs; std::vector m_AngleSqOutputs; + std::vector m_Sigmas; + OutputImagePointer m_Sigma; /** The two quadric functions defining the object support. */ RQIType::Pointer m_QuadricIn; @@ -143,6 +153,9 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : /** Use binning as in Collins-Fekete */ bool m_WeightsCF; + + /** Output sigma map */ + bool m_SigmaMap; }; } // end namespace pct diff --git a/pctProtonPairsToDistanceDrivenProjection.txx b/pctProtonPairsToDistanceDrivenProjection.txx index 0f17e75..0b31b1b 100644 --- a/pctProtonPairsToDistanceDrivenProjection.txx +++ b/pctProtonPairsToDistanceDrivenProjection.txx @@ -29,6 +29,10 @@ ProtonPairsToDistanceDrivenProjection m_AnglesVectors.resize( this->GetInput()->GetLargestPossibleRegion().GetNumberOfPixels() ); m_AnglesSq.resize( this->GetNumberOfWorkUnits() ); } + if(m_SigmaMap) + { + m_Sigmas.resize( this->GetNumberOfWorkUnits() ); + } if(m_QuadricOut.GetPointer()==NULL) m_QuadricOut = m_QuadricIn; @@ -48,15 +52,7 @@ ProtonPairsToDistanceDrivenProjection ::ThreadedGenerateData( const OutputImageRegionType& itkNotUsed(outputRegionForThread), rtk::ThreadIdType threadId) { // Create MLP depending on type - pct::MostLikelyPathFunction::Pointer mlp; - if(m_MostLikelyPathType == "polynomial") - mlp = pct::ThirdOrderPolynomialMLPFunction::New(); - else if (m_MostLikelyPathType == "schulte") - mlp = pct::SchulteMLPFunction::New(); - else - { - itkGenericExceptionMacro("MLP must either be schulte or polynomial, not [" << m_MostLikelyPathType << ']'); - } + pct::SchulteMLPFunction::Pointer mlp = pct::SchulteMLPFunction::New(); // Create thread image and corresponding stack to count events m_Counts[threadId] = CountImageType::New(); @@ -77,6 +73,14 @@ ProtonPairsToDistanceDrivenProjection m_AnglesSq[threadId]->FillBuffer(0); } + if( m_SigmaMap ) + { + m_Sigmas[threadId] = OutputImageType::New(); + m_Sigmas[threadId]->SetRegions(this->GetInput()->GetLargestPossibleRegion()); + m_Sigmas[threadId]->Allocate(); + m_Sigmas[threadId]->FillBuffer(0); + } + if(threadId==0) { m_Outputs[0] = this->GetOutput(); @@ -86,6 +90,10 @@ ProtonPairsToDistanceDrivenProjection m_Angle = m_Angles[0]; m_AngleSq = m_AnglesSq[0]; } + if(m_SigmaMap) + { + m_Sigma = m_Sigmas[0]; + } } else { @@ -109,11 +117,16 @@ ProtonPairsToDistanceDrivenProjection typename OutputImageType::PixelType *imgData = m_Outputs[threadId]->GetBufferPointer(); double *imgCountData = m_Counts[threadId]->GetBufferPointer(); float *imgAngleData = NULL, *imgAngleSqData = NULL; + float *imgSigmaData = NULL; if(m_ComputeScattering && !m_Robust) { imgAngleData = m_Angles[threadId]->GetBufferPointer(); imgAngleSqData = m_AnglesSq[threadId]->GetBufferPointer(); } + if(m_SigmaMap) + { + imgSigmaData = m_Sigmas[threadId]->GetBufferPointer(); + } itk::Vector imgSpacingInv; for(unsigned int i=0; i<3; i++) @@ -217,6 +230,12 @@ ProtonPairsToDistanceDrivenProjection if(pSOut[2]pOut[2]) pSOut = pOut + dOut * farDistOut; } + else + { + pSOut = pOut - dOut * pOut[2]; + pSIn[2] = 0; + pSOut[2] = 0; + } } // Normalize direction with respect to z @@ -227,13 +246,17 @@ ProtonPairsToDistanceDrivenProjection dOut[1] /= dOut[2]; //dOut[2] = 1.; SR: implicit in the following + itk::Matrix mlp_error; + mlp_error.Fill(0); // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); - + if(m_SigmaMap) + mlp->SetTrackerInfo(200,200,10,0.15,0.005); std::vector offsets; for(unsigned int k=0; k { mlp->Evaluate(zmm[k], xx, yy); } + if(m_SigmaMap) + { + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); + std_error=mlp_error(0,0); + } // Source at (0,0,args_info.source_arg), mag then to voxel xx = (xx*zmag[k] - imgOrigin[0]) * imgSpacingInv[0]; @@ -275,6 +303,11 @@ ProtonPairsToDistanceDrivenProjection imgCountData[ idx ]++; } + if(m_SigmaMap && std_error==std_error) + { + imgSigmaData[ idx ] += std_error; + } + if(m_ComputeScattering) { if(m_Robust) @@ -367,6 +400,42 @@ ProtonPairsToDistanceDrivenProjection ++itCOut; } + if(m_SigmaMap) + { + ImageIteratorType itSigmaOut(m_Sigmas[0], m_Outputs[0]->GetLargestPossibleRegion()); + + // Merge the projection computed in each thread to the first one + for(unsigned int i=1; iGetNumberOfWorkUnits(); i++) + { + if(m_Outputs[i].GetPointer() == NULL) + continue; + ImageIteratorType itSigmaOutThread(m_Sigmas[i], m_Outputs[i]->GetLargestPossibleRegion()); + + while(!itSigmaOut.IsAtEnd()) + { + itSigmaOut.Set(itSigmaOut.Get()+itSigmaOutThread.Get()); + ++itSigmaOutThread; + ++itSigmaOut; + } + + itSigmaOut.GoToBegin(); + } + + // Set count image information + m_Sigma->SetSpacing( this->GetOutput()->GetSpacing() ); + m_Sigma->SetOrigin( this->GetOutput()->GetOrigin() ); + + itCOut.GoToBegin(); + // Normalize variance with proton count (average) + while(!itCOut.IsAtEnd()) + { + if(itCOut.Get()) + itSigmaOut.Set(itSigmaOut.Get()/itCOut.Get()); + ++itSigmaOut; + ++itCOut; + } + } + if(m_ComputeScattering) { typedef itk::ImageRegionIterator ImageAngleIteratorType; @@ -450,6 +519,7 @@ ProtonPairsToDistanceDrivenProjection m_Angles.resize( 0 ); m_AnglesSq.resize( 0 ); m_AnglesVectors.resize( 0 ); + m_Sigmas.resize( 0 ); } } diff --git a/pctSchulteMLPFunction.h b/pctSchulteMLPFunction.h index e7f33cb..3904b8e 100644 --- a/pctSchulteMLPFunction.h +++ b/pctSchulteMLPFunction.h @@ -148,6 +148,13 @@ class ITK_EXPORT SchulteMLPFunction: /** Evaluate the error (x,y) (equation 27) at depth z. */ void EvaluateError( const double u1, itk::Matrix &error); + /** Set tracker parameters. */ + void SetTrackerInfo(const double dentry, const double dexit, const double dT, const double sigmap, const double xOverX0 ); + + /** Evaluate error including tracker uncertainties. */ + void EvaluateErrorWithTrackerUncertainty(const double u1, const double eIn, const double eOut, itk::Matrix &error); + + #ifdef MLP_TIMING /** Print timing information */ virtual void PrintTiming(std::ostream& os); @@ -197,6 +204,13 @@ class ITK_EXPORT SchulteMLPFunction: double m_IntForSigmaSqTTheta2; double m_IntForSigmaSqT2; + // Tracker info + double m_dentry; + double m_dexit; + double m_sigmap; + double m_dT; + double m_xOverX0; + #ifdef MLP_TIMING TimeProbe m_EvaluateProbe1; TimeProbe m_EvaluateProbe2; diff --git a/pctSchulteMLPFunction.txx b/pctSchulteMLPFunction.txx index c4bdc1b..4f5f58d 100644 --- a/pctSchulteMLPFunction.txx +++ b/pctSchulteMLPFunction.txx @@ -69,14 +69,20 @@ SchulteMLPFunction 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); + if(u1 == m_u0) + m_Sigma1 *= 0; + else + m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(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::GetValue(u1,m_u2); + if(u1 == m_u2) + m_Sigma2 *=0; + else + m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1,m_u2); #ifdef MLP_TIMING m_EvaluateProbe1.Stop(); @@ -85,21 +91,23 @@ SchulteMLPFunction // x and y, equation 24 // common calculations - InverseMatrix(m_Sigma1); - InverseMatrix(m_Sigma2); - itk::Matrix Sigma1Inv_R0 = m_Sigma1 * m_R0; - itk::Matrix R1T_Sigma2Inv = m_R1T * m_Sigma2; - itk::Matrix part(m_Sigma1 + R1T_Sigma2Inv * m_R1); - InverseMatrix(part); + itk::Matrix Sigma1_R1T = m_Sigma1 * m_R1.GetTranspose(); + itk::Matrix R1_Sigma1 = m_R1 * m_Sigma1; + + InverseMatrix(m_R1); + itk::Matrix part1 = m_R1 * m_Sigma2 + Sigma1_R1T; + InverseMatrix(part1); + itk::Matrix part2 = R1_Sigma1 + m_Sigma2 * m_R1.GetTranspose(); + InverseMatrix(part2); // x itk::Vector xMLP; - xMLP = part * (Sigma1Inv_R0 * m_x0 + R1T_Sigma2Inv * m_x2); + xMLP = m_R1 * m_Sigma2 * part1 * m_R0 * m_x0 + m_Sigma1 * part2 * m_x2; x = xMLP[0]; // y itk::Vector yMLP; - yMLP = part * (Sigma1Inv_R0 * m_y0 + R1T_Sigma2Inv * m_y2); + yMLP = m_R1 * m_Sigma2 * part1 * m_R0 * m_y0 + m_Sigma1 * part2 * m_y2; y = yMLP[0]; #ifdef MLP_TIMING @@ -111,11 +119,18 @@ void SchulteMLPFunction ::EvaluateError( const double u, itk::Matrix &error ) { - double x, y; - Evaluate(u,x,y); - error = m_Sigma1 + m_R1T * m_Sigma2 * m_R1; - InverseMatrix(error); - error *= 2.; + const double u1 = u - m_uOrigin; + if(u1 > m_u0 && u1 < m_u2) + { + double x, y; + Evaluate(u,x,y); + itk::Matrix Sigma2_InvR1T = m_Sigma2 * m_R1.GetTranspose(); + InverseMatrix(m_R1); + itk::Matrix part = Sigma2_InvR1T + m_R1 * m_Sigma1; + InverseMatrix(part); + InverseMatrix(m_R1); + error = m_Sigma1 * part * m_Sigma2 * m_R1.GetTranspose(); + } } #ifdef MLP_TIMING @@ -142,4 +157,88 @@ SchulteMLPFunction mat *= det; } +void +SchulteMLPFunction +::SetTrackerInfo(const double dentry, const double dexit, const double dT, const double sigmap, const double xOverX0 ) +{ + m_dentry = dentry* CLHEP::mm; + m_dexit = dexit * CLHEP::mm; + m_dT = dT * CLHEP::cm; + m_sigmap = sigmap * CLHEP::mm; + m_xOverX0 = xOverX0; +} + +void +SchulteMLPFunction +::EvaluateErrorWithTrackerUncertainty(const double u, const double eIn, const double eOut, itk::Matrix &error) +{ + double u1 = u - m_uOrigin; + itk::Matrix SD; + SD(0,0) = 1; + SD(0,1) = 0; + SD(1,0) = 0; + SD(1,1) = 1; + + if(u1 > m_u2) + { + SD(0,1) = std::abs(u1 - m_u2); + u1 = m_u2; + } + else if(u1 < m_u0) + { + SD(0,1) = std::abs(m_u0 - u1); + u1 = m_u0; + } + + double proton_mass_c2 = 938.272013 * CLHEP::MeV; + double x, y; + Evaluate(u1 + m_uOrigin,x,y); + + itk::Matrix Sin; + itk::Matrix Sout; + itk::Matrix SigmaIn; + itk::Matrix SigmaOut; + + double betapIn = (eIn + 2*proton_mass_c2)* eIn / (eIn + proton_mass_c2); + double betapOut = (eOut + 2*proton_mass_c2)* eOut / (eOut + proton_mass_c2); + + double sigmascIn = 13.6*CLHEP::MeV / betapIn * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); + double sigmascOut = 13.6*CLHEP::MeV / betapOut * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); + Sin(0,0) = 1; + Sin(1,1) = 1; + Sin(1,0) = 0; + Sout = Sin; + Sin(0,1) = m_dentry; + Sout(0,1) = m_dexit; + SigmaIn(1,0) = -1./m_dT; + SigmaIn(1,1) = 1./m_dT; + SigmaOut = SigmaIn; + SigmaIn(0,0) = 0; + SigmaIn(0,1) = 1; + SigmaOut(0,0) = 1; + SigmaOut(0,1) = 0; + SigmaIn = SigmaIn * SigmaIn.GetTranspose() * m_sigmap * m_sigmap; + SigmaIn(1,1) += sigmascIn * sigmascIn; + SigmaOut = SigmaOut * SigmaOut.GetTranspose() * m_sigmap * m_sigmap; + SigmaOut(1,1) += sigmascOut * sigmascOut; + + InverseMatrix(Sout); + itk::Matrix C1 = m_R0 * Sin * SigmaIn * Sin.GetTranspose() * m_R0T + m_Sigma1; + itk::Matrix C2 = m_R1 * Sout * SigmaOut * Sout.GetTranspose() * m_R1.GetTranspose() + m_R1 * m_Sigma2 * m_R1.GetTranspose(); + + itk::Matrix C1C2= C1 + C2; + InverseMatrix(C1C2); + error = C1 * C1C2 * C2; + if(u - m_uOrigin > m_u2) + { + error = SD * error * SD.GetTranspose(); + } + else if(u - m_uOrigin < m_u0) + { + InverseMatrix(SD); + error = SD * error * SD.GetTranspose(); + } + +} + } diff --git a/pctbackprojectionbinning.cxx b/pctbackprojectionbinning.cxx index 8f9f3d6..af88cb5 100644 --- a/pctbackprojectionbinning.cxx +++ b/pctbackprojectionbinning.cxx @@ -84,6 +84,7 @@ int main(int argc, char * argv[]) projection->SetIonizationPotential( args_info.ionpot_arg * CLHEP::eV ); projection->SetDisableRotation( args_info.norotation_flag ); projection->SetWeightsCF( args_info.weightsCF_flag ); + projection->SetSigmaMap( args_info.SigmaMap_flag ); // Geometry if(args_info.verbose_flag) diff --git a/pctbackprojectionbinning.ggo b/pctbackprojectionbinning.ggo index 3eeea30..e9928a7 100644 --- a/pctbackprojectionbinning.ggo +++ b/pctbackprojectionbinning.ggo @@ -17,6 +17,7 @@ option "bpVal" - "Input backprojection image values" option "bpCount" - "Input backprojection image counts" string no option "norotation" - "Bin in parallel coordinate system" flag off option "weightsCF" - "Use Collins-Fekete binning" flag off +option "SigmaMap" - "Output is an MLP error map" flag off section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no diff --git a/pctbinning.cxx b/pctbinning.cxx index 6c70fb2..0b0270e 100644 --- a/pctbinning.cxx +++ b/pctbinning.cxx @@ -45,6 +45,7 @@ int main(int argc, char * argv[]) projection->SetRobust( args_info.robust_flag ); projection->SetComputeScattering( args_info.scatwepl_given ); projection->SetWeightsCF( args_info.weightsCF_flag ); + projection->SetSigmaMap( args_info.SigmaMap_given ); if(args_info.quadricIn_given) { @@ -124,6 +125,16 @@ int main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION( cwriter->Update() ) } + if(args_info.SigmaMap_given) + { + // Write + typedef itk::ImageFileWriter< ProjectionFilter::OutputImageType > SigmaWriterType; + SigmaWriterType::Pointer sigmawriter = SigmaWriterType::New(); + sigmawriter->SetFileName( args_info.SigmaMap_arg ); + sigmawriter->SetInput( projection->GetSigma() ); + TRY_AND_EXIT_ON_ITK_EXCEPTION( sigmawriter->Update() ) + } + if(args_info.scatwepl_given) { // Write diff --git a/pctbinning.ggo b/pctbinning.ggo index dfa31b8..9622174 100644 --- a/pctbinning.ggo +++ b/pctbinning.ggo @@ -17,6 +17,7 @@ option "mlptype" - "Type of most likely path (schulte or polynomial)" 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 "weightsCF" - "Use Collins-Fekete binning" flag off +option "SigmaMap" - "Image of MLP uncertainty" string no section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no From 33943674823633a7dafed70b1c60299fa5e1bfe7 Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Mon, 15 Feb 2021 16:41:09 +0100 Subject: [PATCH 7/8] corrections to MLP Krah + add rotation before intersection for non circular objects --- AddTrackerUncertainty.py | 3 +- pctProtonPairsToBackProjection.txx | 42 +++---- pctProtonPairsToDistanceDrivenProjection.h | 25 +++- pctProtonPairsToDistanceDrivenProjection.txx | 116 +++++++++++++++---- pctSchulteMLPFunction.h | 8 +- pctSchulteMLPFunction.txx | 30 +++-- pctbinning.cxx | 9 ++ pctbinning.ggo | 3 + 8 files changed, 176 insertions(+), 60 deletions(-) diff --git a/AddTrackerUncertainty.py b/AddTrackerUncertainty.py index b365fb6..5098643 100755 --- a/AddTrackerUncertainty.py +++ b/AddTrackerUncertainty.py @@ -64,9 +64,8 @@ def GetSigmaSc(energy, xOverX0): pairs[:,2,0:2] += dYuncertEntry[:,1,:] # Entrance direction X/Y pairs[:,1,0:2] += dYuncertExit[:,0,:] # Exit position X/Y pairs[:,3,0:2] += dYuncertExit[:,1,:] # Exit direction X/Y - + itk.imwrite(itk.GetImageFromArray(pairs, is_vector=True), output) if __name__ == '__main__': AddTrackerUncertainty() - diff --git a/pctProtonPairsToBackProjection.txx b/pctProtonPairsToBackProjection.txx index b1c4f5b..796f7fd 100644 --- a/pctProtonPairsToBackProjection.txx +++ b/pctProtonPairsToBackProjection.txx @@ -192,10 +192,26 @@ ProtonPairsToBackProjection VectorType pSIn = pIn; VectorType pSOut = pOut; double nearDistIn, nearDistOut, farDistIn, farDistOut; + + VectorType pRotIn (0.); + VectorType dRotIn (0.); + VectorType pRotOut(0.); + VectorType dRotOut (0.); + for(unsigned int i=0; i<3; i++) + { + for(unsigned int j=0; j<3; j++) + { + pRotIn[i] += rotMat[i][j] * pIn[j]; + dRotIn[i] += rotMat[i][j] * dIn[j]; + pRotOut[i] += rotMat[i][j] * pOut[j]; + dRotOut[i] += rotMat[i][j] * dOut[j]; + } + } + if(m_QuadricIn.GetPointer()!=NULL) { - if(m_QuadricIn->IsIntersectedByRay(pIn,dIn,nearDistIn,farDistIn) && - m_QuadricOut->IsIntersectedByRay(pOut,dOut,nearDistOut,farDistOut)) + if(m_QuadricIn->IsIntersectedByRay(pRotIn,dRotIn,nearDistIn,farDistIn) && + m_QuadricOut->IsIntersectedByRay(pRotOut,dRotOut,nearDistOut,farDistOut)) { pSIn = pIn + dIn * nearDistIn; if(pSIn[2]pOut[2]) @@ -228,13 +244,11 @@ ProtonPairsToBackProjection // Init MLP before mm to voxel conversion mlp->Init(pSIn, pSOut, dIn, dOut); - mlp->SetTrackerInfo(200,200,10,0.15,0.005); int totalLength = 0.; std::vector offsets; for(unsigned int k=0; k pCurr[0] = pIn[0]+z*dIn[0]; pCurr[1] = pIn[1]+z*dIn[1]; dCurr = dIn; - //mlp->EvaluateError(zmm[k],mlp_error); - mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); - std_error=mlp_error(0,0); ///2 std::sqrt( - //std_error=0; } else if(pCurr[2]>=pSOut[2]) //after exit { @@ -253,10 +263,6 @@ ProtonPairsToBackProjection pCurr[0] = pSOut[0]+z*dOut[0]; pCurr[1] = pSOut[1]+z*dOut[1]; dCurr = dOut; - //mlp->EvaluateError(zmm[k],mlp_error); - mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); - std_error=mlp_error(0,0); ///2 std::sqrt( - //std_error=0; } else //MLP { @@ -266,10 +272,6 @@ ProtonPairsToBackProjection mlp->Evaluate(zmm[k], pCurr[0], pCurr[1]); dCurr[0] = pCurr[0] - dCurr[0]; dCurr[1] = pCurr[1] - dCurr[1]; - //mlp->EvaluateError(zmm[k],mlp_error); - mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); - std_error=mlp_error(0,0); ///2 std::sqrt( - //std_error=mlp_error(1,1); ///2 std::sqrt( } dCurr[2] *= -1.; pCurr[2] *= -1.; @@ -328,14 +330,6 @@ ProtonPairsToBackProjection offsets.push_back(this->GetOutput()->ComputeOffset(idx)); totalLength++; } - else if(m_SigmaMap && std_error==std_error) - { - typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); - m_Mutex.lock(); - imgData[ offset ] += std_error; - imgCountData[ offset ]++; - m_Mutex.unlock(); - } else { typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx); diff --git a/pctProtonPairsToDistanceDrivenProjection.h b/pctProtonPairsToDistanceDrivenProjection.h index 1fe3385..7bf5fad 100644 --- a/pctProtonPairsToDistanceDrivenProjection.h +++ b/pctProtonPairsToDistanceDrivenProjection.h @@ -26,7 +26,7 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : typedef itk::Image ProtonPairsImageType; typedef ProtonPairsImageType::Pointer ProtonPairsImagePointer; - typedef itk::Image CountImageType; + typedef itk::Image CountImageType; typedef CountImageType::Pointer CountImagePointer; typedef itk::Image AngleImageType; @@ -37,7 +37,8 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : typedef typename OutputImageType::RegionType OutputImageRegionType; typedef rtk::QuadricShape RQIType; - + typedef rtk::ThreeDCircularProjectionGeometry GeometryType; + typedef typename GeometryType::Pointer GeometryPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); @@ -97,6 +98,17 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : itkSetMacro(SigmaMap, bool); itkBooleanMacro(SigmaMap); + itkGetMacro(Geometry, GeometryPointer); + itkSetMacro(Geometry, GeometryPointer); + + itkGetMacro(Index, int); + itkSetMacro(Index, int); + + /** Get / Set the boolean to use MLP Krah. Default is off. */ + itkGetMacro(MLPKrah, bool); + itkSetMacro(MLPKrah, bool); + itkBooleanMacro(MLPKrah); + protected: ProtonPairsToDistanceDrivenProjection(); virtual ~ProtonPairsToDistanceDrivenProjection() {} @@ -156,6 +168,15 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection : /** Output sigma map */ bool m_SigmaMap; + + /** RTK geometry object */ + GeometryPointer m_Geometry; + + /** Index of the projection in the geometry object */ + int m_Index; + + /** Use MLP of Krah */ + bool m_MLPKrah; }; } // end namespace pct diff --git a/pctProtonPairsToDistanceDrivenProjection.txx b/pctProtonPairsToDistanceDrivenProjection.txx index 0b31b1b..8d83004 100644 --- a/pctProtonPairsToDistanceDrivenProjection.txx +++ b/pctProtonPairsToDistanceDrivenProjection.txx @@ -115,7 +115,7 @@ ProtonPairsToDistanceDrivenProjection const unsigned long npixelsPerSlice = imgSize[0] * imgSize[1]; typename OutputImageType::PixelType *imgData = m_Outputs[threadId]->GetBufferPointer(); - double *imgCountData = m_Counts[threadId]->GetBufferPointer(); + float *imgCountData = m_Counts[threadId]->GetBufferPointer(); float *imgAngleData = NULL, *imgAngleSqData = NULL; float *imgSigmaData = NULL; if(m_ComputeScattering && !m_Robust) @@ -148,9 +148,12 @@ ProtonPairsToDistanceDrivenProjection zmag[i] = (m_SourceDistance==0.)?1:(zPlaneOutInMM-m_SourceDistance)/(zmm[i]-m_SourceDistance); } - // Process pairs - while(!it.IsAtEnd()) - { + GeometryType::ThreeDHomogeneousMatrixType rotMat; + rotMat = m_Geometry->GetRotationMatrices()[m_Index].GetInverse(); + +// Process pairs + while(!it.IsAtEnd() ) // + { if(threadId==0 && it.GetIndex()[1]%10000==0) { std::cout << '\r' @@ -158,7 +161,6 @@ ProtonPairsToDistanceDrivenProjection << 100*it.GetIndex()[1]/region.GetSize(1) << "%) in thread 1" << std::flush; } - VectorType pIn = it.Get(); ++it; VectorType pOut = it.Get(); @@ -198,6 +200,7 @@ ProtonPairsToDistanceDrivenProjection const double eIn = it.Get()[0]; const double eOut = it.Get()[1]; + double value = 0.; if(eIn==0.) value = eOut; // Directly read WEPL @@ -213,28 +216,53 @@ ProtonPairsToDistanceDrivenProjection nucInfo = it.Get(); ++it; } + std::vector offsets; + + VectorType pRotIn (0.); + VectorType dRotIn (0.); + VectorType pRotOut(0.); + VectorType dRotOut (0.); + for(unsigned int i=0; i<3; i++) + { + for(unsigned int j=0; j<3; j++) + { + pRotIn[i] += rotMat[i][j] * pIn[j]; + dRotIn[i] += rotMat[i][j] * dIn[j]; + pRotOut[i] += rotMat[i][j] * pOut[j]; + dRotOut[i] += rotMat[i][j] * dOut[j]; + } + } // Move straight to entrance and exit shapes VectorType pSIn = pIn; VectorType pSOut = pOut; double nearDistIn, nearDistOut, farDistIn, farDistOut; + double dHullIn, dHullOut; + bool SLP=0; if(m_QuadricIn.GetPointer()!=NULL) { - if(m_QuadricIn->IsIntersectedByRay(pIn,dIn,nearDistIn,farDistIn) && - m_QuadricOut->IsIntersectedByRay(pOut,dOut,nearDistOut,farDistOut)) + if(m_QuadricIn->IsIntersectedByRay(pRotIn,dRotIn,nearDistIn,farDistIn) && + m_QuadricOut->IsIntersectedByRay(pRotOut,dRotOut,nearDistOut,farDistOut)) { pSIn = pIn + dIn * nearDistIn; + dHullIn=std::abs(nearDistIn); if(pSIn[2]pOut[2]) + { pSIn = pIn + dIn * farDistIn; + dHullIn=std::abs(farDistIn); + } pSOut = pOut + dOut * nearDistOut; + dHullOut=std::abs(nearDistOut); if(pSOut[2]pOut[2]) + { pSOut = pOut + dOut * farDistOut; + dHullOut=std::abs(farDistOut); + } + } else { - pSOut = pOut - dOut * pOut[2]; - pSIn[2] = 0; - pSOut[2] = 0; + SLP = 1; } } @@ -249,37 +277,80 @@ ProtonPairsToDistanceDrivenProjection itk::Matrix mlp_error; mlp_error.Fill(0); // Init MLP before mm to voxel conversion + + double xSIn = pSIn[0]; + double ySIn = pSIn[1]; + double xSOut = pSOut[0]; + double ySOut = pSOut[1]; + double dxSIn, dySIn, dxSOut, dySOut; + VectorType dSIn = dIn; + VectorType dSOut = dOut; mlp->Init(pSIn, pSOut, dIn, dOut); - if(m_SigmaMap) - mlp->SetTrackerInfo(200,200,10,0.15,0.005); - std::vector offsets; + if(m_MLPKrah && !SLP) + { + mlp->SetTrackerInfo(std::abs(dHullIn),std::abs(dHullOut),10,0.065817931,0.005); + mlp->EvaluateErrorWithTrackerUncertainty(pSIn[2],eIn, eOut, mlp_error,xSIn,ySIn, dxSIn, dySIn); + mlp->EvaluateErrorWithTrackerUncertainty(pSOut[2],eIn, eOut, mlp_error,xSOut,ySOut,dxSOut,dySOut); + dSIn[0] = std::tan(dxSIn); + dSIn[1] = std::tan(dySIn); + dSOut[0] = std::tan(dxSOut); + dSOut[1] = std::tan(dySOut); + pSIn[0] = xSIn; + pSIn[1] = ySIn; + pSOut[0] = xSOut; + pSOut[1] = ySOut; + } for(unsigned int k=0; kEvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error,dummy_xx,dummy_yy, dummy_dxx, dummy_dyy); + std_error=mlp_error(0,0)*zmag[k]*zmag[k]; + } + } else if(dk>=pSOut[2]) //after exit { const double z = (dk-pSOut[2]); - xx = pSOut[0]+z*dOut[0]; - yy = pSOut[1]+z*dOut[1]; + xx = pSOut[0]+z*dSOut[0]; + yy = pSOut[1]+z*dSOut[1]; + if(m_MLPKrah) + { + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error,dummy_xx,dummy_yy, dummy_dxx, dummy_dyy); + std_error=mlp_error(0,0)*zmag[k]*zmag[k]; + } } else //MLP { - mlp->Evaluate(zmm[k], xx, yy); - } - if(m_SigmaMap) - { - mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error); - std_error=mlp_error(0,0); + if(m_MLPKrah) + { + mlp->EvaluateErrorWithTrackerUncertainty(zmm[k],eIn, eOut, mlp_error,xx,yy, dummy_dxx, dummy_dyy); + std_error=mlp_error(0,0)*zmag[k]*zmag[k]; + } + else + { + mlp->Evaluate(zmm[k], xx, yy); + mlp->EvaluateError(zmm[k],mlp_error); + std_error=mlp_error(0,0)*zmag[k]*zmag[k]; + } } + } // Source at (0,0,args_info.source_arg), mag then to voxel xx = (xx*zmag[k] - imgOrigin[0]) * imgSpacingInv[0]; @@ -327,6 +398,7 @@ ProtonPairsToDistanceDrivenProjection } } } + if(m_WeightsCF) { std::vector channels(offsets); diff --git a/pctSchulteMLPFunction.h b/pctSchulteMLPFunction.h index 3904b8e..4417a88 100644 --- a/pctSchulteMLPFunction.h +++ b/pctSchulteMLPFunction.h @@ -55,7 +55,11 @@ class ConstantPartOfIntegrals // Constant out of integral, use [Schulte, 2008] and not [Williams, 2004] static const double c = 13.6*CLHEP::MeV * 13.6*CLHEP::MeV / (36.1*CLHEP::cm); static const double invX0 = 1./(36.1*CLHEP::cm); - double v = 1.+0.038*std::log((uy-ux)*invX0); + double v; + if((uy-ux)<1e-3) + v = 1.+0.038*std::log(1e-3*invX0); + else + v = 1.+0.038*std::log((uy-ux)*invX0); return c*v*v; } }; @@ -152,7 +156,7 @@ class ITK_EXPORT SchulteMLPFunction: void SetTrackerInfo(const double dentry, const double dexit, const double dT, const double sigmap, const double xOverX0 ); /** Evaluate error including tracker uncertainties. */ - void EvaluateErrorWithTrackerUncertainty(const double u1, const double eIn, const double eOut, itk::Matrix &error); + void EvaluateErrorWithTrackerUncertainty(const double u1, const double eIn, const double eOut, itk::Matrix &error,double &x, double&y, double &dxt, double &dyt); #ifdef MLP_TIMING diff --git a/pctSchulteMLPFunction.txx b/pctSchulteMLPFunction.txx index 4f5f58d..76f134b 100644 --- a/pctSchulteMLPFunction.txx +++ b/pctSchulteMLPFunction.txx @@ -71,7 +71,7 @@ SchulteMLPFunction m_Sigma1(0,0) = u1 * ( 2*m_Sigma1(0,1) - u1*m_Sigma1(1,1) ) + intForSigmaSqT1/* - m_IntForSigmaSqT0*/; if(u1 == m_u0) m_Sigma1 *= 0; - else + else m_Sigma1 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(m_u0,u1); // Construct Sigma2 (equations 15-18) @@ -81,7 +81,7 @@ SchulteMLPFunction m_Sigma2(0,0) = m_u2 * ( 2*m_Sigma2(0,1) - m_u2*m_Sigma2(1,1) ) + m_IntForSigmaSqT2 - intForSigmaSqT1; if(u1 == m_u2) m_Sigma2 *=0; - else + else m_Sigma2 *= Functor::SchulteMLP::ConstantPartOfIntegrals::GetValue(u1,m_u2); #ifdef MLP_TIMING @@ -115,6 +115,7 @@ SchulteMLPFunction #endif } + void SchulteMLPFunction ::EvaluateError( const double u, itk::Matrix &error ) @@ -131,6 +132,8 @@ SchulteMLPFunction InverseMatrix(m_R1); error = m_Sigma1 * part * m_Sigma2 * m_R1.GetTranspose(); } + else + error = m_Sigma1 * 0; } #ifdef MLP_TIMING @@ -170,7 +173,7 @@ SchulteMLPFunction void SchulteMLPFunction -::EvaluateErrorWithTrackerUncertainty(const double u, const double eIn, const double eOut, itk::Matrix &error) +::EvaluateErrorWithTrackerUncertainty(const double u, const double eIn, const double eOut, itk::Matrix &error,double &xt, double&yt, double &dxt, double &dyt) { double u1 = u - m_uOrigin; itk::Matrix SD; @@ -184,7 +187,7 @@ SchulteMLPFunction SD(0,1) = std::abs(u1 - m_u2); u1 = m_u2; } - else if(u1 < m_u0) + else if(u1 < m_u0) { SD(0,1) = std::abs(m_u0 - u1); u1 = m_u0; @@ -202,8 +205,9 @@ SchulteMLPFunction double betapIn = (eIn + 2*proton_mass_c2)* eIn / (eIn + proton_mass_c2); double betapOut = (eOut + 2*proton_mass_c2)* eOut / (eOut + proton_mass_c2); - double sigmascIn = 13.6*CLHEP::MeV / betapIn * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); - double sigmascOut = 13.6*CLHEP::MeV / betapOut * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); + double sigmascIn = 13.6*CLHEP::MeV / betapIn * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); + double sigmascOut = 13.6*CLHEP::MeV / betapOut * std::sqrt(m_xOverX0) * (1 + 0.038 * std::log(m_xOverX0)); + Sin(0,0) = 1; Sin(1,1) = 1; Sin(1,0) = 0; @@ -218,7 +222,7 @@ SchulteMLPFunction SigmaOut(0,0) = 1; SigmaOut(0,1) = 0; SigmaIn = SigmaIn * SigmaIn.GetTranspose() * m_sigmap * m_sigmap; - SigmaIn(1,1) += sigmascIn * sigmascIn; + SigmaIn(1,1) += sigmascIn * sigmascIn; SigmaOut = SigmaOut * SigmaOut.GetTranspose() * m_sigmap * m_sigmap; SigmaOut(1,1) += sigmascOut * sigmascOut; @@ -233,11 +237,21 @@ SchulteMLPFunction { error = SD * error * SD.GetTranspose(); } - else if(u - m_uOrigin < m_u0) + else if(u - m_uOrigin < m_u0) { InverseMatrix(SD); error = SD * error * SD.GetTranspose(); } + // x + itk::Vector xMLP; + xMLP = C2 * C1C2 * m_R0 * m_x0 + C1 * C1C2 * m_R1 * m_x2; + xt = xMLP[0]; + dxt = xMLP[1]; + // y + itk::Vector yMLP; + yMLP = C2 * C1C2 * m_R0 * m_y0 + C1 * C1C2 * m_R1 * m_y2; + yt = yMLP[0]; + dyt= yMLP[1]; } diff --git a/pctbinning.cxx b/pctbinning.cxx index 0b0270e..1552956 100644 --- a/pctbinning.cxx +++ b/pctbinning.cxx @@ -3,6 +3,7 @@ #include #include #include +#include #include "pctProtonPairsToDistanceDrivenProjection.h" #include "SmallHoleFiller.h" @@ -46,6 +47,14 @@ int main(int argc, char * argv[]) projection->SetComputeScattering( args_info.scatwepl_given ); projection->SetWeightsCF( args_info.weightsCF_flag ); projection->SetSigmaMap( args_info.SigmaMap_given ); + projection->SetIndex( args_info.Index_arg ); + projection->SetMLPKrah( args_info.krah_flag ); + + rtk::ThreeDCircularProjectionGeometryXMLFileReader::Pointer geometryReader; + geometryReader = rtk::ThreeDCircularProjectionGeometryXMLFileReader::New(); + geometryReader->SetFilename(args_info.geometry_arg); + TRY_AND_EXIT_ON_ITK_EXCEPTION( geometryReader->GenerateOutputInformation() ) + projection->SetGeometry( geometryReader->GetOutputObject() ); if(args_info.quadricIn_given) { diff --git a/pctbinning.ggo b/pctbinning.ggo index 9622174..544c16f 100644 --- a/pctbinning.ggo +++ b/pctbinning.ggo @@ -18,6 +18,7 @@ option "ionpot" - "Ionization potential used in the reconstruction in eV" option "fill" - "Fill holes, i.e. pixels that were not hit by protons" flag off option "weightsCF" - "Use Collins-Fekete binning" flag off option "SigmaMap" - "Image of MLP uncertainty" string no +option "krah" - "Use MLP of Krah" flag off section "Projections parameters" option "origin" - "Origin (default=centered)" double multiple no @@ -25,3 +26,5 @@ option "dimension" - "Dimension" int multiple no default="256 option "spacing" - "Spacing" double multiple no default="1" option "direction" - "Direction" double multiple no option "like" - "Copy information from this image (origin, dimension, spacing, direction)" string no +option "geometry" - "XML geometry file name" string yes +option "Index" - "Index of projection in the geometry file" int yes From c42e139a312f938fe75584d09900e39604eaf604 Mon Sep 17 00:00:00 2001 From: Feriel Khellaf Date: Fri, 19 Feb 2021 11:48:22 +0100 Subject: [PATCH 8/8] Python code for backprojection-then-filtering reconstruction (Poludniowski 2014) --- BPF.py | 80 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100755 BPF.py diff --git a/BPF.py b/BPF.py new file mode 100755 index 0000000..0184d3d --- /dev/null +++ b/BPF.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python + +import click +import SimpleITK as sitk +import numpy as np +import matplotlib.pyplot as plt +import scipy.signal as signal +import scipy.special as special + + +@click.command() +@click.option("-i", "--input", required=True, help="Input file name") +@click.option("-o", "--output", required=True, help="Output file name") +@click.option("-m", "--m", required=True, default=500, help="Size of the reconstruction image matrix mxm") +@click.option("-k", "--k",required=True, default=500, help="Size of the region used for offset correction kxk") + +def BPF(input, output, m, k): + + def BPFfilter(N, dx): + x = y = np.arange(N) * dx - (N / 2 - 1) *dx #filter sampled such that x=y=0 is a sample otherwise artifacts appear + X , Y = np.meshgrid (x, y) + k = np.zeros(X.shape) + R = np.sqrt(X**2 + Y**2) + not_zero = R!=0 + var_k = np.pi * R[not_zero] / dx + k[not_zero] = 1./(4 * np.pi**2 * R[not_zero]**3) * (var_k**2 * special.jn(1,var_k) - np.pi * var_k/2 * (special.jn(1,var_k) * special.struve(0,var_k) - special.jn(0,var_k) * special.struve(1,var_k))) + k[R==0] = np.pi / (12 * dx**3) + return k + + corrSize = int(k) + reconSize = int(m) + imgReader = sitk.ImageFileReader() + imgReader.SetFileName(input) + sino = imgReader.Execute() + spacing = sino.GetSpacing() + offset = sino.GetOrigin() + sino = sitk.GetArrayFromImage(sino) + bpSize = sino.shape[1] + nb_proj = sino.shape[0] + dphi = np.pi / nb_proj + proj_axis = 0 + dx = spacing[-1] + + #Compute backprojection + bp = np.sum(sino, axis = proj_axis) * dphi + sino = [] + + #Compute filter on very large matrix + k = BPFfilter(corrSize, dx) + a = int((corrSize - bpSize) / 2) #to resize filter to BP size + + #Calculations for offset correction + xbp = np.arange(bpSize) * dx - (bpSize - 1) / 2 * dx + Xbp, Ybp = np.meshgrid (xbp, xbp) + xc = np.arange(corrSize) * dx - (corrSize - 1) / 2 * dx + Xc, Yc = np.meshgrid (xc, xc) + idx = np.in1d(Xc, Xbp) * np.in1d(Yc, Ybp) + idx = np.reshape(idx, Xc.shape) + xcind = np.arange(corrSize) + Xcind, Ycind = np.meshgrid(xcind, xcind) + SecondIntegral = (dx**2 * np.sum(k[corrSize - 1 - Xcind[~idx], corrSize - 1 - Ycind[~idx]] / (np.sqrt(Xc[~idx]**2 + Yc[~idx]**2)))) + b = int((bpSize - reconSize) / 2) #resize bp to recon size + + reconstruction = np.zeros(bp.shape) + for z in np.arange(bp.shape[1]): + reconstruction[:,z,:] = signal.fftconvolve(bp[:,z,:], k[a:-a,a:-a], mode='same') * dx**2 + delta = (dx**2 * np.sum(reconstruction[b:-b,z,b:-b])) * SecondIntegral + reconstruction[:,z,:] += delta + + offset_recon = -(reconSize - 1) / 2 * dx + reconstruction = reconstruction[b:-b,:,b:-b] + reconimg = sitk.GetImageFromArray(reconstruction) + reconimg.SetSpacing(spacing[:-1]) + reconimg.SetOrigin([offset_recon, offset[1],offset_recon]) + writer = sitk.ImageFileWriter() + writer.SetFileName(output) + writer.Execute(reconimg) + +if __name__ == '__main__': + BPF()