diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp index 53a38e5e80..8e1100202f 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp @@ -1,340 +1,392 @@ /*=================================================================== 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()); + } + + void set_reference_measuremnt(const double & x) + { + log_S_0 = x; + } + + double log_S_0; + vnl_vector measurements; + vnl_vector b; + + int N; + + lestSquaresFunction(unsigned int number_of_measurements) : + vnl_least_squares_function(2, number_of_measurements, no_gradient) + { + N = get_number_of_residuals(); + } + + void f(const vnl_vector& x, vnl_vector& fx) { + + double D = x[0]; + double K = x[1]; + + 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); for(BValueMap::const_iterator it = m_B_ValueMap.begin();it != m_B_ValueMap.end(); it++) { if((*it).first == 0.0) continue; // if any #ShellDirection < 6 --> itkException (No interpolation possible) if((*it).second.size() < 6){ MITK_INFO << "Abort: No interpolation possible Shell-" << (*it).first << " has less than 6 directions."; itkExceptionMacro(<<"No interpolation possible"); } } m_ShellInterpolationMatrixVector.reserve(m_B_ValueMap.size()-1); - vnl_matrix lsfParameterMatrix(m_B_ValueMap.size()-1,2); + + m_bValueVector = vnl_vector(2); // 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); - - // lsf Matrix A=|b_1 1/6| - // |. | - // |b_n 1/6| - - lsfParameterMatrix.put(shellIndex, 0, it->first); - lsfParameterMatrix.put(shellIndex, 1, 1./6.); - - + m_bValueVector.put(shellIndex,it->first); ++shellIndex; } - MITK_INFO << "lsfParameteR(1)"; - MITK_INFO << lsfParameterMatrix; - - vnl_matrix_inverse A_A(lsfParameterMatrix.transpose() * lsfParameterMatrix); - MITK_INFO << "lsfParameteR(2)"; - MITK_INFO << A_A; - - m_lsfParameterMatrix = A_A.inverse() * lsfParameterMatrix.transpose(); - MITK_INFO << "lsfParameteR(3)"; - MITK_INFO << m_lsfParameterMatrix; 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 = OutputImageType::New(); 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); // ** walking over each Voxel while(!iit.IsAtEnd()) { InputPixelType b = iit.Get(); BZeroAverage=0.0; for(unsigned int i = 0 ; i < BZeroIndices.size(); i++){ //MITK_INFO << "BValues("<second; vnl_vector InterpVector(currentShell.size()); // - get raw Signal for currente shell for(unsigned int i = 0 ; i < currentShell.size(); i++) InterpVector.put(i,b[currentShell[i]]); - //- normalization of the raw Signal divided by Reference Signal - S_S0Normalization(InterpVector, BZeroAverage); + logarithm(InterpVector); //- interpolate the Signal if necessary using corresponding interpolationSHBasis SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; SignalMatrix.set_column(shellIndex, SignalVector); - shellIterator++; shellIndex++; } // row_i = {D, D^2*K} vnl_matrix lsfCoeffs(m_allDirectionsSize , 2); - calculateLsfCoeffs(lsfCoeffs,SignalMatrix); - + calculateCoeffs(lsfCoeffs,SignalMatrix, m_bValueVector, std::log(BZeroAverage)); calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs,m_TargetB_Value,BZeroAverage); for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); oit.Set(out); - MITK_INFO << out; + // MITK_INFO << out; ++oit; ++iit; } } template void MultiShellRadialAdcKurtosisImageFilter -::S_S0Normalization( vnl_vector & vec, const double & S0 ) +::logarithm( vnl_vector & vec) { for(unsigned int i = 0; i < vec.size(); i++){ - vec[i] = std::log( vec[i]/S0 ); - if(vec[i]>1.0) vec[i] = 1.; - if(vec[i]<0.0) vec[i] = 0.; - + if(vec[i] < 0.0) vec[i] = 0.00001; + vec[i] = std::log( vec[i] ); } + // MITK_INFO << vec; } template void MultiShellRadialAdcKurtosisImageFilter -::calculateLsfCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix) +::calculateCoeffs(vnl_matrix &lsfCoeffs, const vnl_matrix & SignalMatrix, const vnl_vector & bValueVector, const double & reference_b_value) { + lestSquaresFunction model(bValueVector.size()); + vnl_vector initalGuess(2); + + for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++) { - // x = (A' A)^-1 A' b - vnl_vector lsfCoeffsVector(m_lsfParameterMatrix * SignalMatrix.get_row(i)); - lsfCoeffs.set_row(i, lsfCoeffsVector); + model.set_measurements(SignalMatrix.get_row(i)); + model.set_bvalues(bValueVector); + model.set_reference_measuremnt(reference_b_value); + + initalGuess.put(0, 0.8); + initalGuess.put(1, 1.0); + + // 3. perform least squares minimization of model results by adapting model parameter + vnl_levenberg_marquardt minimizer(model); + minimizer.set_f_tolerance( 1e-15 ); + minimizer.set_max_function_evals( 10000 ); + bool status = minimizer.minimize(initalGuess); + if(!status) + { + MITK_INFO<< "Minimizer f Error: " << minimizer.get_f_tolerance(); + MITK_INFO<< "Minimizer end Error: " << minimizer.get_end_error(); + } + + lsfCoeffs.set_row(i, initalGuess); } + + MITK_INFO << lsfCoeffs; } 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 &bValue, const double & referenceSignal) { for(unsigned int i = 0 ; i < lsfCoeffs.rows();i++){ // S = S0 * e^(-b*D + 1/6*D^2*K) double D = lsfCoeffs(i,0); - double K = lsfCoeffs(i,1) / (D*D); - MITK_INFO << D << " " << K; - vec[i] = referenceSignal * exp((-bValue) * D + 1./6. * D* D * K); + double K = lsfCoeffs(i,1); + //MITK_INFO << D << " " << K; + vec[i] = referenceSignal * exp((-bValue) * D + 1./6. * bValue * bValue * D*D*K); } } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h index 9e6d22c997..50192e1e19 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h @@ -1,120 +1,120 @@ /*=================================================================== 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 S_S0Normalization( vnl_vector & vec, const double & S0 ); - void calculateLsfCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix); + 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); 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_matrix m_lsfParameterMatrix; + 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