From d1a361a0e7d412ff1a3f0e44b8d745685e0393fa Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Wed, 15 Nov 2023 23:04:15 +0100 Subject: [PATCH] Clean up TOD estimator. --- .../odf/tod_estimator/CMakeLists.txt | 45 +- .../odf/tod_estimator/animaTODEstimator.cxx | 71 +- .../animaTODEstimatorImageFilter.h | 176 ++- .../animaTODEstimatorImageFilter.hxx | 1114 ++++++----------- 4 files changed, 547 insertions(+), 859 deletions(-) diff --git a/Anima/diffusion/odf/tod_estimator/CMakeLists.txt b/Anima/diffusion/odf/tod_estimator/CMakeLists.txt index 5b8e31076..3effe2ce5 100644 --- a/Anima/diffusion/odf/tod_estimator/CMakeLists.txt +++ b/Anima/diffusion/odf/tod_estimator/CMakeLists.txt @@ -1,39 +1,30 @@ if(BUILD_TOOLS AND USE_NLOPT) -project(animaTODEstimator) + project(animaTODEstimator) -## ############################################################################# -## List Sources -## ############################################################################# + # ############################################################################ + # List Sources + # ############################################################################ -list_source_files(${PROJECT_NAME} - ${CMAKE_CURRENT_SOURCE_DIR} - ) + list_source_files(${PROJECT_NAME} ${CMAKE_CURRENT_SOURCE_DIR}) + # ############################################################################ + # add executable + # ############################################################################ -## ############################################################################# -## add executable -## ############################################################################# + add_executable(${PROJECT_NAME} ${${PROJECT_NAME}_CFILES}) -add_executable(${PROJECT_NAME} - ${${PROJECT_NAME}_CFILES} - ) + # ############################################################################ + # Link + # ############################################################################ -## ############################################################################# -## Link -## ############################################################################# + target_link_libraries(${PROJECT_NAME} ${ITKIO_LIBRARIES} AnimaDataIO + AnimaSHTools) -target_link_libraries(${PROJECT_NAME} - ${VTK_LIBRARIES} - ${ITKIO_LIBRARIES} - AnimaSHTools - AnimaDataIO - ) + # ############################################################################ + # install + # ############################################################################ -## ############################################################################# -## install -## ############################################################################# - -set_exe_install_rules(${PROJECT_NAME}) + set_exe_install_rules(${PROJECT_NAME}) endif() diff --git a/Anima/diffusion/odf/tod_estimator/animaTODEstimator.cxx b/Anima/diffusion/odf/tod_estimator/animaTODEstimator.cxx index 2f991ddb7..4adcb43b1 100644 --- a/Anima/diffusion/odf/tod_estimator/animaTODEstimator.cxx +++ b/Anima/diffusion/odf/tod_estimator/animaTODEstimator.cxx @@ -1,54 +1,61 @@ +#include #include + #include -#include -#include +#include -void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData) +void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData) { - itk::ProcessObject * processObject = (itk::ProcessObject*) caller; - std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"<GetProgress() * 100) << "%" << std::flush; } int main(int argc, char **argv) { - TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION); - - TCLAP::ValueArg inArg("i","input","Input tractography image (.vtk or .vtp)",true,"","input tractography image",cmd); - TCLAP::ValueArg resArg("o","outputfile","Result TOD image",true,"","result TOD image",cmd); - TCLAP::ValueArg refArg("g","geometry","Output image geometry",true,"","output geometry",cmd); - - TCLAP::SwitchArg normArg("N", "Normalize", "Normalize TOD", cmd, false); - - TCLAP::ValueArg orderArg("k","order","Order of spherical harmonics basis (default 4)",false,4,"Order of SH basis",cmd); - - TCLAP::ValueArg nbpArg("p","numberofthreads","Number of threads to run on (default: all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd); + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg inArg( + "i", "input-file", + "A string specifying the name of a file storing the input tractography image. Supported formats are `.vtk`, `.vtp` or `.fds`.", + true, "", "input tractography image", cmd); + TCLAP::ValueArg outArg( + "o", "output-file", + "A string specifying the name of a file storing the output TOD image.", + true, "", "output TOD image", cmd); + TCLAP::ValueArg refArg( + "g", "geometry-file", + "A string specifying the name of a file storing the reference geometry image.", + true, "", "reference geometry image", cmd); + + TCLAP::SwitchArg normArg( + "N", "normalize-tod", + "A switch to turn on TOD normalization.", + cmd, false); + + TCLAP::ValueArg nbpArg( + "T", "nb-threads", + "An integer value specifying the number of threads to run on (default: all cores).", + false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(), "number of threads", cmd); try { - cmd.parse(argc,argv); + cmd.parse(argc, argv); } - catch (TCLAP::ArgException& e) + catch (TCLAP::ArgException &e) { std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; return EXIT_FAILURE; } - typedef anima::TODEstimatorImageFilter FilterType; + using FilterType = anima::TODEstimatorImageFilter; + using InputImageType = FilterType::InputImageType; FilterType::Pointer mainFilter = FilterType::New(); - -// if (orderArg.getValue() % 2 == 0) -// mainFilter->SetLOrder(orderArg.getValue()); -// else -// mainFilter->SetLOrder(orderArg.getValue() - 1); - typedef FilterType::InputImageType InputImageType; mainFilter->SetInput(anima::readImage(refArg.getValue())); - - mainFilter->SetLOrder(orderArg.getValue()); mainFilter->SetInputFileName(inArg.getValue()); - mainFilter->SetRefFileName(refArg.getValue()); - mainFilter->SetNormalize(normArg.getValue()); + mainFilter->SetReferenceFileName(refArg.getValue()); + mainFilter->SetUseNormalization(normArg.getValue()); mainFilter->SetNumberOfWorkUnits(nbpArg.getValue()); itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New(); @@ -57,14 +64,12 @@ int main(int argc, char **argv) itk::TimeProbe tmpTime; tmpTime.Start(); - mainFilter->Update(); - tmpTime.Stop(); - std::cout << std::endl << "Execution Time: " << tmpTime.GetTotal() << std::endl; + std::cout << "\nExecution Time: " << tmpTime.GetTotal() << "s" << std::endl; - anima::writeImage (resArg.getValue(),mainFilter->GetOutput()); + anima::writeImage(outArg.getValue(), mainFilter->GetOutput()); return EXIT_SUCCESS; } diff --git a/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.h b/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.h index 29abdc6f4..d84896712 100644 --- a/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.h +++ b/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.h @@ -1,138 +1,112 @@ #pragma once -#include -#include -#include -#include -#include - #include #include -#include -#include #include -#include namespace anima { -//template -class TODEstimatorImageFilter : -public anima::NumberedThreadImageToImageFilter < itk::Image, itk::VectorImage > -{ -public: - - typedef TODEstimatorImageFilter Self; -// typedef itk::Vector pointType; - typedef itk::Point PointType; -// typedef PointType DirType; - typedef itk::Vector DirType; - typedef std::vector DirVectorType; - typedef std::vector FiberType; - - typedef double MathScalarType; - - typedef itk::Matrix Matrix3DType; - typedef itk::Vector Vector3DType; - typedef itk::VariableLengthVector VectorType; - - typedef anima::ODFSphericalHarmonicBasis baseSH; - typedef std::complex complexType; - - typedef itk::VectorImage TOutputImage; - typedef itk::Image TRefImage; - - typedef anima::NumberedThreadImageToImageFilter Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - - itkNewMacro(Self) - - itkTypeMacro(TODEstimatorImageFilter, anima::NumberedThreadImageToImageFilter); - - typedef typename TOutputImage::Pointer OutputImagePointer; - - typedef typename Superclass::InputImageRegionType InputImageRegionType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - - itkSetMacro(InputFileName,std::string); - itkSetMacro(RefFileName,std::string); - itkSetMacro(LOrder,unsigned int); - itkSetMacro(Normalize, bool); - - -protected: - TODEstimatorImageFilter() - : Superclass() + template + class TODEstimatorImageFilter : public anima::NumberedThreadImageToImageFilter, itk::VectorImage> { - } + public: + /** Standard class typedefs. */ + using Self = TODEstimatorImageFilter; + using InputImageType = itk::Image; + using OutputImageType = itk::VectorImage; + using ReferenceImageType = itk::Image; + using Superclass = anima::NumberedThreadImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; - virtual ~TODEstimatorImageFilter() - { + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods) */ + itkTypeMacro(TODEstimatorImageFilter, anima::NumberedThreadImageToImageFilter); - } + /** Superclass typedefs. */ + using InputImageRegionType = typename Superclass::InputImageRegionType; + using OutputImageRegionType = typename Superclass::OutputImageRegionType; -// void GenerateData() ITK_OVERRIDE; + using InputImagePointerType = typename InputImageType::Pointer; + using OutputImagePointerType = typename OutputImageType::Pointer; + using ReferenceImagePointerType = typename ReferenceImageType::Pointer; + using InputImagePixelType = typename InputImageType::PixelType; + using OutputImagePixelType = typename OutputImageType::PixelType; + using ReferenceImagePixelType = ReferenceImageType::PixelType; - void BeforeThreadedGenerateData() ITK_OVERRIDE; - void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE; - void AfterThreadedGenerateData() ITK_OVERRIDE; + using PointType = itk::Point; + using DirType = itk::Vector; + using DirVectorType = std::vector; + using FiberType = std::vector; + using Matrix3DType = itk::Matrix; + using Vector3DType = itk::Vector; + using MatrixType = vnl_matrix; + using BasisType = anima::ODFSphericalHarmonicBasis; + using ComplexType = std::complex; + itkSetMacro(InputFileName, std::string); + itkSetMacro(ReferenceFileName, std::string); + itkSetMacro(LOrder, unsigned int); + itkSetMacro(UseNormalization, bool); - FiberType readFiber(vtkIdType numberOfPoints, const vtkIdType *indices, vtkPoints *points); - PointType getCenterVoxel(int index, FiberType &fiber); - DirType getFiberDirection(int index, FiberType &fiber); - void getSHCoefs(DirType dir, VectorType &resSH, baseSH &basis); - void ComputeCoefs(); - void processFiber(FiberType &fiber, baseSH &basis); + protected: + TODEstimatorImageFilter() {} - void getMainDirections(DirVectorType inDirs, DirVectorType &mainDirs); - double getEuclideanDistance(DirType dir1, DirType dir2); - DirType getNewClusterAverage(int numCluster, DirVectorType &dirs, std::vector &cluster); + virtual ~TODEstimatorImageFilter() {} - void getSHCoef(DirType dir, VectorType &coefs); + void BeforeThreadedGenerateData() ITK_OVERRIDE; + void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE; - void precomputeSH(); - void discretizeODF(VectorType ODFCoefs, std::vector &ODFDiscret); - VectorType getSquareRootODF(std::vector ODFDiscret); - VectorType getSquareODF(std::vector ODFDiscret); - void getAverageCoefs(std::vector &vecCoefs, VectorType &avgCoef); - - void averageODFs(std::vector &vecCoefs, VectorType &resOdf); + FiberType ReadFiber(vtkIdType numberOfPoints, const vtkIdType *indices, vtkPoints *points); + PointType GetCenterVoxel(int index, FiberType &fiber); + DirType GetFiberDirection(int index, FiberType &fiber); + void GetSHCoefs(DirType dir, OutputImagePixelType &resSH, BasisType &basis); + void ComputeCoefs(); + void ProcessFiber(FiberType &fiber, BasisType &basis); - vnl_matrix GetRotationMatrix(DirType dir1, DirType dir2); + void GetMainDirections(DirVectorType inDirs, DirVectorType &mainDirs); + double GetEuclideanDistance(DirType dir1, DirType dir2); + DirType GetNewClusterAverage(int numCluster, DirVectorType &dirs, std::vector &cluster); -// void GenerateOutputInformation() ITK_OVERRIDE; + void GetSHCoef(DirType dir, OutputImagePixelType &coefs); -private: - ITK_DISALLOW_COPY_AND_ASSIGN(TODEstimatorImageFilter); + void PrecomputeSH(); + void DiscretizeODF(OutputImagePixelType ODFCoefs, std::vector &ODFDiscret); + OutputImagePixelType GetSquareRootODF(std::vector ODFDiscret); + OutputImagePixelType GetSquareODF(std::vector ODFDiscret); + void GetAverageCoefs(std::vector &vecCoefs, OutputImagePixelType &avgCoef); - std::string m_InputFileName; - std::string m_RefFileName; + void AverageODFs(std::vector &vecCoefs, OutputImagePixelType &resOdf); - DirType m_CstDir; + MatrixType GetRotationMatrix(DirType dir1, DirType dir2); - unsigned int m_LOrder; - int m_VectorLength; + private: + ITK_DISALLOW_COPY_AND_ASSIGN(TODEstimatorImageFilter); - double m_NbSample; + std::string m_InputFileName; + std::string m_ReferenceFileName; - bool m_Normalize; + DirType m_CstDir; + unsigned int m_LOrder; + int m_VectorLength; - VectorType m_GaussCoefs; + double m_NbSample; - std::vector> m_SphereSampl; - vnl_matrix m_SpherHarm; + bool m_UseNormalization; - anima::ODFSphericalHarmonicBasis *m_ODFSHBasis; - std::vector> m_ImgDir; + OutputImagePixelType m_GaussCoefs; - itk::Image::Pointer test; + std::vector> m_SphereSampl; + MatrixType m_SpherHarm; -}; + anima::ODFSphericalHarmonicBasis *m_ODFSHBasis; + std::vector m_ImgDir; + }; } // end namespace anima #include "animaTODEstimatorImageFilter.hxx" diff --git a/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.hxx b/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.hxx index e83a81817..e06d13200 100644 --- a/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.hxx +++ b/Anima/diffusion/odf/tod_estimator/animaTODEstimatorImageFilter.hxx @@ -1,850 +1,568 @@ #pragma once #include "animaTODEstimatorImageFilter.h" - -#include -#include -#include -#include -#include -#include - -#include +#include +#include #include -#include #include -#include -#include #include -#include "animaLogExpMapsUnitSphere.h" + +#include #include namespace anima { + template + void + TODEstimatorImageFilter::BeforeThreadedGenerateData() + { + Superclass::BeforeThreadedGenerateData(); + + ReferenceImagePointerType refImg = anima::readImage(m_ReferenceFileName); + m_LOrder = 8; + unsigned int vectorLength = (m_LOrder + 1) * (m_LOrder + 2) / 2; + m_VectorLength = vectorLength; + OutputImageType *output = this->GetOutput(); + output->SetVectorLength(vectorLength); + output->SetRegions(refImg->GetLargestPossibleRegion()); + output->SetSpacing(refImg->GetSpacing()); + output->SetDirection(refImg->GetDirection()); + output->SetOrigin(refImg->GetOrigin()); + output->Allocate(); + + OutputImagePixelType tmp; + tmp.SetSize((m_LOrder + 1) * (m_LOrder + 2) / 2); + tmp.Fill(0); + output->FillBuffer(tmp); -void TODEstimatorImageFilter::BeforeThreadedGenerateData() -{ - Superclass::BeforeThreadedGenerateData(); - - TRefImage::Pointer refImg = anima::readImage(m_RefFileName); - - unsigned int vectorLength = (m_LOrder + 1)*(m_LOrder + 2)/2; - m_VectorLength = vectorLength; - TOutputImage *output = this->GetOutput(); - output->SetVectorLength(vectorLength); - output->SetRegions(refImg->GetLargestPossibleRegion()); - output->SetSpacing(refImg->GetSpacing()); - output->SetDirection(refImg->GetDirection()); - output->SetOrigin(refImg->GetOrigin()); - output->Allocate(); - - VectorType tmp; - tmp.SetSize((m_LOrder+1)*(m_LOrder+2)/2); - tmp.Fill(0); - output->FillBuffer(tmp); - - test = itk::Image::New(); - test->SetRegions(refImg->GetLargestPossibleRegion()); - test->SetSpacing(refImg->GetSpacing()); - test->SetDirection(refImg->GetDirection()); - test->SetOrigin(refImg->GetOrigin()); - test->Allocate(); - - m_CstDir[0] = 0; - m_CstDir[1] = 0; - m_CstDir[2] = 1; - - m_NbSample = 200; - anima::GetSphereEvenSampling(m_SphereSampl, m_NbSample); - - m_SpherHarm.set_size(m_NbSample, vectorLength); - - m_ODFSHBasis = new anima::ODFSphericalHarmonicBasis(m_LOrder); - this->precomputeSH(); - - const int Xmax = output->GetLargestPossibleRegion().GetSize()[0]; - const int Ymax = output->GetLargestPossibleRegion().GetSize()[1]; - const int Zmax = output->GetLargestPossibleRegion().GetSize()[2]; - - const int Xmin = output->GetLargestPossibleRegion().GetIndex()[0]; - const int Ymin = output->GetLargestPossibleRegion().GetIndex()[1]; - const int Zmin = output->GetLargestPossibleRegion().GetIndex()[2]; + m_CstDir[0] = 0; + m_CstDir[1] = 0; + m_CstDir[2] = 1; - const int dirX = output->GetDirection()[0][0]; - const int dirY = output->GetDirection()[1][1]; - const int dirZ = output->GetDirection()[2][2]; + m_NbSample = 200; + anima::GetSphereEvenSampling(m_SphereSampl, m_NbSample); - anima::ShapesReader trackReader; - trackReader.SetFileName(m_InputFileName); - trackReader.Update(); + m_SpherHarm.set_size(m_NbSample, vectorLength); - this->ComputeCoefs(); + m_ODFSHBasis = new anima::ODFSphericalHarmonicBasis(m_LOrder); + this->PrecomputeSH(); + this->ComputeCoefs(); - m_ImgDir.resize(Xmax * Ymax * Zmax); + const int Xmax = output->GetLargestPossibleRegion().GetSize()[0]; + const int Ymax = output->GetLargestPossibleRegion().GetSize()[1]; + const int Zmax = output->GetLargestPossibleRegion().GetSize()[2]; + m_ImgDir.resize(Xmax * Ymax * Zmax); - vtkSmartPointer tracks = trackReader.GetOutput(); - vtkSmartPointer cells = tracks->GetLines(); - vtkSmartPointer points = tracks->GetPoints(); + anima::ShapesReader trackReader; + trackReader.SetFileName(m_InputFileName); + trackReader.Update(); - const vtkIdType *indices; - vtkIdType numberOfPoints; + vtkSmartPointer tracks = trackReader.GetOutput(); + vtkSmartPointer cells = tracks->GetLines(); + vtkSmartPointer points = tracks->GetPoints(); - float nbLines = tracks->GetNumberOfLines(); - float lineCount = 0; + const vtkIdType *indices; + vtkIdType numberOfPoints; - std::cout << "Number of fibers + test : " << nbLines << std::endl; + float nbLines = tracks->GetNumberOfLines(); + float lineCount = 0; - anima::ODFSphericalHarmonicBasis basis(m_LOrder); + std::cout << "Number of fibers + test : " << nbLines << std::endl; - auto iter = vtk::TakeSmartPointer(cells->NewIterator()); - for (iter->GoToFirstCell(); !iter->IsDoneWithTraversal(); iter->GoToNextCell()) - { - iter->GetCurrentCell(numberOfPoints, indices); - // std::cout << "Progress : " << (int) (((float)lineCount / nbLines)*100) << "%\r"; - // fflush (stdout); - lineCount++; + anima::ODFSphericalHarmonicBasis basis(m_LOrder); - FiberType fiber = readFiber(numberOfPoints, indices, points); - int nbPoints = fiber.size(); + auto iter = vtk::TakeSmartPointer(cells->NewIterator()); + for (iter->GoToFirstCell(); !iter->IsDoneWithTraversal(); iter->GoToNextCell()) + { + iter->GetCurrentCell(numberOfPoints, indices); + lineCount++; + FiberType fiber = ReadFiber(numberOfPoints, indices, points); + int nbPoints = fiber.size(); - for(int i = 0; i < nbPoints; i++) - { - DirType dir = getFiberDirection(i, fiber); - // OutputImagePointer::DirectionType refDir = output->GetDirection(); - - // if(output->GetDirection()[0][0] > 0) - // dir[0] *= -1; - // if(output->GetDirection()[1][1] < 0) - // dir[1] *= -1; - // if(output->GetDirection()[2][2] < 0) - // dir[2] *= -1; - // dir[0] *= -1; - // dir[1] *= -1; - // dir[2] *= -1; - - // anima::Normalize(dir, dir); - PointType point = getCenterVoxel(i, fiber); - itk::Index<3> index, index2; - - output->TransformPhysicalPointToIndex(point, index); - - index[0] *= dirX; - index[1] *= dirY; - index[2] *= dirZ; - - index2[0] = std::floor(point[0]); - index2[1] = std::floor(point[0]); - index2[2] = std::floor(point[0]); - - m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])].push_back(dir); - - // // itk::ImageRegionIterator > testItr(test, test->GetLargestPossibleRegion()); - // // while(!testItr.IsAtEnd()) - // // { - // // itk::Image::IndexType index; - // // index = testItr.GetIndex(); - // // testItr.Set(imgDir[index[0] + Xmax * (index[1] + Ymax * index[2])].size()); - // // ++testItr; - // // } - // // anima::writeImage>("./test.nii.gz", test); + for (int i = 0; i < nbPoints; i++) + { + DirType dir = GetFiberDirection(i, fiber); + PointType point = GetCenterVoxel(i, fiber); + itk::Index<3> index; + output->TransformPhysicalPointToIndex(point, index); + m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])].push_back(dir); + } } } -} - -void TODEstimatorImageFilter::ComputeCoefs() -{ - m_GaussCoefs.SetSize(91); - m_GaussCoefs.Fill(0.0); - - double l1 = 1.71e-4; - double l2 = 2e-5; - double d = l1 - l2; - double spi = sqrt(M_PI); - double pi = M_PI; - double at = M_PI+2*atan((d-l2)/(2*sqrt(d*l2))); - - double c0 = 1/(2*sqrt(pi)); - double c1 = - sqrt(5) * (3*sqrt(l2)*(d+l2)*(pi + atan((d-l2)/(2*sqrt(l2*d)))) - 4*sqrt(d)*(2*d+3*l2)) / (16*sqrt(pi)*pow(d,1.5)); - double c2 = 3 * (4*sqrt(d)*(16*pow(d,2)+115*d*l2+105*pow(l2,2)) - 15*sqrt(l2)*(3*pow(d,2)+10*d*l2+7*pow(l2,2))*(pi+2*atan((d-l2)/(2*sqrt(d*l2))))) / (128*spi*pow(d,2.5)); - double c3 = -105*sqrt(13)*(at*sqrt(l2)*(5*pow(d,3)+35*l2*pow(d,2)+63*pow(l2,2)*d+33*pow(l2,3)) - 4*sqrt(d)*(33*pow(l2,3)+52*d*pow(l2,2)+103*pow(d,2)*l2/5+128*pow(d,3)/105)) / (1024*spi*pow(d,3.5)); - double c4 = -315*sqrt(17) * (at*sqrt(l2)*(35*pow(d,4)+420*pow(d,3)*l2+1386*pow(d,2)*pow(l2,2)+1716*pow(l2,3)*d+715*pow(l2,4)) - 4*sqrt(d)*(715*pow(l2,4)+4433*pow(l2,3)*d/3+957*pow(l2,2)*pow(d,2)+6967*l2*pow(d,3)/35+2048*pow(d,4)/315)) / (16384 * spi * pow(d,4.5)); - double c5 = -24255*sqrt(21) * (at * sqrt(l2)*(9*pow(d,5)+165*pow(d,4)*l2+858*pow(d,3)*pow(l2,2)+12870*pow(d,2)*pow(l2,3)/7+12155*d*pow(l2,4)/7+4199*pow(l2,5)/7) - 4*sqrt(d)*(4199*pow(l2,5)/7+32266*pow(l2,4)*d/21+20691*pow(l2,3)*pow(d,2)/15+24830*pow(l2,2)*pow(d,3)/49+28799*l2*pow(d,4)/441+32768*pow(d,5)/24255)) / (262144*pow(d,5.5)*spi); - - - // m_GaussCoefs[0] = this->m_ODFSHBasis->getNthSHValueAtPosition(0, 0, 0, 0); - // std::vector Dir(3), sphDir(2); - // Dir[0] = 0; - // Dir[1] = 0; - // Dir[2] = 1; - // anima::TransformCartesianToSphericalCoordinates(Dir, sphDir); - // int k = 0; - // for(int l = 0; l <= m_LOrder; l+=2) - // { - // for(int m = -l; m <= l; m++) - // { - // m_GaussCoefs[k] = this->m_ODFSHBasis->getNthSHValueAtPosition(l, m, sphDir[0], sphDir[1]); - // k++; - // } - // } - - m_GaussCoefs[0] = c0; - m_GaussCoefs[3] = c1; - m_GaussCoefs[10] = c2; - m_GaussCoefs[21] = c3; - m_GaussCoefs[36] = c4; - m_GaussCoefs[55] = c5; -} - -void TODEstimatorImageFilter::DynamicThreadedGenerateData(const TODEstimatorImageFilter::OutputImageRegionType &outputRegionForThread) -{ - TOutputImage *output = this->GetOutput(); - itk::ImageRegionIterator outItr(output, outputRegionForThread); - typename TOutputImage::IndexType index, index_out; - - VectorType nullVector; - nullVector.SetSize((m_LOrder+1)*(m_LOrder+2)/2); - nullVector.Fill(0); - - const int Xmax = output->GetLargestPossibleRegion().GetSize()[0]; - const int Ymax = output->GetLargestPossibleRegion().GetSize()[1]; - const int Zmax = output->GetLargestPossibleRegion().GetSize()[2]; - const int dirX = output->GetDirection()[0][0]; - const int dirY = output->GetDirection()[1][1]; - const int dirZ = output->GetDirection()[2][2]; - - while(!outItr.IsAtEnd()) + template + void + TODEstimatorImageFilter::ComputeCoefs() { + m_GaussCoefs.SetSize(91); + m_GaussCoefs.Fill(0.0); + + double l1 = 1.71e-4; + double l2 = 2e-5; + double d = l1 - l2; + double spi = sqrt(M_PI); + double pi = M_PI; + double at = M_PI + 2 * atan((d - l2) / (2 * sqrt(d * l2))); + + double c0 = 1 / (2 * sqrt(pi)); + double c1 = -sqrt(5) * (3 * sqrt(l2) * (d + l2) * (pi + atan((d - l2) / (2 * sqrt(l2 * d)))) - 4 * sqrt(d) * (2 * d + 3 * l2)) / (16 * sqrt(pi) * pow(d, 1.5)); + double c2 = 3 * (4 * sqrt(d) * (16 * pow(d, 2) + 115 * d * l2 + 105 * pow(l2, 2)) - 15 * sqrt(l2) * (3 * pow(d, 2) + 10 * d * l2 + 7 * pow(l2, 2)) * (pi + 2 * atan((d - l2) / (2 * sqrt(d * l2))))) / (128 * spi * pow(d, 2.5)); + double c3 = -105 * sqrt(13) * (at * sqrt(l2) * (5 * pow(d, 3) + 35 * l2 * pow(d, 2) + 63 * pow(l2, 2) * d + 33 * pow(l2, 3)) - 4 * sqrt(d) * (33 * pow(l2, 3) + 52 * d * pow(l2, 2) + 103 * pow(d, 2) * l2 / 5 + 128 * pow(d, 3) / 105)) / (1024 * spi * pow(d, 3.5)); + double c4 = -315 * sqrt(17) * (at * sqrt(l2) * (35 * pow(d, 4) + 420 * pow(d, 3) * l2 + 1386 * pow(d, 2) * pow(l2, 2) + 1716 * pow(l2, 3) * d + 715 * pow(l2, 4)) - 4 * sqrt(d) * (715 * pow(l2, 4) + 4433 * pow(l2, 3) * d / 3 + 957 * pow(l2, 2) * pow(d, 2) + 6967 * l2 * pow(d, 3) / 35 + 2048 * pow(d, 4) / 315)) / (16384 * spi * pow(d, 4.5)); + double c5 = -24255 * sqrt(21) * (at * sqrt(l2) * (9 * pow(d, 5) + 165 * pow(d, 4) * l2 + 858 * pow(d, 3) * pow(l2, 2) + 12870 * pow(d, 2) * pow(l2, 3) / 7 + 12155 * d * pow(l2, 4) / 7 + 4199 * pow(l2, 5) / 7) - 4 * sqrt(d) * (4199 * pow(l2, 5) / 7 + 32266 * pow(l2, 4) * d / 21 + 20691 * pow(l2, 3) * pow(d, 2) / 15 + 24830 * pow(l2, 2) * pow(d, 3) / 49 + 28799 * l2 * pow(d, 4) / 441 + 32768 * pow(d, 5) / 24255)) / (262144 * pow(d, 5.5) * spi); + + m_GaussCoefs[0] = c0; + m_GaussCoefs[3] = c1; + m_GaussCoefs[10] = c2; + m_GaussCoefs[21] = c3; + m_GaussCoefs[36] = c4; + m_GaussCoefs[55] = c5; + } + template + void + TODEstimatorImageFilter::DynamicThreadedGenerateData(const TODEstimatorImageFilter::OutputImageRegionType &outputRegionForThread) + { + OutputImageType *output = this->GetOutput(); + itk::ImageRegionIterator outItr(output, outputRegionForThread); + typename OutputImageType::IndexType index, index_out; - index = outItr.GetIndex(); + OutputImagePixelType nullVector; + nullVector.SetSize((m_LOrder + 1) * (m_LOrder + 2) / 2); + nullVector.Fill(0); - // index[0] *= -1; - // index[1] *= -1; - // index[2] *= -1; + const int Xmax = output->GetLargestPossibleRegion().GetSize()[0]; + const int Ymax = output->GetLargestPossibleRegion().GetSize()[1]; + const int Zmax = output->GetLargestPossibleRegion().GetSize()[2]; + while (!outItr.IsAtEnd()) + { - // index[0] *= dirX; - // index[1] *= dirY; - // index[2] *= dirZ; + index = outItr.GetIndex(); - std::vector vecCoefs; - unsigned int tmpSize = m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])].size(); - if (tmpSize != 0) - { - DirVectorType mainDirs; - DirVectorType localDir(tmpSize); - localDir = m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])]; - this->getMainDirections(localDir, mainDirs); - - // std::cout << "----------------------------------------" << std::endl << std::endl; - // for(int i = 0; i < tmpSize; i++) - // std::cout << localDir[i] << ", " << std::endl; - // std::cout << std::endl << std::endl; - // for(int i = 0; i < mainDirs.size(); i++) - // std::cout << mainDirs[i] << ", " << std::endl; - - // std::cout << std::endl << std::endl; - - VectorType tmpCoefs(m_VectorLength); - tmpCoefs.Fill(0); - for(int i = 0; i < mainDirs.size(); i++) + std::vector vecCoefs; + int tmpSize = m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])].size(); + if (tmpSize != 0) { - if(!isnan(mainDirs[i][0])) + DirVectorType mainDirs; + DirVectorType localDir(tmpSize); + localDir = m_ImgDir[index[0] + Xmax * (index[1] + Ymax * index[2])]; + this->GetMainDirections(localDir, mainDirs); + + OutputImagePixelType tmpCoefs(m_VectorLength); + tmpCoefs.Fill(0); + for (int i = 0; i < mainDirs.size(); i++) { - this->getSHCoef(mainDirs[i], tmpCoefs); - // std::vector Dir(3), sphDir(2); - // for(int k = 0; k < 3; k++) - // Dir[k] = mainDirs[i][k]; - // anima::TransformCartesianToSphericalCoordinates(Dir, sphDir); - // int j = 0; - // for(int l = 0; l <= m_LOrder; l+=2) - // { - // for(int m = -l; m <= l; m++) - // { - // tmpCoefs[j] = this->m_ODFSHBasis->getNthSHValueAtPosition(l, m, sphDir[0], sphDir[1]); - // j++; - // } - // } - - vecCoefs.push_back(tmpCoefs); + if (!isnan(mainDirs[i][0])) + { + this->GetSHCoef(mainDirs[i], tmpCoefs); + + vecCoefs.push_back(tmpCoefs); + } } - } - VectorType resCoefs(m_VectorLength); - resCoefs.Fill(0); - if(vecCoefs.size() != 0) - { - // this->averageODFs(vecCoefs, resCoefs); - // for(int k = 0; k < vecCoefs.size(); k++) - // { - // for(int i = 0; i < m_VectorLength; i++) - // resCoefs[i] += (vecCoefs[k][i] / vecCoefs.size()); - // } - for(int k = 0; k < vecCoefs.size(); k++) + OutputImagePixelType resCoefs(m_VectorLength); + resCoefs.Fill(0); + if (vecCoefs.size() != 0) { - for(int i = 0; i < m_VectorLength; i++) - resCoefs[i] += vecCoefs[k][i]; + resCoefs = vecCoefs[0]; + this->AverageODFs(vecCoefs, resCoefs); + if (m_UseNormalization) + { + for (int i = 0; i < m_VectorLength; i++) + resCoefs[i] /= vecCoefs.size(); + } + outItr.Set(resCoefs); } - if(m_Normalize) + + else + outItr.Set(nullVector); + + std::vector testSh(m_VectorLength); + for (int i = 0; i < m_VectorLength; i++) { - for(int i = 0; i < m_VectorLength; i++) - resCoefs[i] /= vecCoefs.size(); + testSh[i] = resCoefs[i]; } - outItr.Set(resCoefs); + localDir.clear(); + localDir.shrink_to_fit(); } - else outItr.Set(nullVector); - std::vector testSh(m_VectorLength); - for(int i = 0; i < m_VectorLength; i++) - { - testSh[i] = resCoefs[i]; - } - localDir.clear(); - localDir.shrink_to_fit(); - } - else - outItr.Set(nullVector); + ++outItr; - ++outItr; - - this->IncrementNumberOfProcessedPoints(); + this->IncrementNumberOfProcessedPoints(); + } } -} - -void TODEstimatorImageFilter::AfterThreadedGenerateData() -{ - -} + template + void + TODEstimatorImageFilter::GetMainDirections(DirVectorType inDirs, DirVectorType &mainDirs) + { -//void TODEstimatorImageFilter::GenerateData() -//{ + if (inDirs.size() == 0) + { + mainDirs.resize(0); + return; + } + const int numIt = 200; + int numDirs = inDirs.size(), crit = 0, numClusters = 0; + if (numDirs < 4) + { + numClusters = numDirs; + mainDirs.resize(numClusters); + for (int i = 0; i < numClusters; i++) + { + mainDirs[i] = inDirs[i]; + } + return; + } + else + numClusters = 4; + DirVectorType clustersAverage(numClusters); + std::vector clusters(numDirs); + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<> distrib(0, numDirs); + for (int i = 0; i < numClusters; i++) + { + int k = distrib(gen); + clustersAverage[i] = inDirs[k]; + } + int It = 0; + bool stop = false; + double dist = 0; + while (It < numIt && !stop) + { + std::vector clustersOld(numDirs); + for (int i = 0; i < numDirs; i++) + clustersOld[i] = clusters[i]; -//} + for (int i = 0; i < numDirs; i++) + { + double minDist = 1e10; + for (int j = 0; j < numClusters; j++) + { + dist = std::sqrt((clustersAverage[j][0] - inDirs[i][0]) * (clustersAverage[j][0] - inDirs[i][0]) + (clustersAverage[j][1] - inDirs[i][1]) * (clustersAverage[j][1] - inDirs[i][1]) + (clustersAverage[j][2] - inDirs[i][2]) * (clustersAverage[j][2] - inDirs[i][2])); + if (dist < minDist) + { + minDist = dist; + clusters[i] = j; + } + } + } + crit = 0; + for (int i = 0; i < numClusters; i++) + { + clustersAverage[i] = GetNewClusterAverage(i, inDirs, clusters); + } + for (int i = 0; i < numDirs; i++) + crit += clusters[i] - clustersOld[i]; + if (crit == 0) + stop = true; -void TODEstimatorImageFilter::getMainDirections(DirVectorType inDirs, DirVectorType &mainDirs) -{ + ++It; + } - if(inDirs.size() == 0) - { - mainDirs.resize(0); - return; + mainDirs.resize(numClusters); + for (int i = 0; i < numClusters; i++) + mainDirs[i] = clustersAverage[i]; } - const int numIt = 200; - int numDirs = inDirs.size(), crit = 0, numClusters = 0; - if(numDirs < 4) + template + typename TODEstimatorImageFilter::DirType + TODEstimatorImageFilter::GetNewClusterAverage(int numCluster, DirVectorType &indirs, std::vector &cluster) { - numClusters = numDirs; - mainDirs.resize(numClusters); - for(int i = 0; i < numClusters; i++) + DirType sum, sumNorm; + DirVectorType dirs; + dirs = indirs; + for (int j = 0; j < 3; j++) + sum[j] = 0; + + int size = 0; + for (int i = 0; i < dirs.size(); i++) { - mainDirs[i] = inDirs[i]; + if (cluster[i] == numCluster) + { + for (int j = 0; j < 3; j++) + sum[j] = sum[j] + dirs[i][j]; + ++size; + } } - return; - } - else - numClusters = 4; - - DirVectorType clustersAverage(numClusters); - std::vector clusters(numDirs); - - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_int_distribution<> distrib(0, numDirs); + for (int j = 0; j < 3; j++) + sum[j] = sum[j] / size; + anima::Normalize(sum, sumNorm); + return sumNorm; + } - for(int i = 0; i < numClusters; i++) + template + double + TODEstimatorImageFilter::GetEuclideanDistance(DirType dir1, DirType dir2) { - int k = distrib(gen); - clustersAverage[i] = inDirs[k]; + return std::sqrt(std::pow((dir2[0] - dir1[0]), 2) + std::pow((dir2[1] - dir1[1]), 2) + std::pow((dir2[2] - dir1[2]), 2)); } - int It = 0; - bool stop = false; - double dist = 0; - while(It < numIt && !stop) + template + void + TODEstimatorImageFilter::GetSHCoef(DirType tmpDir, OutputImagePixelType &coefs) { - std::vector clustersOld(numDirs); - for(int i = 0; i < numDirs; i++) - clustersOld[i] = clusters[i]; - - for(int i = 0; i < numDirs; i++) + int pos = 0, T = this->GetOutput()->GetVectorLength(); + OutputImagePixelType tmp, resSH, rotatedModel(m_VectorLength); + tmp.SetSize(T); + DirType dir; + this->GetOutput()->TransformLocalVectorToPhysicalVector(tmpDir, dir); + for (int k = 0; k <= m_LOrder; k += 2) { - double minDist = 1e10; - for(int j = 0; j < numClusters; j++) + for (int m = -k; m <= k; m++) { - // double dist = getEuclideanDistance(inDirs[i], clustersAverage[j]); - // dist = std::sqrt(std::pow((clustersAverage[j][0] - inDirs[i][0]),2) + std::pow((clustersAverage[j][1] - inDirs[i][1]),2) + std::pow((clustersAverage[j][2] - inDirs[i][2]),2)); - dist = std::sqrt((clustersAverage[j][0] - inDirs[i][0]) * (clustersAverage[j][0] - inDirs[i][0]) + (clustersAverage[j][1] - inDirs[i][1]) * (clustersAverage[j][1] - inDirs[i][1]) + (clustersAverage[j][2] - inDirs[i][2]) * (clustersAverage[j][2] - inDirs[i][2])); - if(dist < minDist) - { - minDist = dist; - clusters[i] = j; - } + tmp[pos] = m_GaussCoefs[pos]; + pos++; } } - crit = 0; - for(int i = 0; i < numClusters; i++) - { - clustersAverage[i] = getNewClusterAverage(i, inDirs, clusters); - } - - for(int i = 0; i < numDirs; i++) - crit += clusters[i] - clustersOld[i]; - - if(crit == 0) - stop = true; - - ++It; - } - - mainDirs.resize(numClusters); - for(int i = 0; i < numClusters; i++) - mainDirs[i] = clustersAverage[i]; - + itk::Matrix itkRotationMatrix = anima::GetRotationMatrixFromVectors(m_CstDir, dir); + vnl_matrix rotationMatrix; + rotationMatrix.set_size(3, 3); + for (int x = 0; x < 3; ++x) + for (int y = 0; y < 3; ++y) + rotationMatrix.put(x, y, itkRotationMatrix[y][x]); -} + std::vector eulerAngles; + anima::GetEulerAnglesFromRotationMatrix(rotationMatrix, eulerAngles); + double a = rotationMatrix(2, 2); + rotatedModel = tmp; -typename TODEstimatorImageFilter::DirType TODEstimatorImageFilter ::getNewClusterAverage(int numCluster, DirVectorType &indirs, std::vector &cluster) -{ - DirType sum, sumNorm; - DirVectorType dirs; - dirs = indirs; - for(int j = 0; j < 3; j++) - sum[j] = 0; - - int size = 0; - for(int i = 0 ; i < dirs.size(); i++) - { - if(cluster[i] == numCluster) + vnl_matrix ODFRotationMatrix; + for (unsigned int l = 0; l <= m_LOrder; l += 2) { - for(int j = 0; j < 3; j++) - sum[j] = sum[j] + dirs[i][j]; - ++size; - } - } - for(int j = 0; j < 3; j++) - sum[j] = sum[j]/size; - - - anima::Normalize(sum, sumNorm); - return sumNorm; -} - - -double TODEstimatorImageFilter::getEuclideanDistance(DirType dir1, DirType dir2) -{ - return std::sqrt(std::pow((dir2[0] - dir1[0]),2) + std::pow((dir2[1] - dir1[1]),2) + std::pow((dir2[2] - dir1[2]),2)); -} + anima::EstimateLocalODFRotationMatrix(ODFRotationMatrix, l, eulerAngles[0], eulerAngles[1], eulerAngles[2]); + unsigned int mBaseInd = (l * l + l + 2) / 2 - l - 1; + for (unsigned int m = 0; m <= 2 * l; ++m) + { + rotatedModel[mBaseInd + m] = 0; -void TODEstimatorImageFilter::getSHCoef(DirType tmpDir, VectorType &coefs) -{ - int pos = 0, T = this->GetOutput()->GetVectorLength(); - VectorType tmp, resSH, rotatedModel; - tmp.SetSize(T); - DirType dir; - this->GetOutput()->TransformLocalVectorToPhysicalVector(tmpDir, dir); - // this->GetOutput()->TransformPhysicalVectorToLocalVector(tmpDir, dir); - // dir[0] *= -1; - // dir[1] *= -1; - // dir[2] *= -1; - - for(int k = 0; k <= m_LOrder; k+=2) - { - for(int m = -k; m <= k; m ++) - { - tmp[pos] = m_GaussCoefs[pos]; - pos++; + for (unsigned int mp = 0; mp <= 2 * l; ++mp) + rotatedModel[mBaseInd + m] += ODFRotationMatrix(m, mp) * tmp[mBaseInd + mp]; + } } + coefs = rotatedModel; } - // if(this->GetOutput()->GetDirection()[0][0] > 0) - // dir[0] *= -1; - // if(this->GetOutput()->GetDirection()[1][1] < 0) - // dir[1] *= -1; - // if(this->GetOutput()->GetDirection()[2][2] < 0) - // dir[2] *= -1; - // itk::Matrix itkRotationMatrix = anima::GetRotationMatrixFromVectors(m_CstDir, dir); - itk::Matrix itkRotationMatrix = anima::GetRotationMatrixFromVectors(dir, m_CstDir); - - vnl_matrix rotationMatrix; - rotationMatrix.set_size(3,3); - for(int x = 0; x < 3; ++x) - for(int y = 0; y < 3; ++y) - rotationMatrix.put(x, y, itkRotationMatrix[y][x]); - - // vnl_matrix rotationMatrix = this->GetRotationMatrix(CstDir, dir); - std::vector eulerAngles; - anima::GetEulerAnglesFromRotationMatrix(rotationMatrix, eulerAngles); - - rotatedModel = tmp; - - vnl_matrix ODFRotationMatrix; - for (unsigned int l = 0;l <= m_LOrder;l += 2) + template + void + TODEstimatorImageFilter::AverageODFs(std::vector &vecCoefs, OutputImagePixelType &resOdf) { - anima::EstimateLocalODFRotationMatrix(ODFRotationMatrix, l, eulerAngles[0],eulerAngles[1],eulerAngles[2]); + int numImages = vecCoefs.size(); + + std::vector> odfHistos(numImages); + std::vector vecSqrtCoefs(numImages); + std::vector avgSqrtHisto(m_NbSample); + OutputImagePixelType avgSqrtCoefs(m_VectorLength); - unsigned int mBaseInd = (l*l + l + 2)/2 - l - 1; - for (unsigned int m = 0;m <= 2*l;++m) + for (int i = 0; i < numImages; i++) { - rotatedModel[mBaseInd + m] = 0; + odfHistos[i].resize(m_NbSample); - for (unsigned int mp = 0;mp <= 2*l;++mp) - rotatedModel[mBaseInd + m] += ODFRotationMatrix(m,mp)*tmp[mBaseInd + mp]; + this->DiscretizeODF(vecCoefs[i], odfHistos[i]); + vecSqrtCoefs[i] = GetSquareRootODF(odfHistos[i]); } - } - coefs = rotatedModel; -} - -void TODEstimatorImageFilter::averageODFs(std::vector &vecCoefs, VectorType &resOdf) -{ - int numImages = vecCoefs.size(); - - std::vector> odfHistos(numImages); - std::vector vecSqrtCoefs(numImages); - std::vector avgSqrtHisto(m_NbSample); - VectorType avgSqrtCoefs(m_VectorLength); - - for(int i = 0; i < numImages; i++) - { - odfHistos[i].resize(m_NbSample); - - this->discretizeODF(vecCoefs[i], odfHistos[i]); - vecSqrtCoefs[i] = getSquareRootODF(odfHistos[i]); + this->GetAverageCoefs(vecSqrtCoefs, avgSqrtCoefs); + this->DiscretizeODF(avgSqrtCoefs, avgSqrtHisto); + resOdf = this->GetSquareODF(avgSqrtHisto); } - this->getAverageCoefs(vecSqrtCoefs, avgSqrtCoefs); - this->discretizeODF(avgSqrtCoefs, avgSqrtHisto); - resOdf = this->getSquareODF(avgSqrtHisto); -} - -void TODEstimatorImageFilter::precomputeSH() -{ - int k = 0; - - std::vector tmpDir(2), revDir(3); - for(int i = 0; i < m_NbSample; i++) + template + void + TODEstimatorImageFilter::PrecomputeSH() { - anima::TransformCartesianToSphericalCoordinates(m_SphereSampl[i], tmpDir); - int c = 0; - for(double l = 0; l <= m_LOrder; l+=2) + int k = 0; + + std::vector tmpDir(2), revDir(3); + for (int i = 0; i < m_NbSample; i++) { - for(double m = -l; m <= l; m++) + anima::TransformCartesianToSphericalCoordinates(m_SphereSampl[i], tmpDir); + int c = 0; + for (double l = 0; l <= m_LOrder; l += 2) { - m_SpherHarm(k, c++) = m_ODFSHBasis->getNthSHValueAtPosition(l, m, tmpDir[0], tmpDir[1]); + for (double m = -l; m <= l; m++) + { + m_SpherHarm(k, c++) = m_ODFSHBasis->getNthSHValueAtPosition(l, m, tmpDir[0], tmpDir[1]); + } } + k++; } - k++; } -} - -void TODEstimatorImageFilter::discretizeODF(VectorType ODFCoefs, std::vector &ODFDiscret) -{ - int k = 0; - double resVal2 = 0; - - std::vector tmpDir(2), revDir(3); - for(int i = 0; i < m_NbSample; i++) + template + void + TODEstimatorImageFilter::DiscretizeODF(OutputImagePixelType ODFCoefs, std::vector &ODFDiscret) { - anima::TransformCartesianToSphericalCoordinates(m_SphereSampl[i], tmpDir); - resVal2 = m_ODFSHBasis->getValueAtPosition(ODFCoefs, tmpDir[0], tmpDir[1]); - - ODFDiscret[i] = resVal2; - k++; - } -} - - - -typename TODEstimatorImageFilter::VectorType TODEstimatorImageFilter::getSquareRootODF(std::vector ODFDiscret) -{ - vnl_matrix Odf(m_NbSample, 1); - vnl_matrix Coef(m_VectorLength, 1); - - for(int i = 0; i < m_NbSample; ++i) - Odf(i, 0) = std::sqrt(ODFDiscret[i]); + int k = 0; + double resVal2 = 0; - Coef = vnl_matrix_inverse(m_SpherHarm.transpose() * m_SpherHarm).as_matrix() * m_SpherHarm.transpose()*Odf; + std::vector tmpDir(2), revDir(3); + for (int i = 0; i < m_NbSample; i++) + { + anima::TransformCartesianToSphericalCoordinates(m_SphereSampl[i], tmpDir); + resVal2 = m_ODFSHBasis->getValueAtPosition(ODFCoefs, tmpDir[0], tmpDir[1]); - VectorType ModelValue(m_VectorLength); - for(int i = 0; i < m_VectorLength; ++i) - { - ModelValue[i] = Coef(i, 0); + ODFDiscret[i] = resVal2; + k++; + } } - return ModelValue; -} - - -typename TODEstimatorImageFilter::VectorType -TODEstimatorImageFilter::getSquareODF(std::vector ODFDiscret) -{ - vnl_matrix Odf(m_NbSample, 1); - vnl_matrix Coef(m_VectorLength, 1); - - for(int i = 0; i < m_NbSample; ++i) - Odf(i, 0) = ODFDiscret[i] * ODFDiscret[i]; - Coef = vnl_matrix_inverse(m_SpherHarm.transpose() * m_SpherHarm).as_matrix() * m_SpherHarm.transpose()*Odf; - - VectorType ModelValue(m_VectorLength); - for(int i = 0; i < m_VectorLength; ++i) + template + typename TODEstimatorImageFilter::OutputImagePixelType + TODEstimatorImageFilter::GetSquareRootODF(std::vector ODFDiscret) { - ModelValue[i] = Coef(i, 0); - } + vnl_matrix Odf(m_NbSample, 1); + vnl_matrix Coef(m_VectorLength, 1); - ModelValue[0] = 1/(2*sqrt(M_PI)); - return ModelValue; -} - - -void -TODEstimatorImageFilter::getAverageCoefs(std::vector &vecCoefs, VectorType &avgCoef) -{ - int numImage = vecCoefs.size(); + for (int i = 0; i < m_NbSample; ++i) + Odf(i, 0) = std::sqrt(ODFDiscret[i]); - std::vector> arrayCoef(numImage), arrayLogMap(numImage); - std::vector expMap(m_VectorLength), mean(m_VectorLength), nextMean(m_VectorLength), tangent(m_VectorLength); + Coef = vnl_matrix_inverse(m_SpherHarm.transpose() * m_SpherHarm).as_matrix() * m_SpherHarm.transpose() * Odf; - const int T = 100; - - for(int i = 0; i < numImage; i++) - { - arrayCoef[i].resize(m_VectorLength); - arrayLogMap[i].resize(m_VectorLength); - - for(int n = 0; n < m_VectorLength; n++) + OutputImagePixelType ModelValue(m_VectorLength); + for (int i = 0; i < m_VectorLength; ++i) { - arrayCoef[i][n] = vecCoefs[i][n]; - // nextMean[n] += 1/(double)numImage * arrayCoef[i][n]; + ModelValue[i] = Coef(i, 0); } + return ModelValue; } - nextMean = arrayCoef[0]; - - for(int t = 0; t < T; t++) + template + typename TODEstimatorImageFilter::OutputImagePixelType + TODEstimatorImageFilter::GetSquareODF(std::vector ODFDiscret) { - mean = nextMean; + vnl_matrix Odf(m_NbSample, 1); + vnl_matrix Coef(m_VectorLength, 1); - for(int i = 0; i < numImage; i++) - anima::sphere_log_map(arrayCoef[i], mean, arrayLogMap[i]); - - std::fill(tangent.begin(), tangent.end(), 0); - for(int n = 0; n < m_VectorLength; n++) - for(int i = 0; i < numImage; i++) - tangent[n] += 1/(double)numImage * arrayLogMap[i][n]; - - anima::sphere_exp_map(tangent, mean, nextMean); - } + for (int i = 0; i < m_NbSample; ++i) + Odf(i, 0) = ODFDiscret[i] * ODFDiscret[i]; - // nextMean[0] = 2.82094792e-01; - for(int i = 0; i < m_VectorLength; ++i) - avgCoef[i] = nextMean[i]; -} + Coef = vnl_matrix_inverse(m_SpherHarm.transpose() * m_SpherHarm).as_matrix() * m_SpherHarm.transpose() * Odf; + OutputImagePixelType ModelValue(m_VectorLength); + for (int i = 0; i < m_VectorLength; ++i) + { + ModelValue[i] = Coef(i, 0); + } + ModelValue[0] = 1 / (2 * sqrt(M_PI)); + return ModelValue; + } -void TODEstimatorImageFilter::processFiber(FiberType &fiber, baseSH &basis) -{ - const int T = (m_LOrder+1)*(m_LOrder+2)/2; - VectorType tmp, resSH, rotatedModel; - tmp.SetSize(T); - resSH.SetSize(T); - rotatedModel.SetSize(T); - // itk::Vector tmp, resSH; - DirType CstDir; - CstDir[0] = 0; - CstDir[1] = 0; - CstDir[2] = 1; - - - TOutputImage *output = this->GetOutput(); - int nbPoints = fiber.size(); - - for(int i = 0; i < nbPoints; i++) + template + void + TODEstimatorImageFilter::GetAverageCoefs(std::vector &vecCoefs, OutputImagePixelType &avgCoef) { - DirType dir = getFiberDirection(i, fiber); - // OutputImagePointer::DirectionType refDir = output->GetDirection(); - - if(output->GetDirection()[0][0] > 0) - dir[0] *= -1; - if(output->GetDirection()[1][1] < 0) - dir[1] *= -1; - if(output->GetDirection()[2][2] < 0) - dir[2] *= -1; + int numImage = vecCoefs.size(); + std::vector> arrayCoef(numImage), arrayLogMap(numImage); + std::vector expMap(m_VectorLength), mean(m_VectorLength), nextMean(m_VectorLength), tangent(m_VectorLength); + const int T = 100; - int pos = 0; - for(int k = 0; k <= m_LOrder; k+=2) + for (int i = 0; i < numImage; i++) { - for(int m = -k; m <= k; m ++) + arrayCoef[i].resize(m_VectorLength); + arrayLogMap[i].resize(m_VectorLength); + + for (int n = 0; n < m_VectorLength; n++) { - tmp[pos] = m_GaussCoefs[pos]; - pos++; + arrayCoef[i][n] = vecCoefs[i][n]; } } - itk::Matrix itkRotationMatrix = anima::GetRotationMatrixFromVectors(CstDir, dir); - vnl_matrix rotationMatrix; - rotationMatrix.set_size(3,3); - for(int x = 0; x < 3; ++x) - for(int y = 0; y < 3; ++y) - rotationMatrix.put(x, y, itkRotationMatrix[y][x]); - - // vnl_matrix rotationMatrix = this->GetRotationMatrix(CstDir, dir); - std::vector eulerAngles; - anima::GetEulerAnglesFromRotationMatrix(rotationMatrix, eulerAngles); - - rotatedModel = tmp; + nextMean = arrayCoef[0]; - vnl_matrix ODFRotationMatrix; - for (unsigned int l = 0;l <= m_LOrder;l += 2) + for (int t = 0; t < T; t++) { - anima::EstimateLocalODFRotationMatrix(ODFRotationMatrix, l, eulerAngles[0],eulerAngles[1],eulerAngles[2]); + mean = nextMean; - unsigned int mBaseInd = (l*l + l + 2)/2 - l - 1; - for (unsigned int m = 0;m <= 2*l;++m) - { - rotatedModel[mBaseInd + m] = 0; + for (int i = 0; i < numImage; i++) + anima::sphere_log_map(arrayCoef[i], mean, arrayLogMap[i]); - for (unsigned int mp = 0;mp <= 2*l;++mp) - rotatedModel[mBaseInd + m] += ODFRotationMatrix(m,mp)*tmp[mBaseInd + mp]; - } - } - PointType point = getCenterVoxel(i, fiber); - itk::Index<3> index; + std::fill(tangent.begin(), tangent.end(), 0); + for (int n = 0; n < m_VectorLength; n++) + for (int i = 0; i < numImage; i++) + tangent[n] += 1 / (double)numImage * arrayLogMap[i][n]; - output->TransformPhysicalPointToIndex(point, index); - // index[0] *= -1; - // index[1] *= -1; - // resSH = (output->GetPixel(index) + tmp); - // resSH /= resSH[0]; - output->SetPixel(index, rotatedModel); - - resSH.Fill(0); - tmp.Fill(0); - rotatedModel.Fill(0); + anima::sphere_exp_map(tangent, mean, nextMean); + } + // nextMean[0] = 2.82094792e-01; + for (int i = 0; i < m_VectorLength; ++i) + avgCoef[i] = nextMean[i]; } -} - -typename TODEstimatorImageFilter::DirType -TODEstimatorImageFilter::getFiberDirection(int index, FiberType &fiber) -{ - DirType resDir; - if(index == 0) - { - for(int i = 0; i < 3; i++) - resDir[i] = fiber[index+1][i] - fiber[index][i]; - } - else if(index == fiber.size() - 1) + template + typename TODEstimatorImageFilter::DirType + TODEstimatorImageFilter::GetFiberDirection(int index, FiberType &fiber) { - for(int i = 0; i < 3; i++) - resDir[i] = fiber[index][i] - fiber[index-1][i]; - } - else - { - for(int i = 0; i < 3; i++) - resDir[i] = fiber[index-1][i] - fiber[index+1][i]; - } - - return resDir; -} - + DirType resDir; + if (index == fiber.size() - 1) + { + for (int i = 0; i < 3; i++) + resDir[i] = fiber[index][i] - fiber[index - 1][i]; + } + else + { + for (int i = 0; i < 3; i++) + resDir[i] = fiber[index + 1][i] - fiber[index][i]; + } -typename TODEstimatorImageFilter::PointType -TODEstimatorImageFilter::getCenterVoxel(int index, FiberType &fiber) -{ - PointType resPoint; - for(int i = 0; i < 3; i++) - { - resPoint[i] = fiber[index][i]; + anima::Normalize(resDir, resDir); + return resDir; } - return resPoint; -} - - -typename TODEstimatorImageFilter::FiberType TODEstimatorImageFilter::readFiber(vtkIdType numberOfPoints, const vtkIdType *indices, vtkPoints *points) -{ - FiberType fiber; - // fiber.resize(numberOfPoints); - for (vtkIdType i = 0; i < numberOfPoints; i++) + template + typename TODEstimatorImageFilter::PointType + TODEstimatorImageFilter::GetCenterVoxel(int index, FiberType &fiber) { - double tmpPoint[3]; - points->GetPoint(indices[i], tmpPoint); - PointType point; - for(int i = 0; i < 3; i++) - point[i] = tmpPoint[i]; + PointType resPoint; + for (int i = 0; i < 3; i++) + { + resPoint[i] = fiber[index][i]; + } - fiber.push_back(tmpPoint); + return resPoint; } - return fiber; -} - - -vnl_matrix TODEstimatorImageFilter::GetRotationMatrix(DirType dir1, DirType dir2) -{ - DirType nu; - anima::ComputeCrossProduct(dir2, dir1, nu); - - double s = anima::ComputeNorm(nu); - double c = anima::ComputeScalarProduct(dir2, dir1); - - vnl_matrix R; - vnl_matrix MatNu; - vnl_matrix SqrMatNu; - MatNu.set_size(3,3); - MatNu.fill(0.0); - - MatNu.put(0,1,-nu[2]); - MatNu.put(0,2, nu[1]); - MatNu.put(1,0, nu[2]); - MatNu.put(1,2,-nu[0]); - MatNu.put(2,0,-nu[1]); - MatNu.put(2,1, nu[0]); - - - // MatNu[1][0] = -nu[2]; - // MatNu[2][0] = nu[1]; - // MatNu[0][1] = nu[2]; - // MatNu[2][1] = -nu[0]; - // MatNu[0][2] = -nu[1]; - // MatNu[1][2] = nu[0]; - - SqrMatNu = MatNu*MatNu; - - R.set_size(3,3); - R.set_identity(); - R = R + MatNu + SqrMatNu*(1-c)/(s*s); - // R = R + MatNu + SqrMatNu*(1/1+c); + template + typename TODEstimatorImageFilter::FiberType + TODEstimatorImageFilter::ReadFiber(vtkIdType numberOfPoints, const vtkIdType *indices, vtkPoints *points) + { + FiberType fiber; + for (vtkIdType i = 0; i < numberOfPoints; i++) + { + double tmpPoint[3]; + points->GetPoint(indices[i], tmpPoint); + PointType point; + for (int i = 0; i < 3; i++) + point[i] = tmpPoint[i]; - return R; + fiber.push_back(tmpPoint); + } -} + return fiber; + } } // end namespace anima -