diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp index 6343288456..5d67171037 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp @@ -1,397 +1,405 @@ /*=================================================================== 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) { bValueVector.set_size(x.size()); bValueVector.copy_in(x.data_block()); } + void set_reference_measurement(const double & x) + { + S0 = x; + } + vnl_vector measurements; vnl_vector bValueVector; + double S0; int N; lestSquaresFunction(unsigned int number_of_measurements) : - vnl_least_squares_function(3 /*number of unknowns*/, number_of_measurements, no_gradient) + vnl_least_squares_function(2 /*number of unknowns [ADC AKC]*/, 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 , 3 /*Number of unknowns [ADC AKC S0]*/); + lsfCoeffs.set_size(m_allDirectionsSize , 2 /*Number of unknowns [ADC AKC]*/); calculateCoeffs(lsfCoeffs, SignalMatrix, BZeroAverage); - calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs,m_TargetB_Value); + calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs, BZeroAverage, 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 double & S0) { - vnl_vector initalGuess(3); + vnl_vector initalGuess(2); // initialize Least Squres Function // 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 // 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] <- shell b Values model.set_bvalues(m_bValueVector); + model.set_reference_measurement(S0); // 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 + //initalGuess.put(2, S0); // user B=0 Signal of processing Voxel // start Levenberg-Marquardt - bool status = minimizer.minimize_without_gradient(initalGuess); + 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) << initalGuess[0] << ";" // fitted ADC << initalGuess[1] << ";" // fitted AKC - << initalGuess[2] << ";" // fitted S0 value + << S0 << ";" // 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] << ";"; // 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 &b /*target bValue*/) +::calculateSignalFromLsfCoeffs(vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & S0, 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); + //const double & S0 = lsfCoeffs(i,2); 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 0d9d504862..debbebce9e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h @@ -1,118 +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 double & S0); - void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & b); + void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & S0, 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