diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp index e90245596e..6343288456 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp @@ -1,399 +1,397 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_MultiShellDirectionalKurtosisFitImageFilter_cpp_ #define _itk_MultiShellDirectionalKurtosisFitImageFilter_cpp_ #endif #define _USE_MATH_DEFINES #include "itkMultiShellRadialAdcKurtosisImageFilter.h" #include #include #include "mitkDiffusionFunctionCollection.h" #include "vnl/vnl_least_squares_function.h" #include "vnl/algo/vnl_levenberg_marquardt.h" namespace itk { /** baseclass for IVIM fitting algorithms */ struct lestSquaresFunction: public vnl_least_squares_function { void set_measurements(const vnl_vector& x) { measurements.set_size(x.size()); measurements.copy_in(x.data_block()); } void set_bvalues(const vnl_vector& x) { - b.set_size(x.size()); - b.copy_in(x.data_block()); + bValueVector.set_size(x.size()); + bValueVector.copy_in(x.data_block()); } - void set_reference_measuremnt(const double & x) - { - reference_measurement = x; - } - - double reference_measurement; vnl_vector measurements; - vnl_vector b; - + vnl_vector bValueVector; int N; lestSquaresFunction(unsigned int number_of_measurements) : - vnl_least_squares_function(2, number_of_measurements, no_gradient) + vnl_least_squares_function(3 /*number of unknowns*/, number_of_measurements, no_gradient) { N = get_number_of_residuals(); } void f(const vnl_vector& x, vnl_vector& fx) { const double & D = x[0]; const double & K = x[1]; + const double & S0= x[2]; + const vnl_vector & b = bValueVector; for(int s=0; s MultiShellRadialAdcKurtosisImageFilter ::MultiShellRadialAdcKurtosisImageFilter() { this->SetNumberOfRequiredInputs( 1 ); //this->SetNumberOfThreads(1); } template void MultiShellRadialAdcKurtosisImageFilter ::BeforeThreadedGenerateData() { // test whether BvalueMap contains all necessary information if(m_B_ValueMap.size() == 0) { itkWarningMacro(<< "No BValueMap given: create one using GradientDirectionContainer"); GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = m_OriginalGradientDirections->Begin(); gdcit != m_OriginalGradientDirections->End(); ++gdcit) { double bValueKey = int(((m_B_Value * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; m_B_ValueMap[bValueKey].push_back(gdcit.Index()); } } //# BValueMap contains no bZero --> itkException if(m_B_ValueMap.find(0.0) == m_B_ValueMap.end()) { MITK_INFO << "No ReferenceSignal (BZeroImages) found!"; itkExceptionMacro(<< "No ReferenceSignal (BZeroImages) found!"); } // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions m_allDirectionsIndicies = mitk::gradients::GetAllUniqueDirections(m_B_ValueMap, m_OriginalGradientDirections); // [sizeAllDirections] size of GradientContainer cointaining all unique directions m_allDirectionsSize = m_allDirectionsIndicies.size(); m_TargetGradientDirections = mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(m_B_ValueMap,m_OriginalGradientDirections); m_ShellInterpolationMatrixVector.reserve(m_B_ValueMap.size()-1); m_bValueVector = vnl_vector(m_B_ValueMap.size()-1); // for each shell BValueMap::const_iterator it = m_B_ValueMap.begin(); it++; //skip bZeroIndices unsigned int shellIndex = 0; for(;it != m_B_ValueMap.end();++it) { //- calculate maxShOrder const IndicesVector currentShell = it->second; unsigned int SHMaxOrder = 12; while( ((SHMaxOrder+1)*(SHMaxOrder+2)/2) > currentShell.size() && ((SHMaxOrder+1)*(SHMaxOrder+2)/2) >= 4 ) { SHMaxOrder -= 2 ; } //- get TragetSHBasis using allDirectionsContainer vnl_matrix sphericalCoordinates; sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(m_allDirectionsIndicies, m_OriginalGradientDirections); vnl_matrix TargetSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); //MITK_INFO << "TargetSHBasis " << TargetSHBasis.rows() << " x " << TargetSHBasis.cols(); //- get ShellSHBasis using currentShellDirections sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(currentShell, m_OriginalGradientDirections); vnl_matrix ShellSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); //MITK_INFO << "ShellSHBasis " << ShellSHBasis.rows() << " x " << ShellSHBasis.cols(); //- calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] vnl_matrix_inverse invShellSHBasis(ShellSHBasis); vnl_matrix shellInterpolationMatrix = TargetSHBasis * invShellSHBasis.pinverse(); //MITK_INFO << "shellInterpolationMatrix " << shellInterpolationMatrix.rows() << " x " << shellInterpolationMatrix.cols(); //- save interpolationSHBasis m_ShellInterpolationMatrixVector.push_back(shellInterpolationMatrix); m_bValueVector.put(shellIndex,it->first); ++shellIndex; } m_WeightsVector.reserve(m_B_ValueMap.size()-1); BValueMap::const_iterator itt = m_B_ValueMap.begin(); itt++; // skip ReferenceImages //- calculate Weights [Weigthing = shell_size / max_shell_size] unsigned int maxShellSize = 0; for(;itt != m_B_ValueMap.end(); itt++){ if(itt->second.size() > maxShellSize) maxShellSize = itt->second.size(); } itt = m_B_ValueMap.begin(); itt++; // skip ReferenceImages for(;itt != m_B_ValueMap.end(); itt++){ m_WeightsVector.push_back(itt->second.size() / (double)maxShellSize); MITK_INFO << m_WeightsVector.back(); } // initialize output image typename OutputImageType::Pointer outImage = static_cast(ProcessObject::GetOutput(0)); outImage->SetSpacing( this->GetInput()->GetSpacing() ); outImage->SetOrigin( this->GetInput()->GetOrigin() ); outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction using bZeroDirection+AllDirectionsContainer outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetVectorLength( 1+m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); BValueMap::iterator ittt = m_B_ValueMap.begin(); ittt++; // skip bZeroImages corresponding to 0-bValue m_TargetB_Value = 0; while(ittt!=m_B_ValueMap.end()) { m_TargetB_Value += ittt->first; ittt++; } m_TargetB_Value /= (double)(m_B_ValueMap.size()-1); MITK_INFO << "Input:" << std::endl << std::endl << " GradientDirections: " << m_OriginalGradientDirections->Size() << std::endl << " Shells: " << (m_B_ValueMap.size() - 1) << std::endl << " ReferenceImages: " << m_B_ValueMap[0.0].size() << std::endl; MITK_INFO << "Output:" << std::endl << std::endl << " OutImageVectorLength: " << outImage->GetVectorLength() << std::endl << " TargetDirections: " << m_allDirectionsSize << std::endl << " TargetBValue: " << m_TargetB_Value << std::endl << std::endl; } template void MultiShellRadialAdcKurtosisImageFilter ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/) { // Get input gradient image pointer typename InputImageType::Pointer inputImage = static_cast< InputImageType * >(ProcessObject::GetInput(0)); // ImageRegionIterator for the input image ImageRegionIterator< InputImageType > iit(inputImage, outputRegionForThread); iit.GoToBegin(); // Get output gradient image pointer typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(ProcessObject::GetOutput(0)); // ImageRegionIterator for the output image ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); // calculate target bZero-Value [b0_t] const IndicesVector BZeroIndices = m_B_ValueMap[0.0]; double BZeroAverage = 0.0; unsigned int numberOfShells = m_B_ValueMap.size()-1; // create empty nxm SignalMatrix containing n->signals/directions (in case of interpolation ~ sizeAllDirections otherwise the size of any shell) for m->shells vnl_matrix SignalMatrix(m_allDirectionsSize, numberOfShells); // create nx1 targetSignalVector vnl_vector SignalVector(m_allDirectionsSize); OutputPixelType out; InputPixelType b; BValueMap::const_iterator shellIterator; vnl_matrix lsfCoeffs; vnl_vector InterpVector; // ** walking over each Voxel while(!iit.IsAtEnd()) { b = iit.Get(); BZeroAverage=0.0; for(unsigned int i = 0 ; i < BZeroIndices.size(); i++){ BZeroAverage += b[BZeroIndices[i]]; } BZeroAverage /= (double)BZeroIndices.size(); out = oit.Get(); out.Fill(0.0); out.SetElement(0,BZeroAverage); shellIterator = m_B_ValueMap.begin(); shellIterator++; //skip bZeroImages unsigned int shellIndex = 0; while(shellIterator != m_B_ValueMap.end()) { // reset Data SignalVector.fill(0.0); // - get the RawSignal const IndicesVector & currentShell = shellIterator->second; InterpVector.set_size(currentShell.size()); // - get raw Signal for currente shell for(unsigned int i = 0 ; i < currentShell.size(); i++) InterpVector.put(i,b[currentShell[i]]); //- interpolate the Signal using corresponding interpolationSHBasis SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; SignalMatrix.set_column(shellIndex, SignalVector); shellIterator++; shellIndex++; } - lsfCoeffs.set_size(m_allDirectionsSize , 2); - calculateCoeffs(lsfCoeffs,SignalMatrix, m_bValueVector, BZeroAverage); - calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs,m_TargetB_Value,BZeroAverage); + lsfCoeffs.set_size(m_allDirectionsSize , 3 /*Number of unknowns [ADC AKC S0]*/); + calculateCoeffs(lsfCoeffs, SignalMatrix, BZeroAverage); + calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs,m_TargetB_Value); for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); oit.Set(out); ++oit; ++iit; } } template void MultiShellRadialAdcKurtosisImageFilter ::logarithm( vnl_vector & vec) { for(unsigned int i = 0; i < vec.size(); i++){ if(vec[i] <= 0.0) vec[i] = 0.00001; vec[i] = std::log( vec[i] ); } } template void MultiShellRadialAdcKurtosisImageFilter -::calculateCoeffs(vnl_matrix &lsfCoeffs, const vnl_matrix & SignalMatrix, const vnl_vector & bValueVector, const double & reference_b_value) +::calculateCoeffs(vnl_matrix &lsfCoeffs, const vnl_matrix & SignalMatrix, const double & S0) { - vnl_vector initalGuess(2); + vnl_vector initalGuess(3); // initialize Least Squres Function - lestSquaresFunction model(bValueVector.size()); + // SignalMatrix.cols() defines the number of shells/measurement points + lestSquaresFunction model(SignalMatrix.cols()); + // initialize Levenberg Marquardt vnl_levenberg_marquardt minimizer(model); minimizer.set_max_function_evals(50000); // Iterations minimizer.set_f_tolerance(1e-10); // Function tolerance - minimizer.set_epsilon_function(1e-10); // Epsilon // for each Direction calculate LSF Coeffs ADC & AKC for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++) { // set Signal Vector model.set_measurements(SignalMatrix.get_row(i)); - // set BValue Vector e.g.: [1000, 2000, 3000] - model.set_bvalues(bValueVector); - // set BZero Value - model.set_reference_measuremnt(reference_b_value); + // set BValue Vector e.g.: [1000, 2000, 3000] <- shell b Values + model.set_bvalues(m_bValueVector); - // Start Vector [ADC, AKC] + // Start Vector [ADC, AKC, S0] initalGuess.put(0, 0.f); initalGuess.put(1, 0.f); + initalGuess.put(2, S0); // user B=0 Signal of processing Voxel // start Levenberg-Marquardt bool status = minimizer.minimize_without_gradient(initalGuess); if(!status) { MITK_INFO<< "______________________________"; MITK_INFO<< "Minimizer f Error: " << minimizer.get_f_tolerance(); MITK_INFO<< "Minimizer end Error: " << minimizer.get_end_error(); MITK_INFO<< "Minimizer Evaluations: " << minimizer.get_num_evaluations(); MITK_INFO<< "Minimizer Iterations: " << minimizer.get_num_iterations(); MITK_INFO<< "______________________________"; }else{ std::cout << std::scientific << std::setprecision(5) - << reference_b_value << ";" // BZero - << initalGuess[0] << ";" // ADC - << initalGuess[1] << ";"; // AKC + << initalGuess[0] << ";" // fitted ADC + << initalGuess[1] << ";" // fitted AKC + << initalGuess[2] << ";" // fitted S0 value + << minimizer.get_end_error() << ";"; // End error for(unsigned int j = 0; j < SignalMatrix.get_row(i).size(); j++ ) std::cout << std::scientific << std::setprecision(5) - << SignalMatrix.get_row(i)[j] << ";"; + << SignalMatrix.get_row(i)[j] << ";"; // S_n Values corresponding to shell 1 to shell n std::cout << std::endl; } // Set Coeffs for each gradient direction lsfCoeffs.set_row(i, initalGuess); } } + template void MultiShellRadialAdcKurtosisImageFilter -::calculateSignalFromLsfCoeffs(vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double &bValue, const double & referenceSignal) +::calculateSignalFromLsfCoeffs(vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double &b /*target bValue*/) { // For each Direction for(unsigned int i = 0 ; i < lsfCoeffs.rows();i++){ const double & D = lsfCoeffs(i,0); const double & K = lsfCoeffs(i,1); + const double & S0 = lsfCoeffs(i,2); - vec[i] = referenceSignal * exp( -bValue * D + 1./6.* bValue * bValue * D * D * K); + vec[i] = S0 * exp( -b * D + 1./6.* b * b * D * D * K); // std::cout << "OUTVAL = " << i << std::endl; } } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h index 9d02cce3df..0d9d504862 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h @@ -1,120 +1,118 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _itk_MultiShellRadialAdcKurtosisImageFilter_h_ #define _itk_MultiShellRadialAdcKurtosisImageFilter_h_ #include #include #include namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class MultiShellRadialAdcKurtosisImageFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef MultiShellRadialAdcKurtosisImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MultiShellAdcAverageReconstructionImageFilter, ImageToImageFilter) typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef std::vector IndicesVector; typedef std::map BValueMap; GradientDirectionContainerType::Pointer GetOriginalGradientDirections(){return m_OriginalGradientDirections;} void SetOriginalGradientDirections(GradientDirectionContainerType::Pointer ptr){m_OriginalGradientDirections = ptr;} GradientDirectionContainerType::Pointer GetTargetGradientDirections(){return m_TargetGradientDirections;} double GetTargetB_Value(){return m_TargetB_Value;} double GetB_Value(){return m_B_Value;} void SetB_Value(double val){m_B_Value = val;} void SetOriginalBValueMap(BValueMap inp){m_B_ValueMap = inp;} protected: MultiShellRadialAdcKurtosisImageFilter(); ~MultiShellRadialAdcKurtosisImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); void logarithm( vnl_vector & vec); - void calculateCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix, const vnl_vector & bValueVector, const double & reference_b_value); - void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & bValue, const double & referenceSignal); - - + void calculateCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix, const double & S0); + void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & b); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions BValueMap m_B_ValueMap; double m_B_Value; double m_TargetB_Value; std::vector m_WeightsVector; std::vector > m_ShellInterpolationMatrixVector; vnl_vector m_bValueVector; std::vector m_bZeroIndicesSplitVectors; IndicesVector m_allDirectionsIndicies; unsigned int m_allDirectionsSize; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellRadialAdcKurtosisImageFilter.cpp" #endif #endif