diff --git a/Anima/diffusion/mcm/CMakeLists.txt b/Anima/diffusion/mcm/CMakeLists.txt index 4168a54e9..e35bae16d 100644 --- a/Anima/diffusion/mcm/CMakeLists.txt +++ b/Anima/diffusion/mcm/CMakeLists.txt @@ -24,6 +24,7 @@ add_library(${PROJECT_NAME} target_link_libraries(${PROJECT_NAME} AnimaMCMBase AnimaSpecialFunctions + AnimaStatisticalDistributions ITKOptimizers ITKCommon ${TinyXML2_LIBRARY} diff --git a/Anima/diffusion/mcm/animaNODDICompartment.cxx b/Anima/diffusion/mcm/animaNODDICompartment.cxx index b96299c38..793e95d7b 100644 --- a/Anima/diffusion/mcm/animaNODDICompartment.cxx +++ b/Anima/diffusion/mcm/animaNODDICompartment.cxx @@ -494,7 +494,8 @@ void NODDICompartment::UpdateKappaValues() double dawsonValue = anima::EvaluateDawsonIntegral(std::sqrt(kappa), true); m_Tau1 = (1.0 / dawsonValue - 1.0) / (2.0 * kappa); m_Tau1Deriv = (1.0 - (1.0 - dawsonValue * (2.0 * kappa - 1.0)) / (2.0 * dawsonValue * dawsonValue)) / (2.0 * kappa * kappa); - anima::GetStandardWatsonSHCoefficients(kappa,m_WatsonSHCoefficients,m_WatsonSHCoefficientDerivatives); + m_WatsonDistribution.SetConcentrationParameter(kappa); + m_WatsonDistribution.GetStandardWatsonSHCoefficients(m_WatsonSHCoefficients, m_WatsonSHCoefficientDerivatives); m_ModifiedConcentration = false; } diff --git a/Anima/diffusion/mcm/animaNODDICompartment.h b/Anima/diffusion/mcm/animaNODDICompartment.h index 85667743a..8178e4723 100644 --- a/Anima/diffusion/mcm/animaNODDICompartment.h +++ b/Anima/diffusion/mcm/animaNODDICompartment.h @@ -2,6 +2,7 @@ #include #include +#include namespace anima { @@ -84,6 +85,12 @@ class ANIMAMCM_EXPORT NODDICompartment : public BaseCompartment m_CurrentBValue = -1.0; m_CurrentGradient.fill(0.0); + + itk::Vector meanAxis; + meanAxis[0] = 0.0; + meanAxis[1] = 0.0; + meanAxis[2] = 1.0; + m_WatsonDistribution.SetMeanAxis(meanAxis); } virtual ~NODDICompartment() {} @@ -112,6 +119,8 @@ class ANIMAMCM_EXPORT NODDICompartment : public BaseCompartment double m_Tau1, m_Tau1Deriv; double m_ExtraAxonalSignal, m_IntraAxonalSignal; double m_IntraAngleDerivative, m_IntraKappaDerivative, m_IntraAxialDerivative; + + anima::WatsonDistribution m_WatsonDistribution; }; } //end namespace anima diff --git a/Anima/diffusion/mcm_tools/CMakeLists.txt b/Anima/diffusion/mcm_tools/CMakeLists.txt index 911579ba5..9c086625c 100644 --- a/Anima/diffusion/mcm_tools/CMakeLists.txt +++ b/Anima/diffusion/mcm_tools/CMakeLists.txt @@ -2,6 +2,7 @@ add_subdirectory(dwi_simulation_from_mcm) add_subdirectory(generate_isoradius_ddi_surface) add_subdirectory(get_scalar_map_from_ddi) add_subdirectory(mcm_average_images) +add_subdirectory(mcm_merge_anisotropic_weights) add_subdirectory(mcm_merge_block_images) add_subdirectory(mcm_scalar_maps) add_subdirectory(mt_estimation_validation) diff --git a/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/CMakeLists.txt b/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/CMakeLists.txt new file mode 100644 index 000000000..949697d73 --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/CMakeLists.txt @@ -0,0 +1,37 @@ +if(BUILD_TOOLS) + +project(animaMCMMergeAnisotropicWeights) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + ${ITKIO_LIBRARIES} + AnimaMCM + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/animaMCMMergeAnisotropicWeights.cxx b/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/animaMCMMergeAnisotropicWeights.cxx new file mode 100644 index 000000000..e1acba647 --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/animaMCMMergeAnisotropicWeights.cxx @@ -0,0 +1,106 @@ +#include +#include +#include +#include +#include +#include + +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg inArg("i", "inputfile", "Input MCM image", true, "", "input mcm image", cmd); + TCLAP::ValueArg outArg("o", "outputfile", "Output image with combined anisotropic weights", true, "", "output conbined weight image", cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + using InputImageType = anima::MCMImage ; + using OutputImageType = itk::VectorImage ; + using InputPixelType = InputImageType::PixelType; + using OutputPixelType = OutputImageType::PixelType; + using MCModelType = anima::MultiCompartmentModel; + using MCModelPointer = MCModelType::Pointer; + + // Read input sample + anima::MCMFileReader mcmReader; + mcmReader.SetFileName(inArg.getValue()); + mcmReader.Update(); + InputImageType::Pointer mcmImage = mcmReader.GetModelVectorImage(); + MCModelPointer mcmPtr = mcmImage->GetDescriptionModel()->Clone(); + InputPixelType mcmValue; + unsigned int numCompartments = mcmPtr->GetNumberOfCompartments(); + unsigned int numIsoCompartments = mcmPtr->GetNumberOfIsotropicCompartments(); + + using InputImageIteratorType = itk::ImageRegionConstIterator ; + using OutputImageIteratorType = itk::ImageRegionIterator ; + InputImageIteratorType inItr(mcmImage, mcmImage->GetLargestPossibleRegion()); + + // Initialize output image + unsigned int nbOutputComponents = numIsoCompartments + 1; + OutputImageType::Pointer outputImage = OutputImageType::New(); + outputImage->SetRegions(mcmImage->GetLargestPossibleRegion()); + outputImage->CopyInformation(mcmImage); + outputImage->SetVectorLength(nbOutputComponents); + outputImage->Allocate(); + + OutputPixelType outputValue(nbOutputComponents); + outputValue.Fill(0.0); + outputImage->FillBuffer(outputValue); + + std::cout << "- Number of compartments in input image: " << numCompartments << std::endl; + std::cout << "- Number of compartments in output image: " << nbOutputComponents << std::endl; + + OutputImageIteratorType outItr(outputImage, outputImage->GetLargestPossibleRegion()); + + while (!outItr.IsAtEnd()) + { + mcmValue = inItr.Get(); + + bool backgroundVoxel = true; + for (unsigned int i = 0;i < mcmValue.GetNumberOfElements();++i) + { + if (mcmValue[i] != 0.0) + { + backgroundVoxel = false; + break; + } + } + + if (backgroundVoxel) + { + ++inItr; + ++outItr; + continue; + } + + mcmPtr->SetModelVector(mcmValue); + + for (unsigned int i = 0;i < numIsoCompartments;++i) + outputValue[i] = mcmPtr->GetCompartmentWeight(i); + + double anisoWeight = 0.0; + for (unsigned int i = numIsoCompartments;i < numCompartments;++i) + anisoWeight += mcmPtr->GetCompartmentWeight(i); + + outputValue[numIsoCompartments] = anisoWeight; + + outItr.Set(outputValue); + + ++inItr; + ++outItr; + } + + anima::writeImage (outArg.getValue(), outputImage); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Anima/diffusion/tractography/CMakeLists.txt b/Anima/diffusion/tractography/CMakeLists.txt index 89db5f87a..026a26209 100644 --- a/Anima/diffusion/tractography/CMakeLists.txt +++ b/Anima/diffusion/tractography/CMakeLists.txt @@ -29,6 +29,7 @@ target_link_libraries(${PROJECT_NAME} AnimaDataIO AnimaOptimizers AnimaSHTools + AnimaStatisticalDistributions AnimaMCM AnimaMCMBase ) diff --git a/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.h b/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.h index d7fdea3c3..3755e39ff 100644 --- a/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.h +++ b/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.h @@ -1,5 +1,7 @@ #pragma once +#include + #include #include @@ -221,6 +223,12 @@ class BaseProbabilisticTractographyImageFilter : public itk::ProcessObject //! Computes additional scalar maps that are model dependent to add to the output virtual void ComputeAdditionalScalarMaps() {} + //! Holds an object of class WatsonDistribution for sampling new directions + anima::WatsonDistribution m_WatsonDistribution; + + //! Holds a sample of size 1 of directions + DirectionVectorType m_SampleOfDirections; + private: ITK_DISALLOW_COPY_AND_ASSIGN(BaseProbabilisticTractographyImageFilter); diff --git a/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.hxx b/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.hxx index 0e584e38a..48f34f17d 100644 --- a/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.hxx +++ b/Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.hxx @@ -58,6 +58,8 @@ BaseProbabilisticTractographyImageFilter m_HighestProcessedSeed = 0; m_ProgressReport = 0; + + m_SampleOfDirections.resize(1); } template diff --git a/Anima/diffusion/tractography/animaDTIProbabilisticTractographyImageFilter.cxx b/Anima/diffusion/tractography/animaDTIProbabilisticTractographyImageFilter.cxx index 426c58929..dac205467 100644 --- a/Anima/diffusion/tractography/animaDTIProbabilisticTractographyImageFilter.cxx +++ b/Anima/diffusion/tractography/animaDTIProbabilisticTractographyImageFilter.cxx @@ -4,7 +4,6 @@ #include #include #include -#include #include #include @@ -64,7 +63,10 @@ DTIProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi // else // anima::SampleFromVMFDistribution(concentrationParameter,sampling_direction,resVec,random_generator); - anima::SampleFromWatsonDistribution(concentrationParameter,sampling_direction,resVec,3,random_generator); + m_WatsonDistribution.SetMeanAxis(sampling_direction); + m_WatsonDistribution.SetConcentrationParameter(concentrationParameter); + m_WatsonDistribution.Random(m_SampleOfDirections, random_generator); + resVec = m_SampleOfDirections[0]; if (is2d) { @@ -75,10 +77,14 @@ DTIProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi if (LC > m_ThresholdForProlateTensor) { // log_prior = anima::safe_log(anima::ComputeVMFPdf(resVec, oldDirection, this->GetKappaOfPriorDistribution())); - log_prior = anima::safe_log(anima::EvaluateWatsonPDF(resVec, oldDirection, this->GetKappaOfPriorDistribution())); + m_WatsonDistribution.SetMeanAxis(oldDirection); + m_WatsonDistribution.SetConcentrationParameter(this->GetKappaOfPriorDistribution()); + log_prior = m_WatsonDistribution.GetLogDensity(resVec); // log_proposal = anima::safe_log(anima::ComputeVMFPdf(resVec, sampling_direction, concentrationParameter)); - log_proposal = anima::safe_log(anima::EvaluateWatsonPDF(resVec, sampling_direction, concentrationParameter)); + m_WatsonDistribution.SetMeanAxis(sampling_direction); + m_WatsonDistribution.SetConcentrationParameter(concentrationParameter); + log_proposal = m_WatsonDistribution.GetLogDensity(resVec); } if (anima::ComputeScalarProduct(oldDirection, resVec) < 0) @@ -199,19 +205,23 @@ double DTIProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0 if (noiseValue > 0) concentrationParameter = b0Value / std::sqrt(noiseValue); + m_WatsonDistribution.SetMeanAxis(newDirection); + if (LC > m_ThresholdForProlateTensor) { Vector3DType dtiPrincipalDirection(0.0); this->GetDTIPrincipalDirection(modelValue, dtiPrincipalDirection, is2d); - logLikelihood = std::log(anima::EvaluateWatsonPDF(dtiPrincipalDirection, newDirection, concentrationParameter)); + m_WatsonDistribution.SetConcentrationParameter(concentrationParameter); + logLikelihood = m_WatsonDistribution.GetLogDensity(dtiPrincipalDirection); } else { Vector3DType dtiMinorDirection(0.0); this->GetDTIMinorDirection(modelValue, dtiMinorDirection); - logLikelihood = std::log(anima::EvaluateWatsonPDF(dtiMinorDirection, newDirection, - concentrationParameter)); + m_WatsonDistribution.SetConcentrationParameter(-concentrationParameter); + logLikelihood = m_WatsonDistribution.GetLogDensity(dtiMinorDirection); } double resVal = log_prior + logLikelihood - log_proposal; diff --git a/Anima/diffusion/tractography/animaMCMProbabilisticTractographyImageFilter.cxx b/Anima/diffusion/tractography/animaMCMProbabilisticTractographyImageFilter.cxx index a2f9c0c7a..4ad058bb7 100644 --- a/Anima/diffusion/tractography/animaMCMProbabilisticTractographyImageFilter.cxx +++ b/Anima/diffusion/tractography/animaMCMProbabilisticTractographyImageFilter.cxx @@ -122,7 +122,10 @@ MCMProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi // else // anima::SampleFromVMFDistribution(chosenKappa,sampling_direction,resVec,random_generator); - anima::SampleFromWatsonDistribution(chosenKappa,sampling_direction,resVec,3,random_generator); + m_WatsonDistribution.SetMeanAxis(sampling_direction); + m_WatsonDistribution.SetConcentrationParameter(chosenKappa); + m_WatsonDistribution.Random(m_SampleOfDirections, random_generator); + resVec = m_SampleOfDirections[0]; if (is2d) { @@ -133,14 +136,18 @@ MCMProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi if (effectiveNumDirs > 0) { // log_prior = anima::safe_log(anima::ComputeVMFPdf(resVec, oldDirection, this->GetKappaOfPriorDistribution())); - log_prior = anima::safe_log(anima::EvaluateWatsonPDF(resVec, oldDirection, this->GetKappaOfPriorDistribution())); + m_WatsonDistribution.SetMeanAxis(oldDirection); + m_WatsonDistribution.SetConcentrationParameter(this->GetKappaOfPriorDistribution()); + log_prior = m_WatsonDistribution.GetLogDensity(resVec); log_proposal = 0; double sumWeights = 0; for (unsigned int i = 0;i < effectiveNumDirs;++i) { // log_proposal += mixtureWeights[i] * anima::ComputeVMFPdf(resVec, maximaMCM[i], kappaValues[i]); - log_proposal += mixtureWeights[i] * anima::EvaluateWatsonPDF(resVec, maximaMCM[i], kappaValues[i]); + m_WatsonDistribution.SetMeanAxis(maximaMCM[i]); + m_WatsonDistribution.SetConcentrationParameter(kappaValues[i]); + log_proposal += mixtureWeights[i] * m_WatsonDistribution.GetDensity(resVec); sumWeights += mixtureWeights[i]; } @@ -290,6 +297,9 @@ double MCMProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0 double concentrationParameter = b0Value / std::sqrt(noiseValue); + m_WatsonDistribution.SetMeanAxis(newDirection); + m_WatsonDistribution.SetConcentrationParameter(concentrationParameter); + bool oneTested = false; for (unsigned int i = workModel->GetNumberOfIsotropicCompartments();i < workModel->GetNumberOfCompartments();++i) { @@ -306,7 +316,7 @@ double MCMProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0 anima::Normalize(tmpVec,tmpVec); } - double tmpVal = std::log(anima::EvaluateWatsonPDF(tmpVec, newDirection, concentrationParameter)); + double tmpVal = m_WatsonDistribution.GetLogDensity(tmpVec); if ((tmpVal > logLikelihood)||(!oneTested)) { logLikelihood = tmpVal; diff --git a/Anima/diffusion/tractography/animaODFProbabilisticTractographyImageFilter.cxx b/Anima/diffusion/tractography/animaODFProbabilisticTractographyImageFilter.cxx index 00ba94562..564129738 100644 --- a/Anima/diffusion/tractography/animaODFProbabilisticTractographyImageFilter.cxx +++ b/Anima/diffusion/tractography/animaODFProbabilisticTractographyImageFilter.cxx @@ -109,7 +109,10 @@ ODFProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi // else // anima::SampleFromVMFDistribution(chosenKappa,sampling_direction,resVec,random_generator); - anima::SampleFromWatsonDistribution(chosenKappa,sampling_direction,resVec,3,random_generator); + m_WatsonDistribution.SetMeanAxis(sampling_direction); + m_WatsonDistribution.SetConcentrationParameter(chosenKappa); + m_WatsonDistribution.Random(m_SampleOfDirections, random_generator); + resVec = m_SampleOfDirections[0]; if (is2d) { @@ -120,13 +123,17 @@ ODFProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi if (numDirs > 0) { // log_prior = anima::safe_log( anima::ComputeVMFPdf(resVec, oldDirection, this->GetKappaOfPriorDistribution())); - log_prior = anima::safe_log( anima::EvaluateWatsonPDF(resVec, oldDirection, this->GetKappaOfPriorDistribution())); + m_WatsonDistribution.SetMeanAxis(oldDirection); + m_WatsonDistribution.SetConcentrationParameter(this->GetKappaOfPriorDistribution()); + log_prior = m_WatsonDistribution.GetLogDensity(resVec); log_proposal = 0; for (unsigned int i = 0;i < numDirs;++i) { // log_proposal += mixtureWeights[i] * anima::ComputeVMFPdf(resVec, maximaODF[i], kappaValues[i]); - log_proposal += mixtureWeights[i] * anima::EvaluateWatsonPDF(resVec, maximaODF[i], kappaValues[i]); + m_WatsonDistribution.SetMeanAxis(maximaODF[i]); + m_WatsonDistribution.SetConcentrationParameter(kappaValues[i]); + log_proposal += mixtureWeights[i] * m_WatsonDistribution.GetDensity(resVec); } log_proposal = anima::safe_log(log_proposal); @@ -234,9 +241,12 @@ double ODFProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0 double concentrationParameter = b0Value / std::sqrt(noiseValue); + m_WatsonDistribution.SetMeanAxis(newDirection); + m_WatsonDistribution.SetConcentrationParameter(concentrationParameter); + for (unsigned int i = 0;i < numDirs;++i) { - double tmpVal = std::log(anima::EvaluateWatsonPDF(maximaODF[i], newDirection, concentrationParameter)); + double tmpVal = m_WatsonDistribution.GetLogDensity(maximaODF[i]); if ((tmpVal > logLikelihood)||(i == 0)) logLikelihood = tmpVal; } diff --git a/Anima/math-tools/matrix_operations/animaVectorOperations.hxx b/Anima/math-tools/matrix_operations/animaVectorOperations.hxx index 12ea16eae..203b7c442 100644 --- a/Anima/math-tools/matrix_operations/animaVectorOperations.hxx +++ b/Anima/math-tools/matrix_operations/animaVectorOperations.hxx @@ -1,11 +1,12 @@ #include "animaVectorOperations.h" -#include - #include #include +#include + #include +#include #include namespace anima diff --git a/Anima/math-tools/special_functions/animaBesselFunctions.cxx b/Anima/math-tools/special_functions/animaBesselFunctions.cxx index 44cef8625..3fcd7a25b 100644 --- a/Anima/math-tools/special_functions/animaBesselFunctions.cxx +++ b/Anima/math-tools/special_functions/animaBesselFunctions.cxx @@ -127,7 +127,7 @@ double bessel_ratio_i_derivative_approx(double x, unsigned int N) return firstTerm - secondTerm + thirdTerm; } -double log_bessel_order_derivative_i(double x, unsigned int order, double emc, unsigned int approx_order) +double log_bessel_order_derivative_i(double x, unsigned int order, unsigned int approx_order) { double y = x; if (y > 2795.0) @@ -143,7 +143,7 @@ double log_bessel_order_derivative_i(double x, unsigned int order, double emc, u tmpVal /= std::tgamma(i+1); tmpVal /= std::tgamma(order+i+1); denom += tmpVal; - tmpVal *= anima::psi_function(order+i+1, emc); + tmpVal *= anima::psi_function(order+i+1); num += tmpVal; } diff --git a/Anima/math-tools/special_functions/animaBesselFunctions.h b/Anima/math-tools/special_functions/animaBesselFunctions.h index 93a5d604b..acfc1303b 100644 --- a/Anima/math-tools/special_functions/animaBesselFunctions.h +++ b/Anima/math-tools/special_functions/animaBesselFunctions.h @@ -27,7 +27,7 @@ ANIMASPECIALFUNCTIONS_EXPORT double bessel_ratio_i_derivative(double x, unsigned ANIMASPECIALFUNCTIONS_EXPORT double bessel_ratio_i_derivative_approx(double x, unsigned int N); //! Computes the derivative of the log of modified Bessel function of the first kind w.r.t. its order (emc is the Euler-Mascheroni constant) -ANIMASPECIALFUNCTIONS_EXPORT double log_bessel_order_derivative_i(double x, unsigned int order, double emc, unsigned int approx_order = 50); +ANIMASPECIALFUNCTIONS_EXPORT double log_bessel_order_derivative_i(double x, unsigned int order, unsigned int approx_order = 50); //! Support function for besserl_ratio_i ANIMASPECIALFUNCTIONS_EXPORT double a0r_support(double x, unsigned int N); diff --git a/Anima/math-tools/special_functions/animaGammaFunctions.cxx b/Anima/math-tools/special_functions/animaGammaFunctions.cxx index 703e2b267..89158b066 100644 --- a/Anima/math-tools/special_functions/animaGammaFunctions.cxx +++ b/Anima/math-tools/special_functions/animaGammaFunctions.cxx @@ -2,24 +2,87 @@ #include "animaGammaFunctions.h" #include +#include +#include #include namespace anima { -double psi_function(unsigned int n, double emc) +double psi_function(unsigned int n) { if (n < 1) throw itk::ExceptionObject(__FILE__, __LINE__,"The Psi function is not defined in 0.",ITK_LOCATION); - double resVal = - emc; + double resVal = digamma(1.0); for (unsigned int i = 1;i < n;++i) resVal += 1.0 / ((double)i); return resVal; } +double digamma(const double x) +{ + if (x >= 600.0) + return std::log(x - 0.5); + + if (x < 1.0e-4) + { + double gammaValue = -boost::math::digamma(1.0); + return -1.0 / x - gammaValue; + } + + return boost::math::digamma(x); +} + +double trigamma(const double x) +{ + return boost::math::trigamma(x); +} + +double inverse_digamma(const double x) +{ + if (x > digamma(600.0)) + return std::exp(x) + 0.5; + + if (x < digamma(1.0e-4)) + { + double gammaValue = -boost::math::digamma(1.0); + return -1.0 / (x + gammaValue); + } + + double gammaValue = -digamma(1.0); + double resValue = 0.0; + + if (x < -2.22) + resValue = -1.0 / (x + gammaValue); + else + resValue = std::exp(x) + 0.5; + + bool continueLoop = true; + while (continueLoop) + { + double oldResValue = resValue; + double residualValue = digamma(resValue) - x; + if (std::abs(residualValue) < std::sqrt(std::numeric_limits::epsilon())) + break; + double stepValue = resValue / 2.0; // bisection + double denomValue = trigamma(resValue); + if (std::abs(denomValue) > std::sqrt(std::numeric_limits::epsilon())) + { + double tmpStepValue = residualValue / denomValue; + if (tmpStepValue < resValue) + stepValue = tmpStepValue; + } + + resValue -= stepValue; + continueLoop = std::abs(resValue - oldResValue) > std::sqrt(std::numeric_limits::epsilon()); + } + + return resValue; +} + double gammaHalfPlusN(unsigned int n) { double resVal = std::tgamma(2*n + 1); diff --git a/Anima/math-tools/special_functions/animaGammaFunctions.h b/Anima/math-tools/special_functions/animaGammaFunctions.h index 124f688f4..fce41f3ea 100644 --- a/Anima/math-tools/special_functions/animaGammaFunctions.h +++ b/Anima/math-tools/special_functions/animaGammaFunctions.h @@ -5,9 +5,12 @@ namespace anima { - ANIMASPECIALFUNCTIONS_EXPORT double psi_function(unsigned int n, double emc); + ANIMASPECIALFUNCTIONS_EXPORT double psi_function(unsigned int n); ANIMASPECIALFUNCTIONS_EXPORT double gammaHalfPlusN(unsigned int n); ANIMASPECIALFUNCTIONS_EXPORT double gammaHalfMinusN(unsigned int n); + ANIMASPECIALFUNCTIONS_EXPORT double digamma(const double x); + ANIMASPECIALFUNCTIONS_EXPORT double trigamma(const double x); + ANIMASPECIALFUNCTIONS_EXPORT double inverse_digamma(const double x); } // end of namespace anima diff --git a/Anima/math-tools/statistical_distributions/CMakeLists.txt b/Anima/math-tools/statistical_distributions/CMakeLists.txt index 8d86d9d5d..e7405409b 100644 --- a/Anima/math-tools/statistical_distributions/CMakeLists.txt +++ b/Anima/math-tools/statistical_distributions/CMakeLists.txt @@ -23,6 +23,7 @@ add_library(${PROJECT_NAME} target_link_libraries(${PROJECT_NAME} ITKCommon + AnimaSpecialFunctions ) diff --git a/Anima/math-tools/statistical_distributions/animaBaseDistribution.cxx b/Anima/math-tools/statistical_distributions/animaBaseDistribution.cxx deleted file mode 100644 index ce78e3f9f..000000000 --- a/Anima/math-tools/statistical_distributions/animaBaseDistribution.cxx +++ /dev/null @@ -1,23 +0,0 @@ -#include "animaBaseDistribution.h" - -#include -#include - -namespace anima -{ - -void BaseDistribution::SetShapeParameter(const double val) -{ - if (val < std::numeric_limits::epsilon()) - throw itk::ExceptionObject(__FILE__, __LINE__, "The shape parameter of a statistical distribution should be strictly positive.", ITK_LOCATION); - m_ShapeParameter = val; -} - -void BaseDistribution::SetScaleParameter(const double val) -{ - if (val < std::numeric_limits::epsilon()) - throw itk::ExceptionObject(__FILE__, __LINE__, "The scale parameter of a statistical distribution should be strictly positive.", ITK_LOCATION); - m_ScaleParameter = val; -} - -} // end of namespace anima diff --git a/Anima/math-tools/statistical_distributions/animaBaseDistribution.h b/Anima/math-tools/statistical_distributions/animaBaseDistribution.h index 39117c821..f6fe649df 100644 --- a/Anima/math-tools/statistical_distributions/animaBaseDistribution.h +++ b/Anima/math-tools/statistical_distributions/animaBaseDistribution.h @@ -1,38 +1,31 @@ #pragma once -#include -#include +#include #include +#include #include namespace anima { + template class ANIMASTATISTICALDISTRIBUTIONS_EXPORT BaseDistribution { public: + using SingleValueType = TSingleValueType; + using MultipleValueType = TMultipleValueType; using GeneratorType = std::mt19937; - using VectorType = std::vector; - - BaseDistribution() - { - m_ShapeParameter = 1.0; - m_ScaleParameter = 1.0; - } - - void SetShapeParameter(const double val); - double GetShapeParameter() {return m_ShapeParameter;} - void SetScaleParameter(const double val); - double GetScaleParameter() {return m_ScaleParameter;} + BaseDistribution() {} - virtual double GetDensity(const double &x) = 0; - virtual double GetLogDensity(const double &x) = 0; - virtual void Fit(const VectorType &sample, const std::string &method) = 0; - virtual void Random(VectorType &sample, GeneratorType &generator) = 0; + virtual bool BelongsToSupport(const SingleValueType &x) = 0; + virtual double GetDensity(const SingleValueType &x) = 0; + virtual double GetLogDensity(const SingleValueType &x) = 0; + virtual void Fit(const MultipleValueType &sample, const std::string &method) = 0; + virtual void Random(MultipleValueType &sample, GeneratorType &generator) = 0; + virtual SingleValueType GetMean() = 0; + virtual double GetVariance() = 0; - private: - double m_ShapeParameter; - double m_ScaleParameter; + double GetEpsilon() {return std::sqrt(std::numeric_limits::epsilon());} }; } // end of namespace diff --git a/Anima/math-tools/statistical_distributions/animaDirichletDistribution.cxx b/Anima/math-tools/statistical_distributions/animaDirichletDistribution.cxx new file mode 100644 index 000000000..65e7496ff --- /dev/null +++ b/Anima/math-tools/statistical_distributions/animaDirichletDistribution.cxx @@ -0,0 +1,285 @@ +#include "animaDirichletDistribution.h" +#include + +#include + +#include + +namespace anima +{ + +bool DirichletDistribution::BelongsToSupport(const SingleValueType &x) +{ + unsigned int numParameters = x.size(); + + double sumValue = 0.0; + for (unsigned int i = 0;i < numParameters;++i) + { + double tmpValue = x[i]; + if (tmpValue < this->GetEpsilon() || tmpValue > 1.0 - this->GetEpsilon()) + return false; + sumValue += tmpValue; + } + + if (std::abs(sumValue - 1.0) > this->GetEpsilon()) + return false; + + return true; +} + +void DirichletDistribution::SetConcentrationParameters(const std::vector &val) +{ + unsigned int numParameters = val.size(); + + for (unsigned int i = 0;i < numParameters;++i) + { + if (val[i] < this->GetEpsilon()) + throw itk::ExceptionObject(__FILE__, __LINE__, "The concentration parameters of a statistical distribution should be strictly positive.", ITK_LOCATION); + } + + m_ConcentrationParameters = val; +} + +double DirichletDistribution::GetDensity(const SingleValueType &x) +{ + if (!this->BelongsToSupport(x)) + return 0.0; + + return std::exp(this->GetLogDensity(x)); +} + +double DirichletDistribution::GetLogDensity(const SingleValueType &x) +{ + if (!this->BelongsToSupport(x)) + throw itk::ExceptionObject(__FILE__, __LINE__, "The log-density of the Dirichlet distribution is not defined for elements outside the simplex.", ITK_LOCATION); + + unsigned int numParameters = x.size(); + + if (m_ConcentrationParameters.size() != numParameters) + throw itk::ExceptionObject(__FILE__, __LINE__, "The input argument does not belong to the same simplex as the one on which the Dirichlet distribution is defined.", ITK_LOCATION); + + double logBeta = 0.0; + double alphaZero = 0.0; + double resValue = 0.0; + for (unsigned int i = 0;i < numParameters;++i) + { + logBeta += std::lgamma(m_ConcentrationParameters[i]); + alphaZero += m_ConcentrationParameters[i]; + resValue += (m_ConcentrationParameters[i] - 1.0) * std::log(x[i]); + } + logBeta -= std::lgamma(alphaZero); + resValue -= logBeta; + + return resValue; +} + +void DirichletDistribution::Fit(const MultipleValueType &sample, const std::string &method) +{ + unsigned int numObservations = sample.rows(); + unsigned int numParameters = sample.cols(); + + std::vector alphaParameters(numParameters, 1.0); + + std::vector usefulValues(numObservations, false); + SingleValueType sampleValue(numParameters); + unsigned int numUsefulValues = 0; + for (unsigned int i = 0;i < numObservations;++i) + { + for (unsigned int j = 0;j < numParameters;++j) + sampleValue[j] = sample.get(i, j); + if (this->BelongsToSupport(sampleValue)) + { + usefulValues[i] = true; + numUsefulValues++; + } + } + + if (numUsefulValues < 2) + { + this->SetConcentrationParameters(alphaParameters); + return; + } + + std::vector logParameters(numParameters, 0.0); + for (unsigned int j = 0;j < numParameters;++j) + { + for (unsigned int i = 0;i < numObservations;++i) + { + if (!usefulValues[i]) + continue; + logParameters[j] += std::log(sample.get(i, j)); + } + + logParameters[j] /= static_cast(numUsefulValues); + } + + // Compute initial sum of alpha's + // From https://tminka.github.io/papers/dirichlet/minka-dirichlet.pdf + // there are several options. We pick the one maximizing the first term in Eq. (4) + + // First compute empirical moments of order 1 and 2 + std::vector firstMomentValues(numParameters); + std::vector secondMomentValues(numParameters); + for (unsigned int j = 0;j < numParameters;++j) + { + double firstMoment = 0.0; + double secondMoment = 0.0; + for (unsigned int i = 0;i < numObservations;++i) + { + if (!usefulValues[i]) + continue; + double tmpValue = sample.get(i, j); + firstMoment += tmpValue; + secondMoment += tmpValue * tmpValue; + } + firstMoment /= static_cast(numUsefulValues); + secondMoment /= static_cast(numUsefulValues); + + firstMomentValues[j] = firstMoment; + secondMomentValues[j] = secondMoment; + } + + double alphaSum = 0.0; + + // Then, try Eq. (23) in https://tminka.github.io/papers/dirichlet/minka-dirichlet.pdf + for (unsigned int j = 0;j < numParameters;++j) + { + double sumValue = 0.0; + bool skipCandidate = false; + for (unsigned int k = 0;k < numParameters;++k) + { + if (j == k) + continue; + + double firstMoment = firstMomentValues[k]; + double secondMoment = secondMomentValues[k]; + + double sampleVariance = secondMoment - firstMoment * firstMoment; + if (sampleVariance < this->GetEpsilon()) + { + skipCandidate = true; + break; + } + + double inLogValue = firstMoment * (1.0 - firstMoment) / sampleVariance - 1.0; + if (inLogValue < this->GetEpsilon()) + { + skipCandidate = true; + break; + } + + sumValue += std::log(inLogValue); + } + + if (skipCandidate) + continue; + + double candidateAlphaSum = std::exp(sumValue / (numParameters - 1.0)); + + if (candidateAlphaSum > alphaSum) + alphaSum = candidateAlphaSum; + } + + // Then, try Eq. (21) from https://tminka.github.io/papers/dirichlet/minka-dirichlet.pdf + for (unsigned int j = 0;j < numParameters;++j) + { + double firstMoment = firstMomentValues[j]; + double secondMoment = secondMomentValues[j]; + + double sampleVariance = secondMoment - firstMoment * firstMoment; + if (sampleVariance < this->GetEpsilon()) + continue; + + double candidateAlphaSum = (firstMoment - secondMoment) / sampleVariance; + if (candidateAlphaSum < this->GetEpsilon()) + continue; + + if (candidateAlphaSum > alphaSum) + alphaSum = candidateAlphaSum; + } + + // Finally, if no solutions, output flat Dirichlet distribution + if (alphaSum < this->GetEpsilon()) + { + this->SetConcentrationParameters(alphaParameters); + return; + } + + double digammaValue = digamma(alphaSum); + + // Fixed point iteration method + // Eq. (9) in https://tminka.github.io/papers/dirichlet/minka-dirichlet.pdf + bool continueLoop = true; + while (continueLoop) + { + double oldDigammaValue = digammaValue; + + for (unsigned int i = 0;i < numParameters;++i) + { + double psiNew = digammaValue + logParameters[i]; + alphaParameters[i] = anima::inverse_digamma(psiNew); + } + + alphaSum = std::accumulate(alphaParameters.begin(), alphaParameters.end(), 0.0); + digammaValue = digamma(alphaSum); + continueLoop = std::abs(digammaValue - oldDigammaValue) > std::sqrt(std::numeric_limits::epsilon()); + } + + m_TotalConcentration = alphaSum; + m_MeanValues.resize(numParameters); + for (unsigned int i = 0;i < numParameters;++i) + m_MeanValues[i] = alphaParameters[i] / alphaSum; + + this->SetConcentrationParameters(alphaParameters); +} + +void DirichletDistribution::Random(MultipleValueType &sample, GeneratorType &generator) +{ + unsigned int numObservations = sample.rows(); + unsigned int numParameters = sample.cols(); + + if (m_ConcentrationParameters.size() != numParameters) + throw itk::ExceptionObject(__FILE__, __LINE__, "The requested sample does not belong to the same simplex as the one on which the Dirichlet distribution is defined.", ITK_LOCATION); + + BetaDistributionType betaDist; + UniformDistributionType unifDist(0.0, 1.0); + for (unsigned int i = 0;i < numObservations;++i) + { + double samplePartialSum = 0.0; + for (unsigned int j = 0;j < numParameters - 1;++j) + { + double alphaPartialSum = 0.0; + for (unsigned int k = j + 1;k < numParameters;++k) + alphaPartialSum += m_ConcentrationParameters[k]; + + betaDist = BetaDistributionType(m_ConcentrationParameters[j], alphaPartialSum); + double phiValue = boost::math::quantile(betaDist, unifDist(generator)); + + sample(i, j) = (1.0 - samplePartialSum) * phiValue; + samplePartialSum += sample(i, j); + } + sample(i, numParameters - 1) = 1.0 - samplePartialSum; + } +} + +vnl_matrix DirichletDistribution::GetCovarianceMatrix() +{ + unsigned int numComponents = m_ConcentrationParameters.size(); + vnl_matrix covarianceMatrix(numComponents, numComponents); + + for (unsigned int i = 0;i < numComponents;++i) + { + covarianceMatrix(i, i) = m_MeanValues[i] * (1.0 - m_MeanValues[i]) / (1.0 + m_TotalConcentration); + + for (unsigned int j = i + 1;j < numComponents;++j) + { + double tmpValue = -1.0 * m_MeanValues[i] * m_MeanValues[j] / (1.0 + m_TotalConcentration); + covarianceMatrix(i, j) = tmpValue; + covarianceMatrix(j, i) = tmpValue; + } + } + + return covarianceMatrix; +} + +} // end of namespace anima diff --git a/Anima/math-tools/statistical_distributions/animaDirichletDistribution.h b/Anima/math-tools/statistical_distributions/animaDirichletDistribution.h new file mode 100644 index 000000000..eeab552e2 --- /dev/null +++ b/Anima/math-tools/statistical_distributions/animaDirichletDistribution.h @@ -0,0 +1,44 @@ +#pragma once + +#include + +#include +#include +#include +#include + +namespace anima +{ + class ANIMASTATISTICALDISTRIBUTIONS_EXPORT DirichletDistribution : public BaseDistribution,vnl_matrix> + { + public: + using UniformDistributionType = std::uniform_real_distribution; + using BetaDistributionType = boost::math::beta_distribution; + + DirichletDistribution() + { + m_ConcentrationParameters.clear(); + m_MeanValues.clear(); + m_TotalConcentration = 0.0; + } + + bool BelongsToSupport(const SingleValueType &x); + double GetDensity(const SingleValueType &x); + double GetLogDensity(const SingleValueType &x); + void Fit(const MultipleValueType &sample, const std::string &method); + void Random(MultipleValueType &sample, GeneratorType &generator); + SingleValueType GetMean() {return m_MeanValues;} + double GetVariance() {return vnl_trace(this->GetCovarianceMatrix());} + + void SetConcentrationParameters(const std::vector &val); + std::vector GetConcentrationParameters() {return m_ConcentrationParameters;} + + vnl_matrix GetCovarianceMatrix(); + + private: + std::vector m_ConcentrationParameters; + SingleValueType m_MeanValues; + double m_TotalConcentration; + }; + +} // end of namespace diff --git a/Anima/math-tools/statistical_distributions/animaDistributionSampling.h b/Anima/math-tools/statistical_distributions/animaDistributionSampling.h index a16be1754..1cd354d36 100644 --- a/Anima/math-tools/statistical_distributions/animaDistributionSampling.h +++ b/Anima/math-tools/statistical_distributions/animaDistributionSampling.h @@ -34,22 +34,6 @@ void SampleFromVMFDistribution(const ScalarType &kappa, const VectorType &meanDi template void SampleFromVMFDistributionNumericallyStable(const ScalarType &kappa, const VectorType &meanDirection, VectorType &resVec, std::mt19937 &generator); -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const VectorType &meanDirection, VectorType &resVec, unsigned int DataDimension, std::mt19937 &generator); - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const vnl_vector_fixed < ScalarType, DataDimension > &meanDirection, vnl_vector_fixed < ScalarType, DataDimension > &resVec, std::mt19937 &generator); - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const itk::Point < ScalarType, DataDimension > &meanDirection, itk::Point < ScalarType, DataDimension > &resVec, std::mt19937 &generator); - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const itk::Vector < ScalarType, DataDimension > &meanDirection, itk::Vector < ScalarType, DataDimension > &resVec, std::mt19937 &generator); - } // end of namespace anima #include "animaDistributionSampling.hxx" diff --git a/Anima/math-tools/statistical_distributions/animaDistributionSampling.hxx b/Anima/math-tools/statistical_distributions/animaDistributionSampling.hxx index 84efaffb3..a9e60acd3 100644 --- a/Anima/math-tools/statistical_distributions/animaDistributionSampling.hxx +++ b/Anima/math-tools/statistical_distributions/animaDistributionSampling.hxx @@ -172,203 +172,4 @@ void SampleFromVMFDistributionNumericallyStable(const ScalarType &kappa, const V resVec[i] += rotationMatrix(i,j) * tmpVec[j]; } -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const VectorType &meanDirection, VectorType &resVec, unsigned int DataDimension, std::mt19937 &generator) -{ - /**********************************************************************************************//** - * \fn template - * void - * SampleFromWatsonDistribution(const ScalarType &kappa, - * const VectorType &meanDirection, - * VectorType &resVec, - * unsigned int DataDimension, - * std::mt19937 &generator) - * - * \brief Sample from the Watson distribution using the procedure described in - * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.59. - * - * \author Aymeric Stamm - * \date October 2013 - * - * \param kappa Concentration parameter of the Watson distribution. - * \param meanDirection Mean direction of the Watson distribution. - * \param resVec Resulting sample. - * \param DataDimension Dimension of the sphere + 1. - * \param generator Pseudo-random number generator. - **************************************************************************************************/ - - VectorType tmpVec; - - for (unsigned int i = 0;i < DataDimension;++i) - { - tmpVec[i] = 0; - resVec[i] = 0; - } - - // Code to compute rotation matrix from vectors, taken from registration code - tmpVec[2] = 1; - - VectorType rotationNormal; - anima::ComputeCrossProduct(tmpVec, meanDirection, rotationNormal); - anima::Normalize(rotationNormal, rotationNormal); - - // Now resuming onto sampling around direction [0,0,1] - anima::TransformCartesianToSphericalCoordinates(meanDirection, tmpVec); - - double rotationAngle = tmpVec[0]; - - if (std::abs(tmpVec[2] - 1.0) > 1.0e-6) - throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION); - - double U, V, S; - if (kappa > 1.0e-6) // Bipolar distribution - { - U = SampleFromUniformDistribution(0.0, 1.0, generator); - S = 1.0 + std::log(U + (1.0 - U) * std::exp(-kappa)) / kappa; - - V = SampleFromUniformDistribution(0.0, 1.0, generator); - - if (V > 1.0e-6) - { - while (std::log(V) > kappa * S * (S - 1.0)) - { - U = SampleFromUniformDistribution(0.0, 1.0, generator); - S = 1.0 + std::log(U + (1.0 - U) * std::exp(-kappa)) / kappa; - - V = SampleFromUniformDistribution(0.0, 1.0, generator); - - if (V < 1.0e-6) - break; - } - } - } - else if (kappa < -1.0e-6) // Gridle distribution - { - double C1 = std::sqrt(std::abs(kappa)); - double C2 = std::atan(C1); - U = SampleFromUniformDistribution(0.0, 1.0, generator); - V = SampleFromUniformDistribution(0.0, 1.0, generator); - S = (1.0 / C1) * std::tan(C2 * U); - - double T = kappa * S * S; - while (V > (1.0 - T) * std::exp(T)) - { - U = SampleFromUniformDistribution(0.0, 1.0, generator); - V = SampleFromUniformDistribution(0.0, 1.0, generator); - S = (1.0 / C1) * std::tan(C2 * U); - T = kappa * S * S; - } - } - else - { - // Sampling uniformly on the sphere - S = std::cos(SampleFromUniformDistribution(0.0, M_PI, generator)); - } - - double phi = SampleFromUniformDistribution(0.0, 2.0 * M_PI, generator); - - tmpVec[0] = std::sqrt(1.0 - S*S) * std::cos(phi); - tmpVec[1] = std::sqrt(1.0 - S*S) * std::sin(phi); - tmpVec[2] = S; - - anima::RotateAroundAxis(tmpVec, rotationAngle, rotationNormal, resVec); - - double resNorm = anima::ComputeNorm(resVec); - - if (std::abs(resNorm - 1.0) > 1.0e-4) - { - std::cout << "Sampled direction norm: " << resNorm << std::endl; - std::cout << "Mean direction: " << meanDirection << std::endl; - std::cout << "Concentration parameter: " << kappa << std::endl; - throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson sampler should generate points on the 2-sphere.",ITK_LOCATION); - } - - anima::Normalize(resVec,resVec); -} - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const vnl_vector_fixed < ScalarType, DataDimension > &meanDirection, vnl_vector_fixed < ScalarType, DataDimension > &resVec, std::mt19937 &generator) -{ - /**********************************************************************************************//** - * \fn template - * void - * SampleFromWatsonDistribution(const ScalarType &kappa, - * const vnl_vector_fixed < ScalarType, DataDimension > &meanDirection, - * vnl_vector_fixed < ScalarType, DataDimension > &resVec, - * std::mt19937 &generator) - * - * \brief Sample from the Watson distribution using the procedure described in - * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.59. - * - * \author Aymeric Stamm - * \date October 2013 - * - * \param kappa Concentration parameter of the Watson distribution. - * \param meanDirection Mean direction of the Watson distribution. - * \param resVec Resulting sample. - * \param generator Pseudo-random number generator. - **************************************************************************************************/ - - resVec.fill(0.0); - SampleFromWatsonDistribution(kappa, meanDirection, resVec, DataDimension, generator); -} - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const itk::Point < ScalarType, DataDimension > &meanDirection, itk::Point < ScalarType, DataDimension > &resVec, std::mt19937 &generator) -{ - /**********************************************************************************************//** - * \fn template - * void - * SampleFromWatsonDistribution(const ScalarType &kappa, - * const itk::Point < ScalarType, DataDimension > &meanDirection, - * itk::Point < ScalarType, DataDimension > &resVec, - * std::mt19937 &generator) - * - * \brief Sample from the Watson distribution using the procedure described in - * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.59. - * - * \author Aymeric Stamm - * \date October 2013 - * - * \param kappa Concentration parameter of the Watson distribution. - * \param meanDirection Mean direction of the Watson distribution. - * \param resVec Resulting sample. - * \param generator Pseudo-random number generator. - **************************************************************************************************/ - - resVec.Fill(0.0); - SampleFromWatsonDistribution(kappa, meanDirection, resVec, DataDimension, generator); -} - -template -void -SampleFromWatsonDistribution(const ScalarType &kappa, const itk::Vector < ScalarType, DataDimension > &meanDirection, itk::Vector < ScalarType, DataDimension > &resVec, std::mt19937 &generator) -{ - /**********************************************************************************************//** - * \fn template - * void - * SampleFromWatsonDistribution(const ScalarType &kappa, - * const itk::Vector < ScalarType, DataDimension > &meanDirection, - * itk::Vector < ScalarType, DataDimension > &resVec, - * std::mt19937 &generator) - * - * \brief Sample from the Watson distribution using the procedure described in - * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.59. - * - * \author Aymeric Stamm - * \date October 2013 - * - * \param kappa Concentration parameter of the Watson distribution. - * \param meanDirection Mean direction of the Watson distribution. - * \param resVec Resulting sample. - * \param generator Pseudo-random number generator. - **************************************************************************************************/ - - resVec.Fill(0.0); - SampleFromWatsonDistribution(kappa, meanDirection, resVec, DataDimension, generator); -} - } // end of namespace anima diff --git a/Anima/math-tools/statistical_distributions/animaGammaDistribution.cxx b/Anima/math-tools/statistical_distributions/animaGammaDistribution.cxx index cacaaeeeb..af2afca23 100644 --- a/Anima/math-tools/statistical_distributions/animaGammaDistribution.cxx +++ b/Anima/math-tools/statistical_distributions/animaGammaDistribution.cxx @@ -6,31 +6,48 @@ namespace anima { -double GammaDistribution::GetDensity(const double &x) +bool GammaDistribution::BelongsToSupport(const SingleValueType &x) { - if (x < std::numeric_limits::epsilon()) + return x > std::sqrt(std::numeric_limits::epsilon()); +} + +void GammaDistribution::SetShapeParameter(const double val) +{ + if (val < std::numeric_limits::epsilon()) + throw itk::ExceptionObject(__FILE__, __LINE__, "The shape parameter of the Gamma distribution should be strictly positive.", ITK_LOCATION); + m_ShapeParameter = val; +} + +void GammaDistribution::SetScaleParameter(const double val) +{ + if (val < std::numeric_limits::epsilon()) + throw itk::ExceptionObject(__FILE__, __LINE__, "The scale parameter of the Gamma distribution should be strictly positive.", ITK_LOCATION); + m_ScaleParameter = val; +} + +double GammaDistribution::GetDensity(const SingleValueType &x) +{ + if (!this->BelongsToSupport(x)) return 0.0; return std::exp(this->GetLogDensity(x)); } -double GammaDistribution::GetLogDensity(const double &x) +double GammaDistribution::GetLogDensity(const SingleValueType &x) { - if (x < std::numeric_limits::epsilon()) + if (!this->BelongsToSupport(x)) throw itk::ExceptionObject(__FILE__, __LINE__, "The log-density of the Gamma distribution is not defined for negative or null arguments.", ITK_LOCATION); - double shapeParameter = this->GetShapeParameter(); - double scaleParameter = this->GetScaleParameter(); - double resValue = (shapeParameter - 1.0) * std::log(x); - resValue -= x / scaleParameter; - resValue -= std::lgamma(shapeParameter); - resValue -= shapeParameter * std::log(scaleParameter); + double resValue = (m_ShapeParameter - 1.0) * std::log(x); + resValue -= x / m_ScaleParameter; + resValue -= std::lgamma(m_ShapeParameter); + resValue -= m_ShapeParameter * std::log(m_ScaleParameter); return resValue; } -void GammaDistribution::Fit(const VectorType &sample, const std::string &method) +void GammaDistribution::Fit(const MultipleValueType &sample, const std::string &method) { unsigned int dimValue = sample.size(); double doubleDimValue = static_cast(dimValue); - double epsValue = std::numeric_limits::epsilon(); + double epsValue = std::sqrt(std::numeric_limits::epsilon()); double shapeParameter, scaleParameter; if (method == "mle") @@ -41,7 +58,7 @@ void GammaDistribution::Fit(const VectorType &sample, const std::string &method) for (unsigned int i = 0;i < dimValue;++i) { double tmpValue = sample[i]; - if (tmpValue < std::numeric_limits::epsilon()) + if (tmpValue < epsValue) throw itk::ExceptionObject(__FILE__, __LINE__, "Negative or null values are not allowed in a sample drawn from the Gamma distribution.", ITK_LOCATION); meanValue += tmpValue; meanLogValue += std::log(tmpValue); @@ -52,9 +69,11 @@ void GammaDistribution::Fit(const VectorType &sample, const std::string &method) double sValue = logMeanValue - meanLogValue; shapeParameter = (3.0 - sValue + std::sqrt((sValue - 3.0) * (sValue - 3.0) + 24.0 * sValue)) / (12.0 * sValue); - scaleParameter = 0.0; - if (shapeParameter > epsValue) - scaleParameter = meanValue / shapeParameter; + + if (shapeParameter < epsValue) + throw itk::ExceptionObject(__FILE__, __LINE__, "The shape parameter of a Gamma distribution cannot be negative or null.", ITK_LOCATION); + + scaleParameter = meanValue / shapeParameter; } else if (method == "biased-closed-form" || method == "unbiased-closed-form") { @@ -77,7 +96,7 @@ void GammaDistribution::Fit(const VectorType &sample, const std::string &method) scaleParameter = doubleDimValue * sumLogXValue - sumLogValue * sumValue; if (denValue > epsValue) shapeParameter /= denValue; - scaleParameter /= (doubleDimValue * doubleDimValue); + scaleParameter /= (doubleDimValue * doubleDimValue); if (method == "unbiased-closed-form") { @@ -92,9 +111,9 @@ void GammaDistribution::Fit(const VectorType &sample, const std::string &method) this->SetScaleParameter(scaleParameter); } -void GammaDistribution::Random(VectorType &sample, GeneratorType &generator) +void GammaDistribution::Random(MultipleValueType &sample, GeneratorType &generator) { - DistributionType distributionValue(this->GetShapeParameter(), this->GetScaleParameter()); + DistributionType distributionValue(m_ShapeParameter, m_ScaleParameter); unsigned int nSamples = sample.size(); for (unsigned int i = 0;i < nSamples;++i) sample[i] = distributionValue(generator); diff --git a/Anima/math-tools/statistical_distributions/animaGammaDistribution.h b/Anima/math-tools/statistical_distributions/animaGammaDistribution.h index 516a65cb2..acb7ad1a4 100644 --- a/Anima/math-tools/statistical_distributions/animaGammaDistribution.h +++ b/Anima/math-tools/statistical_distributions/animaGammaDistribution.h @@ -2,16 +2,38 @@ #include +#include + namespace anima { - class ANIMASTATISTICALDISTRIBUTIONS_EXPORT GammaDistribution : public BaseDistribution + class ANIMASTATISTICALDISTRIBUTIONS_EXPORT GammaDistribution : public BaseDistribution> { public: - using DistributionType = std::gamma_distribution<>; - double GetDensity(const double &x); - double GetLogDensity(const double &x); - void Fit(const VectorType &sample, const std::string &method); - void Random(VectorType &sample, GeneratorType &generator); + using DistributionType = std::gamma_distribution; + + GammaDistribution() + { + m_ShapeParameter = 1.0; + m_ScaleParameter = 1.0; + } + + bool BelongsToSupport(const SingleValueType &x); + double GetDensity(const SingleValueType &x); + double GetLogDensity(const SingleValueType &x); + void Fit(const MultipleValueType &sample, const std::string &method); + void Random(MultipleValueType &sample, GeneratorType &generator); + SingleValueType GetMean() {return m_ShapeParameter * m_ScaleParameter;} + double GetVariance() {return m_ShapeParameter * m_ScaleParameter * m_ScaleParameter;} + + void SetShapeParameter(const double val); + double GetShapeParameter() {return m_ShapeParameter;} + + void SetScaleParameter(const double val); + double GetScaleParameter() {return m_ScaleParameter;} + + private: + double m_ShapeParameter; + double m_ScaleParameter; }; } // end of namespace diff --git a/Anima/math-tools/statistical_distributions/animaWatsonDistribution.cxx b/Anima/math-tools/statistical_distributions/animaWatsonDistribution.cxx new file mode 100644 index 000000000..2c4c07821 --- /dev/null +++ b/Anima/math-tools/statistical_distributions/animaWatsonDistribution.cxx @@ -0,0 +1,376 @@ +#include "animaWatsonDistribution.h" + +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace anima +{ + +bool WatsonDistribution::BelongsToSupport(const SingleValueType &x) +{ + return std::abs(x.GetNorm() - 1.0) < std::sqrt(std::numeric_limits::epsilon()); +} + +void WatsonDistribution::SetMeanAxis(const itk::Vector &val) +{ + if (!this->BelongsToSupport(val)) + throw itk::ExceptionObject(__FILE__, __LINE__, "The mean axis parameter of the Watson distribution should be of unit norm.", ITK_LOCATION); + m_MeanAxis = val; +} + +void WatsonDistribution::SetConcentrationParameter(const double &val) +{ + m_ConcentrationParameter = val; +} + +double WatsonDistribution::GetDensity(const SingleValueType &x) +{ + if (!this->BelongsToSupport(x)) + return 0.0; + + // Case 1: k = 0 - Uniform distribution on the 2-sphere + if (std::abs(m_ConcentrationParameter) <= std::numeric_limits::epsilon()) + return 1.0 / (4.0 * M_PI); + + // Case 2: k > 0 - Anisotropic distribution on the 2-sphere + if (m_ConcentrationParameter > std::numeric_limits::epsilon()) + { + double kappaSqrt = std::sqrt(m_ConcentrationParameter); + double c = anima::ComputeScalarProduct(x, m_MeanAxis); + double inExp = m_ConcentrationParameter * (c * c - 1.0); + return kappaSqrt * std::exp(inExp) / (4.0 * M_PI * anima::EvaluateDawsonFunctionNR(kappaSqrt)); + } + + // Case 3: k < 0 - Gridle distribution on the 2-sphere + double Ck = std::sqrt(-m_ConcentrationParameter / M_PI) / (2.0 * M_PI * std::erf(std::sqrt(-m_ConcentrationParameter))); + double c = anima::ComputeScalarProduct(x, m_MeanAxis); + double inExp = m_ConcentrationParameter * c * c; + return Ck * std::exp(inExp); +} + +double WatsonDistribution::GetLogDensity(const SingleValueType &x) +{ + if (!this->BelongsToSupport(x)) + throw itk::ExceptionObject(__FILE__, __LINE__, "The log-density of the Watson distribution is not defined for arguments outside the 2-sphere.", ITK_LOCATION); + + return std::log(this->GetDensity(x)); +} + +double WatsonDistribution::ComputeConcentrationMLE(const double rValue, const double aValue, const double cValue, double &logLik) +{ + double bBoundValue = (rValue * cValue - aValue) / (2.0 * rValue * (1.0 - rValue)); + double rootInValue = 1.0 + (4.0 * (cValue + 1.0) * rValue * (1.0 - rValue)) / (aValue * (cValue - aValue)); + bBoundValue *= (1.0 + std::sqrt(rootInValue)); + + double lBoundValue = (rValue * cValue - aValue) / (rValue * (1.0 - rValue)); + lBoundValue *= (1.0 + (1.0 - rValue) / (cValue - aValue)); + + double uBoundValue = (rValue * cValue - aValue) / (rValue * (1.0 - rValue)); + uBoundValue *= (1.0 + rValue / aValue); + + double concentrationParameter = 0.0; + if (rValue > aValue / cValue) + concentrationParameter = (bBoundValue + lBoundValue) / 2.0; + if (rValue < aValue / cValue) + concentrationParameter = (bBoundValue + uBoundValue) / 2.0; + + logLik = concentrationParameter * (rValue - 1.0) - std::log(anima::GetScaledKummerFunctionValue(concentrationParameter, aValue, cValue)); + + return concentrationParameter; +} + +void WatsonDistribution::Fit(const MultipleValueType &sample, const std::string &method) +{ + /**********************************************************************************************//** + * \fn void Random(std::vector>, std::mt19937 &generator) + * + * \brief Closed-form approximations of the maximum likelihood estimators for the mean axis + * and the concentration parameter of the Watson distribution using the procedure + * described in Sra & Karp, The multivariate Watson distribution: Maximum-likelihood + * estimation and other aspects, Journal of Multivariate Analysis, 2013, Vol. 114, + * pp. 256-69. + * + * \author Aymeric Stamm + * \date October 2023 + * + * \param sample A numeric matrix of shape n x 3 specifying a sample of size n drawn from the + * Watson distribution on the 2-sphere. + * \param method A string specifying the estimation method. Unused here. + **************************************************************************************************/ + + unsigned int numberOfObservations = sample.size(); + itk::Vector meanAxis; + meanAxis.Fill(0.0); + meanAxis[2] = 1.0; + double concentrationParameter = 0.0; + + using MatrixType = itk::Matrix; + MatrixType scatterMatrix; + scatterMatrix.Fill(0.0); + for (unsigned int i = 0;i < numberOfObservations;++i) + { + for (unsigned int j = 0;j < m_AmbientDimension;++j) + { + for (unsigned int k = j;k < m_AmbientDimension;++k) + { + double tmpValue = sample[i][j] * sample[i][k]; + scatterMatrix(j, k) += tmpValue; + if (j != k) + scatterMatrix(k, j) += tmpValue; + } + } + } + + scatterMatrix /= static_cast(numberOfObservations); + + // Eigen decomposition of S + itk::SymmetricEigenAnalysis eigenDecomp(m_AmbientDimension); + SingleValueType eigenVals; + MatrixType eigenVecs; + eigenDecomp.ComputeEigenValuesAndVectors(scatterMatrix, eigenVals, eigenVecs); + + double aValue = 0.5; + double cValue = static_cast(m_AmbientDimension) / 2.0; + + // Compute estimate of mean axis assuming estimated concentration parameter is positive + double bipolarRValue = 0.0; + for (unsigned int i = 0;i < m_AmbientDimension;++i) + { + for (unsigned int j = 0;j < m_AmbientDimension;++j) + bipolarRValue += eigenVecs(2, i) * scatterMatrix(i, j) * eigenVecs(2, j); + } + + double bipolarLogLik; + double bipolarKappa = ComputeConcentrationMLE(bipolarRValue, aValue, cValue, bipolarLogLik); + bool bipolarValid = bipolarKappa > 0; + + // Compute estimate of mean axis assuming estimated concentration parameter is negative + double gridleRValue = 0.0; + for (unsigned int i = 0;i < m_AmbientDimension;++i) + { + for (unsigned int j = 0;j < m_AmbientDimension;++j) + gridleRValue += eigenVecs(0, i) * scatterMatrix(i, j) * eigenVecs(0, j); + } + + double gridleLogLik; + double gridleKappa = ComputeConcentrationMLE(gridleRValue, aValue, cValue, gridleLogLik); + bool gridleValid = gridleKappa < 0; + + if (bipolarValid && gridleValid) + { + if (gridleLogLik > bipolarLogLik) + { + for (unsigned int i = 0;i < m_AmbientDimension;++i) + meanAxis[i] = eigenVecs(0, i); + concentrationParameter = gridleKappa; + m_RValue = gridleRValue; + } + else + { + for (unsigned int i = 0;i < m_AmbientDimension;++i) + meanAxis[i] = eigenVecs(2, i); + concentrationParameter = bipolarKappa; + m_RValue = bipolarRValue; + } + } + else if (bipolarValid) + { + for (unsigned int i = 0;i < m_AmbientDimension;++i) + meanAxis[i] = eigenVecs(2, i); + concentrationParameter = bipolarKappa; + m_RValue = bipolarRValue; + } + else if (gridleValid) + { + for (unsigned int i = 0;i < m_AmbientDimension;++i) + meanAxis[i] = eigenVecs(0, i); + concentrationParameter = bipolarKappa; + m_RValue = gridleRValue; + } + else + { + concentrationParameter = 0.0; + m_RValue = gridleRValue; + } + + meanAxis.Normalize(); + + this->SetMeanAxis(meanAxis); + this->SetConcentrationParameter(concentrationParameter); +} + +void WatsonDistribution::Random(MultipleValueType &sample, GeneratorType &generator) +{ + /**********************************************************************************************//** + * \fn void Random(std::vector>, std::mt19937 &generator) + * + * \brief Sample from the Watson distribution using the procedure described in Fisher et al., + * Statistical Analysis of Spherical Data, Cambridge University Press, 1993, pp. 59. + * + * \author Aymeric Stamm + * \date October 2013 + * + * \param sample A numeric matrix of shape n x 3 storing a sample of size n drawn from the + * Watson distribution on the 2-sphere. + * \param generator A pseudo-random number generator. + **************************************************************************************************/ + + SingleValueType tmpVec, resVec; + tmpVec.Fill(0.0); + tmpVec[2] = 1.0; + + // Compute rotation matrix to bring [0,0,1] on meanAxis + itk::Matrix rotationMatrix = anima::GetRotationMatrixFromVectors(tmpVec, m_MeanAxis); + + // Now resuming onto sampling around direction [0,0,1] + double U, V, S; + UniformDistributionType distributionValue(0.0, 1.0); + unsigned int nSamples = sample.size(); + + for (unsigned int i = 0;i < nSamples;++i) + { + if (m_ConcentrationParameter > std::sqrt(std::numeric_limits::epsilon())) // Bipolar distribution + { + U = distributionValue(generator); + S = 1.0 + std::log(U + (1.0 - U) * std::exp(-m_ConcentrationParameter)) / m_ConcentrationParameter; + + V = distributionValue(generator); + + if (V > 1.0e-6) + { + while (std::log(V) > m_ConcentrationParameter * S * (S - 1.0)) + { + U = distributionValue(generator); + S = 1.0 + std::log(U + (1.0 - U) * std::exp(-m_ConcentrationParameter)) / m_ConcentrationParameter; + + V = distributionValue(generator); + + if (V < 1.0e-6) + break; + } + } + } + else if (m_ConcentrationParameter < -std::sqrt(std::numeric_limits::epsilon())) // Gridle distribution + { + double C1 = std::sqrt(std::abs(m_ConcentrationParameter)); + double C2 = std::atan(C1); + U = distributionValue(generator); + V = distributionValue(generator); + S = (1.0 / C1) * std::tan(C2 * U); + + double T = m_ConcentrationParameter * S * S; + while (V > (1.0 - T) * std::exp(T)) + { + U = distributionValue(generator); + V = distributionValue(generator); + S = (1.0 / C1) * std::tan(C2 * U); + T = m_ConcentrationParameter * S * S; + } + } + else // Uniform distribution + S = std::cos(M_PI * distributionValue(generator)); + + double phi = 2.0 * M_PI * distributionValue(generator); + + tmpVec[0] = std::sqrt(1.0 - S * S) * std::cos(phi); + tmpVec[1] = std::sqrt(1.0 - S * S) * std::sin(phi); + tmpVec[2] = S; + + for (unsigned int j = 0;j < m_AmbientDimension;++j) + { + resVec[j] = 0.0; + for (unsigned int k = 0;k < m_AmbientDimension;++k) + resVec[j] += rotationMatrix(j, k) * tmpVec[k]; + } + + double resNorm = resVec.Normalize(); + + if (std::abs(resNorm - 1.0) > std::sqrt(std::numeric_limits::epsilon())) + { + std::cout << "Sampled direction norm: " << resNorm << std::endl; + std::cout << "Mean direction: " << m_MeanAxis << std::endl; + std::cout << "Concentration parameter: " << m_ConcentrationParameter << std::endl; + throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson sampler should generate points on the 2-sphere.",ITK_LOCATION); + } + + sample[i] = resVec; + } +} + +void WatsonDistribution::GetStandardWatsonSHCoefficients(std::vector &coefficients, + std::vector &derivatives) +{ + // Computes the first 7 non-zero SH coefficients of the standard Watson PDF (multiplied by 4 M_PI). + const unsigned int nbCoefs = 7; + coefficients.resize(nbCoefs); + derivatives.resize(nbCoefs); + + double sqrtPi = std::sqrt(M_PI); + double k = m_ConcentrationParameter; + double dawsonValue = anima::EvaluateDawsonIntegral(std::sqrt(k), true); + double k2 = k * k; + double k3 = k2 * k; + double k4 = k3 * k; + double k5 = k4 * k; + double k6 = k5 * k; + double k7 = k6 * k; + + coefficients[0] = 2.0 * sqrtPi; + derivatives[0] = 0.0; + + if (k < 1.0e-4) // Handles small k values by Taylor series expansion + { + coefficients[1] = k / 3.0; + coefficients[2] = k2 / 35.0; + coefficients[3] = k3 / 693.0; + coefficients[4] = k4 / 19305.0; + coefficients[5] = k5 / 692835.0; + coefficients[6] = k6 / 30421755.0; + derivatives[1] = 1.0 / 3.0; + derivatives[2] = 2.0 * k / 35.0; + derivatives[3] = k2 / 231.0; + derivatives[4] = 4.0 * k3 / 19305.0; + derivatives[5] = k4 / 138567.0; + derivatives[6] = 6.0 * k5 / 30421755.0; + + for (unsigned int i = 1;i < nbCoefs;++i) + { + double tmpVal = std::pow(2.0, i + 1.0) * sqrtPi / std::sqrt(4.0 * i + 1.0); + coefficients[i] *= tmpVal; + derivatives[i] *= tmpVal; + } + + return; + } + + coefficients[1] = (3.0 - (3.0 + 2.0 * k) * dawsonValue) / (2.0 * k); + coefficients[2] = (5.0 * (-21.0 + 2.0 * k) + 3.0 * (35.0 + 4.0 * k * (5.0 + k)) * dawsonValue) / (16.0 * k2); + coefficients[3] = (21.0 * (165.0 + 4.0 * (-5.0 + k) * k) - 5.0 * (693.0 + 378.0 * k + 84.0 * k2 + 8.0 * k3) * dawsonValue) / (64.0 * k3); + coefficients[4] = (3.0 * (-225225.0 + 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + 35.0 * (19305.0 + 8.0 * k * (1287.0 + k * (297.0 + 2.0 * k * (18.0 + k)))) * dawsonValue) / (1024.0 * k4); + coefficients[5] = (11.0 * (3968055.0 + 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) - 63.0 * (692835.0 + 2.0 * k * (182325.0 + 4.0 * k * (10725.0 + 2.0 * k * (715.0 + k * (55.0 + 2.0 * k))))) * dawsonValue) / (4096.0 * k5); + coefficients[6] = (13.0 * (-540571185.0 + 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + 231.0 * (30421755.0 + 4.0 * k * (3968055.0 + k * (944775.0 + 4.0 * k * (33150.0 + k * (2925.0 + 4.0 * k * (39.0 + k)))))) * dawsonValue) / (32768.0 * k6); + + derivatives[1] = 3.0 * (-1.0 + (-1.0 + 2.0 * k) * dawsonValue + 2.0 * dawsonValue * dawsonValue) / (4.0 * k2); + derivatives[2] = 5.0 * ((21.0 - 2.0 * k) + (63.0 + 4.0 * (-11.0 + k) * k) * dawsonValue - 12.0 * (7.0 + 2.0 * k) * dawsonValue * dawsonValue) / (32.0 * k3); + derivatives[3] = 21.0 * ((-165.0 - 4.0 * (-5.0 + k) * k) + (-5.0 + 2.0 * k) * (165.0 + 4.0 * (-3.0 + k) * k) * dawsonValue + 10.0 * (99.0 + 4.0 * k * (9.0 + k)) * dawsonValue * dawsonValue) / (128.0 * k4); + derivatives[4] = 3.0 * ((225225.0 - 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + (1576575.0 + 8.0 * k * (-75075.0 + k * (10395.0 + 2.0 * k * (-978.0 + 31.0 * k)))) * dawsonValue - 840.0 * (2145.0 + 2.0 * k * (429.0 + 66.0 * k + 4.0 * k2)) * dawsonValue * dawsonValue) / (2048.0 * k5); + derivatives[5] = 11.0 * ((-3968055.0 - 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) + (-35712495.0 + 2.0 * k * (5917275.0 + 8.0 * k * (-118755.0 + k * (21060.0 + k * (-965.0 + 58.0 * k))))) * dawsonValue + 630.0 * (62985.0 + 8.0 * k * (3315.0 + k * (585.0 + 2.0 * k * (26.0 + k)))) * dawsonValue * dawsonValue) / (8192.0 * k6); + derivatives[6] = 13.0 * ((540571185.0 - 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + (5946283035.0 + 4.0 * k * (-446558805.0 + k * (79910523.0 + 4.0 * k * (-3322242.0 + k * (187341.0 + 4.0 * k * (-3765.0 + 61.0 * k)))))) * dawsonValue - 2772.0 * (2340135.0 + 2.0 * k * (508725.0 + 4.0 * k * (24225.0 + 2.0 * k * (1275.0 + k * (75.0 + 2.0 * k))))) * dawsonValue * dawsonValue) / (65536.0 * k7); + + for (unsigned int i = 1;i < nbCoefs;++i) + { + double sqrtVal = std::sqrt(1.0 + 4.0 * i); + coefficients[i] *= (sqrtPi * sqrtVal / dawsonValue); + derivatives[i] *= (sqrtPi * sqrtVal / (dawsonValue * dawsonValue)); + } +} + +} // end of namespace anima diff --git a/Anima/math-tools/statistical_distributions/animaWatsonDistribution.h b/Anima/math-tools/statistical_distributions/animaWatsonDistribution.h index 099279bae..d58907c83 100644 --- a/Anima/math-tools/statistical_distributions/animaWatsonDistribution.h +++ b/Anima/math-tools/statistical_distributions/animaWatsonDistribution.h @@ -1,27 +1,50 @@ #pragma once -#include -#include +#include + #include namespace anima -{ - - template - double EvaluateWatsonPDF(const VectorType &v, const VectorType &meanAxis, const ScalarType &kappa); - - template - double EvaluateWatsonPDF(const vnl_vector_fixed &v, const vnl_vector_fixed &meanAxis, const ScalarType &kappa); - - template - double EvaluateWatsonPDF(const itk::Point &v, const itk::Point &meanAxis, const ScalarType &kappa); - - template - double EvaluateWatsonPDF(const itk::Vector &v, const itk::Vector &meanAxis, const ScalarType &kappa); - - template - void GetStandardWatsonSHCoefficients(const ScalarType k, std::vector &coefficients, std::vector &derivatives); +{ + class ANIMASTATISTICALDISTRIBUTIONS_EXPORT WatsonDistribution : public BaseDistribution,std::vector>> + { + public: + using UniformDistributionType = std::uniform_real_distribution; + + WatsonDistribution() + { + m_MeanAxis[0] = 0; + m_MeanAxis[1] = 0; + m_MeanAxis[2] = 1; + m_ConcentrationParameter = 1.0; + m_RValue = 0.0; + } + + bool BelongsToSupport(const SingleValueType &x); + double GetDensity(const SingleValueType &x); + double GetLogDensity(const SingleValueType &x); + void Fit(const MultipleValueType &sample, const std::string &method); + void Random(MultipleValueType &sample, GeneratorType &generator); + SingleValueType GetMean() {return m_MeanAxis;} + double GetVariance() {return 1.0 - m_RValue;} + + void SetMeanAxis(const itk::Vector &x); + SingleValueType GetMeanAxis() {return m_MeanAxis;} + + void SetConcentrationParameter(const double &x); + double GetConcentrationParameter() {return m_ConcentrationParameter;} + + void GetStandardWatsonSHCoefficients( + std::vector &coefficients, + std::vector &derivatives + ); + + private: + itk::Vector m_MeanAxis; + double m_ConcentrationParameter; + double m_RValue; + double ComputeConcentrationMLE(const double rValue, const double aValue, const double cValue, double &logLik); + const unsigned int m_AmbientDimension = 3; + }; } // end of namespace anima - -#include "animaWatsonDistribution.hxx" diff --git a/Anima/math-tools/statistical_distributions/animaWatsonDistribution.hxx b/Anima/math-tools/statistical_distributions/animaWatsonDistribution.hxx deleted file mode 100644 index ea40ae519..000000000 --- a/Anima/math-tools/statistical_distributions/animaWatsonDistribution.hxx +++ /dev/null @@ -1,149 +0,0 @@ -#pragma once -#include - -#include "animaWatsonDistribution.h" -#include -#include - -#include -#include - -#include - -namespace anima -{ - -template -double EvaluateWatsonPDF(const VectorType &v, const VectorType &meanAxis, const ScalarType &kappa) -{ - /************************************************************************************************ - * \fn template - * double - * EvaluateWatsonPDF(const VectorType &v, - * const VectorType &meanAxis, - * const ScalarType &kappa) - * - * \brief Evaluate the Watson probability density function using the definition of - * Fisher et al., Statistical Analysis of Spherical Data, 1993, p.89. - * - * \author Aymeric Stamm - * \date July 2014 - * - * \param v Sample on which evaluating the PDF. - * \param meanAxis Mean axis of the Watson distribution. - * \param kappa Concentration parameter of the Watson distribution. - **************************************************************************************************/ - - if (std::abs(kappa) < 1.0e-6) - return 1.0 / (4.0 * M_PI); - else if (kappa > 0) - { - double kappaSqrt = std::sqrt(kappa); - double c = anima::ComputeScalarProduct(v, meanAxis); - double inExp = kappa * (c * c - 1.0); - return kappaSqrt * std::exp(inExp) / (4.0 * M_PI * anima::EvaluateDawsonFunctionNR(kappaSqrt)); - } - else - { - double Ck = std::sqrt(-kappa / M_PI) / (2.0 * M_PI * std::erf(std::sqrt(-kappa))); - double c = anima::ComputeScalarProduct(v, meanAxis); - double inExp = kappa * c * c; - return Ck * std::exp(inExp); - } -} - -template -double EvaluateWatsonPDF(const vnl_vector_fixed &v, const vnl_vector_fixed &meanAxis, const ScalarType &kappa) -{ - if (std::abs(v.squared_magnitude() - 1.0) > 1.0e-6 || std::abs(meanAxis.squared_magnitude() - 1.0) > 1.0e-6) - throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION); - - return EvaluateWatsonPDF, ScalarType>(v,meanAxis,kappa); -} - -template -double EvaluateWatsonPDF(const itk::Point &v, const itk::Point &meanAxis, const ScalarType &kappa) -{ - if (std::abs(v.GetVnlVector().squared_magnitude() - 1.0) > 1.0e-6 || std::abs(meanAxis.GetVnlVector().squared_magnitude() - 1.0) > 1.0e-6) - throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION); - - return EvaluateWatsonPDF, ScalarType>(v,meanAxis,kappa); -} - -template -double EvaluateWatsonPDF(const itk::Vector &v, const itk::Vector &meanAxis, const ScalarType &kappa) -{ - if (std::abs(v.GetSquaredNorm() - 1.0) > 1.0e-6 || std::abs(meanAxis.GetSquaredNorm() - 1.0) > 1.0e-6) - throw itk::ExceptionObject(__FILE__, __LINE__,"The Watson distribution is on the 2-sphere.",ITK_LOCATION); - - return EvaluateWatsonPDF, ScalarType>(v,meanAxis,kappa); -} - -template - void GetStandardWatsonSHCoefficients(const ScalarType k, std::vector &coefficients, std::vector &derivatives) -{ - // Computes the first 7 non-zero SH coefficients of the standard Watson PDF (multiplied by 4 M_PI). - const unsigned int nbCoefs = 7; - coefficients.resize(nbCoefs); - derivatives.resize(nbCoefs); - - double sqrtPi = std::sqrt(M_PI); - double dawsonValue = anima::EvaluateDawsonIntegral(std::sqrt(k), true); - double k2 = k * k; - double k3 = k2 * k; - double k4 = k3 * k; - double k5 = k4 * k; - double k6 = k5 * k; - double k7 = k6 * k; - - coefficients[0] = 2.0 * sqrtPi; - derivatives[0] = 0.0; - - if (k < 1.0e-4) // Handles small k values by Taylor series expansion - { - coefficients[1] = k / 3.0; - coefficients[2] = k2 / 35.0; - coefficients[3] = k3 / 693.0; - coefficients[4] = k4 / 19305.0; - coefficients[5] = k5 / 692835.0; - coefficients[6] = k6 / 30421755.0; - derivatives[1] = 1.0 / 3.0; - derivatives[2] = 2.0 * k / 35.0; - derivatives[3] = k2 / 231.0; - derivatives[4] = 4.0 * k3 / 19305.0; - derivatives[5] = k4 / 138567.0; - derivatives[6] = 6.0 * k5 / 30421755.0; - - for (unsigned int i = 1;i < nbCoefs;++i) - { - double tmpVal = std::pow(2.0, i + 1.0) * sqrtPi / std::sqrt(4.0 * i + 1.0); - coefficients[i] *= tmpVal; - derivatives[i] *= tmpVal; - } - - return; - } - - coefficients[1] = (3.0 - (3.0 + 2.0 * k) * dawsonValue) / (2.0 * k); - coefficients[2] = (5.0 * (-21.0 + 2.0 * k) + 3.0 * (35.0 + 4.0 * k * (5.0 + k)) * dawsonValue) / (16.0 * k2); - coefficients[3] = (21.0 * (165.0 + 4.0 * (-5.0 + k) * k) - 5.0 * (693.0 + 378.0 * k + 84.0 * k2 + 8.0 * k3) * dawsonValue) / (64.0 * k3); - coefficients[4] = (3.0 * (-225225.0 + 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + 35.0 * (19305.0 + 8.0 * k * (1287.0 + k * (297.0 + 2.0 * k * (18.0 + k)))) * dawsonValue) / (1024.0 * k4); - coefficients[5] = (11.0 * (3968055.0 + 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) - 63.0 * (692835.0 + 2.0 * k * (182325.0 + 4.0 * k * (10725.0 + 2.0 * k * (715.0 + k * (55.0 + 2.0 * k))))) * dawsonValue) / (4096.0 * k5); - coefficients[6] = (13.0 * (-540571185.0 + 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + 231.0 * (30421755.0 + 4.0 * k * (3968055.0 + k * (944775.0 + 4.0 * k * (33150.0 + k * (2925.0 + 4.0 * k * (39.0 + k)))))) * dawsonValue) / (32768.0 * k6); - - derivatives[1] = 3.0 * (-1.0 + (-1.0 + 2.0 * k) * dawsonValue + 2.0 * dawsonValue * dawsonValue) / (4.0 * k2); - derivatives[2] = 5.0 * ((21.0 - 2.0 * k) + (63.0 + 4.0 * (-11.0 + k) * k) * dawsonValue - 12.0 * (7.0 + 2.0 * k) * dawsonValue * dawsonValue) / (32.0 * k3); - derivatives[3] = 21.0 * ((-165.0 - 4.0 * (-5.0 + k) * k) + (-5.0 + 2.0 * k) * (165.0 + 4.0 * (-3.0 + k) * k) * dawsonValue + 10.0 * (99.0 + 4.0 * k * (9.0 + k)) * dawsonValue * dawsonValue) / (128.0 * k4); - derivatives[4] = 3.0 * ((225225.0 - 2.0 * k * (15015.0 + 2.0 * k * (-1925.0 + 62.0 * k))) + (1576575.0 + 8.0 * k * (-75075.0 + k * (10395.0 + 2.0 * k * (-978.0 + 31.0 * k)))) * dawsonValue - 840.0 * (2145.0 + 2.0 * k * (429.0 + 66.0 * k + 4.0 * k2)) * dawsonValue * dawsonValue) / (2048.0 * k5); - derivatives[5] = 11.0 * ((-3968055.0 - 8.0 * k * (-69615.0 + 2.0 * k * (9828.0 + k * (-468.0 + 29.0 * k)))) + (-35712495.0 + 2.0 * k * (5917275.0 + 8.0 * k * (-118755.0 + k * (21060.0 + k * (-965.0 + 58.0 * k))))) * dawsonValue + 630.0 * (62985.0 + 8.0 * k * (3315.0 + k * (585.0 + 2.0 * k * (26.0 + k)))) * dawsonValue * dawsonValue) / (8192.0 * k6); - derivatives[6] = 13.0 * ((540571185.0 - 2.0 * k * (39171825.0 + 4.0 * k * (-2909907.0 + 2.0 * k * (82467.0 + k * (-7469.0 + 122.0 * k))))) + (5946283035.0 + 4.0 * k * (-446558805.0 + k * (79910523.0 + 4.0 * k * (-3322242.0 + k * (187341.0 + 4.0 * k * (-3765.0 + 61.0 * k)))))) * dawsonValue - 2772.0 * (2340135.0 + 2.0 * k * (508725.0 + 4.0 * k * (24225.0 + 2.0 * k * (1275.0 + k * (75.0 + 2.0 * k))))) * dawsonValue * dawsonValue) / (65536.0 * k7); - - for (unsigned int i = 1;i < nbCoefs;++i) - { - double sqrtVal = std::sqrt(1.0 + 4.0 * i); - coefficients[i] *= (sqrtPi * sqrtVal / dawsonValue); - derivatives[i] *= (sqrtPi * sqrtVal / (dawsonValue * dawsonValue)); - } -} - -} // end namespace anima diff --git a/Anima/math-tools/statistical_distributions/animaSampleImageFromDistributionImageFilter.h b/Anima/math-tools/statistical_distributions/sample_image_from_distribution/animaSampleImageFromDistributionImageFilter.h similarity index 100% rename from Anima/math-tools/statistical_distributions/animaSampleImageFromDistributionImageFilter.h rename to Anima/math-tools/statistical_distributions/sample_image_from_distribution/animaSampleImageFromDistributionImageFilter.h diff --git a/Anima/math-tools/statistical_distributions/animaSampleImageFromDistributionImageFilter.hxx b/Anima/math-tools/statistical_distributions/sample_image_from_distribution/animaSampleImageFromDistributionImageFilter.hxx similarity index 100% rename from Anima/math-tools/statistical_distributions/animaSampleImageFromDistributionImageFilter.hxx rename to Anima/math-tools/statistical_distributions/sample_image_from_distribution/animaSampleImageFromDistributionImageFilter.hxx diff --git a/Anima/math-tools/statistical_distributions/watson_sh_test/CMakeLists.txt b/Anima/math-tools/statistical_distributions/watson_sh_test/CMakeLists.txt index 26b35c12f..f91fcde00 100644 --- a/Anima/math-tools/statistical_distributions/watson_sh_test/CMakeLists.txt +++ b/Anima/math-tools/statistical_distributions/watson_sh_test/CMakeLists.txt @@ -25,6 +25,7 @@ add_executable(${PROJECT_NAME} target_link_libraries(${PROJECT_NAME} AnimaSpecialFunctions + AnimaStatisticalDistributions ) ## ############################################################################# diff --git a/Anima/math-tools/statistical_distributions/watson_sh_test/animaWatsonSHTest.cxx b/Anima/math-tools/statistical_distributions/watson_sh_test/animaWatsonSHTest.cxx index c5b81c08b..c914bd3dc 100644 --- a/Anima/math-tools/statistical_distributions/watson_sh_test/animaWatsonSHTest.cxx +++ b/Anima/math-tools/statistical_distributions/watson_sh_test/animaWatsonSHTest.cxx @@ -1,4 +1,6 @@ #include +#include + #include int main(int argc, char **argv) @@ -8,7 +10,14 @@ int main(int argc, char **argv) std::cout << "Testing SH approximation for kappa : " << kappa << std::endl; std::vector coefs, derivatives; - anima::GetStandardWatsonSHCoefficients(kappa,coefs,derivatives); + anima::WatsonDistribution watsonDistr; + itk::Vector meanAxis; + meanAxis[0] = 0.0; + meanAxis[1] = 0.0; + meanAxis[2] = 1.0; + watsonDistr.SetMeanAxis(meanAxis); + watsonDistr.SetConcentrationParameter(kappa); + watsonDistr.GetStandardWatsonSHCoefficients(coefs, derivatives); std::cout << "Dawson integral: " << anima::EvaluateDawsonIntegral(std::sqrt(kappa),false) << std::endl; diff --git a/Anima/math-tools/statistics/CMakeLists.txt b/Anima/math-tools/statistics/CMakeLists.txt index 1e8c08d91..f9b5e598e 100644 --- a/Anima/math-tools/statistics/CMakeLists.txt +++ b/Anima/math-tools/statistics/CMakeLists.txt @@ -4,9 +4,12 @@ add_subdirectory(build_samples) add_subdirectory(gamma_estimation) +add_subdirectory(dirichlet_estimation) if(BUILD_TESTING) + add_subdirectory(dirichlet_estimation_test) add_subdirectory(gamma_estimation_test) + add_subdirectory(watson_estimation_test) endif() add_subdirectory(boot_strap_4d_volume) diff --git a/Anima/math-tools/statistics/dirichlet_estimation/CMakeLists.txt b/Anima/math-tools/statistics/dirichlet_estimation/CMakeLists.txt new file mode 100644 index 000000000..7a93a515f --- /dev/null +++ b/Anima/math-tools/statistics/dirichlet_estimation/CMakeLists.txt @@ -0,0 +1,37 @@ +if(BUILD_TOOLS) + +project(animaDirichletEstimation) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + AnimaStatisticalDistributions + ${ITKIO_LIBRARIES} + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/math-tools/statistics/dirichlet_estimation/animaDirichletEstimation.cxx b/Anima/math-tools/statistics/dirichlet_estimation/animaDirichletEstimation.cxx new file mode 100644 index 000000000..92f03da4a --- /dev/null +++ b/Anima/math-tools/statistics/dirichlet_estimation/animaDirichletEstimation.cxx @@ -0,0 +1,258 @@ +#include +#include +#include +#include + +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg inArg("i", "inputfiles", "Input image list in text file", true, "", "input image list", cmd); + TCLAP::ValueArg maskArg("m", "maskfiles", "Input masks list in text file (mask images should contain only zeros or ones)", false, "", "input masks list", cmd); + TCLAP::ValueArg outArg("o", "outputfile", "Output image with concentration parameters", true, "", "output kappa image", cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + using InputImageType = itk::VectorImage ; + using MaskImageType = itk::Image ; + using OutputImageType = itk::VectorImage ; + using OutputPixelType = OutputImageType::PixelType; + + // Read input sample + std::ifstream imageIn(inArg.getValue()); + char imageN[2048]; + std::vector inputImages; + while (!imageIn.eof()) + { + imageIn.getline(imageN, 2048); + if (strcmp(imageN, "") == 0) + continue; + inputImages.push_back(anima::readImage(imageN)); + } + imageIn.close(); + + unsigned int nbImages = inputImages.size(); + unsigned int nbComponents = inputImages[0]->GetVectorLength(); + + for (unsigned int i = 1;i < nbImages;++i) + { + if (inputImages[i]->GetVectorLength() != nbComponents) + { + std::cerr << "The number of components should be the same for all input images." << std::endl; + return EXIT_FAILURE; + } + } + + std::vector maskImages; + if (maskArg.getValue() != "") + { + // Read input masks + std::ifstream masksIn(maskArg.getValue()); + char maskN[2048]; + while (!masksIn.eof()) + { + masksIn.getline(maskN, 2048); + if (strcmp(maskN, "") == 0) + continue; + maskImages.push_back(anima::readImage(maskN)); + } + masksIn.close(); + + if (maskImages.size() != nbImages) + { + std::cerr << "The number of mask images should match the number of input images." << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "- Number of input images: " << nbImages << std::endl; + std::cout << "- Number of components: " << nbComponents << std::endl; + + using InputImageIteratorType = itk::ImageRegionConstIterator ; + using MaskImageIteratorType = itk::ImageRegionConstIterator ; + using OutputImageIteratorType = itk::ImageRegionIterator ; + std::vector inItrs(nbImages); + std::vector maskItrs(nbImages); + for (unsigned int i = 0;i < nbImages;++i) + { + inItrs[i] = InputImageIteratorType(inputImages[i], inputImages[i]->GetLargestPossibleRegion()); + if (maskArg.getValue() != "") + maskItrs[i] = MaskImageIteratorType(maskImages[i], maskImages[i]->GetLargestPossibleRegion()); + } + + // Initialize output image + OutputPixelType outputValue(nbComponents); + outputValue.Fill(0.0); + + OutputImageType::Pointer outputImage = OutputImageType::New(); + outputImage->SetRegions(inputImages[0]->GetLargestPossibleRegion()); + outputImage->CopyInformation(inputImages[0]); + outputImage->Allocate(); + outputImage->FillBuffer(outputValue); + + OutputImageIteratorType outItr(outputImage, outputImage->GetLargestPossibleRegion()); + vnl_matrix inputValues, correctedInputValues; + anima::DirichletDistribution dirichletDistribution; + std::vector computedOutputValue; + std::vector usefulValues(nbImages, false); + std::vector usefulComponents(nbComponents, false); + double epsValue = std::sqrt(std::numeric_limits::epsilon()); + + while (!outItr.IsAtEnd()) + { + // Discard background voxels or voxels full of zeros + unsigned int nbUsedImages = 0; + std::fill(usefulValues.begin(), usefulValues.end(), false); + if (maskArg.getValue() != "") + { + for (unsigned int i = 0;i < nbImages;++i) + { + if (maskItrs[i].Get()) + { + outputValue = inItrs[i].Get(); + + double sumValue = 0.0; + for (unsigned int j = 0;j < nbComponents;++j) + sumValue += outputValue[j]; + + if (sumValue > epsValue) + { + usefulValues[i] = true; + nbUsedImages++; + } + } + } + } + else + { + for (unsigned int i = 0;i < nbImages;++i) + { + outputValue = inItrs[i].Get(); + + double sumValue = 0.0; + for (unsigned int j = 0;j < nbComponents;++j) + sumValue += outputValue[j]; + + if (sumValue > epsValue) + { + usefulValues[i] = true; + nbUsedImages++; + } + } + } + + if (nbUsedImages == 0) + { + for (unsigned int i = 0;i < nbImages;++i) + { + ++inItrs[i]; + if (maskArg.getValue() != "") + ++maskItrs[i]; + } + ++outItr; + continue; + } + + inputValues.set_size(nbUsedImages, nbComponents); + unsigned int pos = 0; + for (unsigned int i = 0;i < nbImages;++i) + { + if (usefulValues[i]) + { + outputValue = inItrs[i].Get(); + for (unsigned int j = 0;j < nbComponents;++j) + inputValues.put(pos, j, outputValue[j]); + pos++; + } + } + + std::fill(usefulComponents.begin(), usefulComponents.end(), false); + unsigned int nbValidComponents = 0; + for (unsigned int j = 0;j < nbComponents;++j) + { + double meanValue = 0.0; + for (unsigned int i = 0;i < nbUsedImages;++i) + meanValue += inputValues.get(i, j); + meanValue /= static_cast(nbUsedImages); + + double varValue = 0.0; + for (unsigned int i = 0;i < nbUsedImages;++i) + varValue += (inputValues.get(i, j) - meanValue) * (inputValues.get(i, j) - meanValue); + varValue /= (nbUsedImages - 1.0); + + if (varValue > epsValue) + { + usefulComponents[j] = true; + nbValidComponents++; + } + } + + if (nbValidComponents == 0) + { + for (unsigned int i = 0;i < nbImages;++i) + { + ++inItrs[i]; + if (maskArg.getValue() != "") + ++maskItrs[i]; + } + ++outItr; + continue; + } + + correctedInputValues.set_size(nbUsedImages, nbValidComponents); + pos = 0; + for (unsigned int j = 0;j < nbComponents;++j) + { + if (usefulComponents[j]) + { + for (unsigned int i = 0;i < nbUsedImages;++i) + correctedInputValues.put(i, pos, inputValues.get(i, j)); + pos++; + } + } + + dirichletDistribution.Fit(correctedInputValues, "mle"); + computedOutputValue = dirichletDistribution.GetConcentrationParameters(); + + if (!std::isfinite(computedOutputValue[0])) + { + std::cout << inputValues << std::endl; + exit(-1); + } + + outputValue.Fill(1.0); + pos = 0; + for (unsigned int j = 0;j < nbComponents;++j) + { + if (usefulComponents[j]) + { + outputValue[j] = computedOutputValue[pos]; + pos++; + } + } + + outItr.Set(outputValue); + + for (unsigned int i = 0;i < nbImages;++i) + { + ++inItrs[i]; + if (maskArg.getValue() != "") + ++maskItrs[i]; + } + ++outItr; + } + + anima::writeImage (outArg.getValue(), outputImage); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Anima/math-tools/statistics/dirichlet_estimation_test/CMakeLists.txt b/Anima/math-tools/statistics/dirichlet_estimation_test/CMakeLists.txt new file mode 100644 index 000000000..ef813926e --- /dev/null +++ b/Anima/math-tools/statistics/dirichlet_estimation_test/CMakeLists.txt @@ -0,0 +1,36 @@ +if(BUILD_TESTING) + +project(animaDirichletEstimationTest) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + AnimaStatisticalDistributions + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/math-tools/statistics/dirichlet_estimation_test/animaDirichletEstimationTest.cxx b/Anima/math-tools/statistics/dirichlet_estimation_test/animaDirichletEstimationTest.cxx new file mode 100644 index 000000000..12b024e1f --- /dev/null +++ b/Anima/math-tools/statistics/dirichlet_estimation_test/animaDirichletEstimationTest.cxx @@ -0,0 +1,45 @@ +#include +#include +#include +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg aArg("a", "a1", "First concentration parameter.", true, 1, "first concentration", cmd); + TCLAP::ValueArg bArg("b", "a2", "Second concentration parameter.", true, 1, "second concentration", cmd); + TCLAP::ValueArg nSampleArg("n", "nsample", "Sample size.", false, 10000, "sample size", cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + std::vector concentrationParameters(2); + concentrationParameters[0] = aArg.getValue(); + concentrationParameters[1] = bArg.getValue(); + + anima::DirichletDistribution dirichletDistribution; + dirichletDistribution.SetConcentrationParameters(concentrationParameters); + + std::mt19937 generator(1234); + + vnl_matrix sample(nSampleArg.getValue(), 2); + dirichletDistribution.Random(sample, generator); + + dirichletDistribution.Fit(sample, "mle"); + + concentrationParameters = dirichletDistribution.GetConcentrationParameters(); + std::cout << "Concentration parameters: "; + for (unsigned int i = 0;i < sample.cols();++i) + std::cout << concentrationParameters[i] << " "; + std::cout << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Anima/math-tools/statistics/gamma_estimation/animaGammaEstimation.cxx b/Anima/math-tools/statistics/gamma_estimation/animaGammaEstimation.cxx index d2740576a..693eeb189 100644 --- a/Anima/math-tools/statistics/gamma_estimation/animaGammaEstimation.cxx +++ b/Anima/math-tools/statistics/gamma_estimation/animaGammaEstimation.cxx @@ -95,27 +95,37 @@ int main(int argc, char **argv) std::vector inputValues; anima::GammaDistribution gammaDistribution; + double epsValue = std::sqrt(std::numeric_limits::epsilon()); + std::vector usefulValues(nbImages, false); while (!kappaItr.IsAtEnd()) { - inputValues.clear(); - unsigned int nbUsedImages = 0; + std::fill(usefulValues.begin(), usefulValues.end(), false); if (maskArg.getValue() != "") { for (unsigned int i = 0;i < nbImages;++i) { - bool maskValue = maskItrs[i].Get(); - nbUsedImages += maskValue; - if (maskValue) - inputValues.push_back(inItrs[i].Get()); + if (maskItrs[i].Get() && inItrs[i].Get() > epsValue) + { + usefulValues[i] = true; + nbUsedImages++; + } } } else - nbUsedImages = nbImages; + { + for (unsigned int i = 0;i < nbImages;++i) + { + if (inItrs[i].Get() > epsValue) + { + usefulValues[i] = true; + nbUsedImages++; + } + } + } - - if (nbUsedImages == 0) + if (nbUsedImages == 0) { for (unsigned int i = 0;i < nbImages;++i) { @@ -127,6 +137,17 @@ int main(int argc, char **argv) ++thetaItr; continue; } + + inputValues.resize(nbUsedImages); + unsigned int pos = 0; + for (unsigned int i = 0;i < nbImages;++i) + { + if (usefulValues[i]) + { + inputValues[pos] = inItrs[i].Get(); + ++pos; + } + } gammaDistribution.Fit(inputValues, methodArg.getValue()); diff --git a/Anima/math-tools/statistics/watson_estimation_test/CMakeLists.txt b/Anima/math-tools/statistics/watson_estimation_test/CMakeLists.txt new file mode 100644 index 000000000..55036e559 --- /dev/null +++ b/Anima/math-tools/statistics/watson_estimation_test/CMakeLists.txt @@ -0,0 +1,36 @@ +if(BUILD_TESTING) + +project(animaWatsonEstimationTest) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + AnimaStatisticalDistributions + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/math-tools/statistics/watson_estimation_test/animaWatsonEstimationTest.cxx b/Anima/math-tools/statistics/watson_estimation_test/animaWatsonEstimationTest.cxx new file mode 100644 index 000000000..b14da4d6b --- /dev/null +++ b/Anima/math-tools/statistics/watson_estimation_test/animaWatsonEstimationTest.cxx @@ -0,0 +1,52 @@ +#include +#include +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg kappaArg("k", "kappa", "Shape parameter.", true, 1, "shape parameter", cmd); + TCLAP::ValueArg nSampleArg("n", "nsample", "Sample size.", false, 10000, "sample size", cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + double kappa = kappaArg.getValue(); + + anima::WatsonDistribution watsonDistribution; + + using VectorType = itk::Vector; + VectorType meanAxis; + meanAxis.Fill(1.0 / std::sqrt(3.0)); + double normValue = meanAxis.GetNorm(); + meanAxis /= normValue; + watsonDistribution.SetMeanAxis(meanAxis); + watsonDistribution.SetConcentrationParameter(kappa); + + std::mt19937 generator(1234); + + std::vector sample(nSampleArg.getValue()); + watsonDistribution.Random(sample, generator); + + watsonDistribution.Fit(sample, ""); + + std::cout << "----- Mean Axis -----" << std::endl; + std::cout << "- Expected : " << meanAxis << std::endl; + std::cout << "- Estimated: " << watsonDistribution.GetMeanAxis() << std::endl; + + std::cout << std::endl; + + std::cout << "--- Concentration ---" << std::endl; + std::cout << "- Expected : " << kappa << std::endl; + std::cout << "- Estimated: " << watsonDistribution.GetConcentrationParameter() << std::endl; + + return EXIT_SUCCESS; +}