diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index f3e5d43528..85c37a2017 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,194 +1,225 @@ /*=================================================================== 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_MultiShellAdcAverageReconstructionImageFilter_cpp_ #define _itk_MultiShellAdcAverageReconstructionImageFilter_cpp_ #endif #define _USE_MATH_DEFINES #include "itkMultiShellAdcAverageReconstructionImageFilter.h" #include #include #include #include #include "mitkDiffusionFunctionCollection.h" namespace itk { template MultiShellAdcAverageReconstructionImageFilter ::MultiShellAdcAverageReconstructionImageFilter(): m_Interpolation(false) { this->SetNumberOfRequiredInputs( 1 ); } template void MultiShellAdcAverageReconstructionImageFilter ::BeforeThreadedGenerateData() { // test whether BvalueMap contains all necessary information if(m_BValueMap.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_BValue * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; m_BValueMap[bValueKey].push_back(gdcit.Index()); } } //# BValueMap contains no bZero --> itkException if(m_BValueMap.find(0.0) == m_BValueMap.end()) { MITK_INFO << "No ReferenceSignal (BZeroImages) found!"; itkExceptionMacro(<< "No ReferenceSignal (BZeroImages) found!"); } // test whether interpolation is necessary // - Gradeint directions on different shells are different m_Interpolation = mitk::gradients::CheckForDifferingShellDirections(m_BValueMap, m_OriginalGradientDirections.GetPointer()); // if INTERPOLATION necessary if(m_Interpolation) { for(BValueMap::const_iterator it = m_BValueMap.begin();it != m_BValueMap.end(); it++) { if((*it).first == 0.0) continue; // if any #ShellDirection < 15 --> itkException (No interpolation possible) if((*it).second.size() < 15){ MITK_INFO << "Abort: No interpolation possible Shell-" << (*it).first << " has less than 15 directions."; itkExceptionMacro(<<"No interpolation possible"); } } // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions IndicesVector allDirectionsIndicies = mitk::gradients::GetAllUniqueDirections(m_BValueMap, m_OriginalGradientDirections); // [sizeAllDirections] size of GradientContainer cointaining all unique directions - const int allDirectionsSize = allDirections.size(); - - /* for each shell - * - calculate maxShOrder - * - calculate Weights [Weigthing = shell_size / max_shell_size] - * - get TragetSHBasis using allDirectionsContainer - * - get ShellSHBasis using currentShellDirections - * - calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] - * - save interpolationSHBasis - */ + const int allDirectionsSize = allDirectionsIndicies.size(); + std::vector SHMaxOrderVector(m_BValueMap.size()-1); + std::vector WeightsVector(m_BValueMap.size()-1); + std::vector > ShellInterpolationMatrixVector(m_BValueMap.size()-1); + // for Weightings + unsigned int maxShellSize = 0; + + // for each shell + BValueMap::const_iterator it = m_BValueMap.begin(); + it++; //skip bZeroIndices + for(;it != m_BValueMap.end();it++) + { + //- calculate maxShOrder + IndicesVector currentShell = (*it).second; + unsigned int SHMaxOrder = 12; + while( ((SHMaxOrder+1)*(SHMaxOrder+2)/2) > currentShell.size()) + SHMaxOrder -= 2 ; + + //- save shell size + WeightsVector.push_back(currentShell.size()); + if(currentShell.size() > maxShellSize) + maxShellSize = currentShell.size(); + + //- get TragetSHBasis using allDirectionsContainer + vnl_matrix sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(allDirectionsIndicies, m_OriginalGradientDirections); + vnl_matrix TargetSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); + //- get ShellSHBasis using currentShellDirections + sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(currentShell, m_OriginalGradientDirections); + vnl_matrix ShellSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); + //- calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] + vnl_matrix_inverse invShellSHBasis(ShellSHBasis); + vnl_matrix shellInterpolationMatrix = TargetSHBasis * invShellSHBasis.inverse(); + ShellInterpolationMatrixVector.push_back(shellInterpolationMatrix); + //- save interpolationSHBasis + + } + + //- calculate Weights [Weigthing = shell_size / max_shell_size] + for(int i = 0 ; i < WeightsVector.size(); i++) + WeightsVector.at(i) /= maxShellSize; + } // calculate average b-Value for target b-Value [bVal_t] // calculate target bZero-Value [b0_t] MITK_INFO << "Input:" << std::endl << "GradientDirections: " << m_OriginalGradientDirections->Size() << std::endl << "Shells: " << (m_BValueMap.size() - 1) << std::endl << "ReferenceImages: " << m_BValueMap.at(0.0).size() << std::endl << "Interpolation: " << m_Interpolation; } template void MultiShellAdcAverageReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int /*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(); const int numShells = m_BValueMap.size()-1; BValueMap::iterator it = m_BValueMap.begin(); //std::vector adcVec = new Vector(numShells); */ // create empty nxm SignalMatrix containing n->signals/directions (in case of interpolation ~ sizeAllDirections otherwise the size of any shell) for m->shells // create nx1 targetSignalVector // ** walking over each Voxel /* for each shell in this Voxel * - get the RawSignal * - interpolate the Signal if necessary using corresponding interpolationSHBasis * - save the (interpolated) ShellSignalVector as the ith column in the SignalMatrix */ // normalize the signals in SignalMatrix [bZeroAverage????] [S/S0 = E] /* for each row in the SignalMatrix * - calculate for each rowentry the ADC [ln(E)/-b = ADC] * - weight the ith ADC entry with the corresponding shellWeighting * - average the Signal over the row [ADC_t] * - calculate target Signal using target b-Value [S_t = b0_t * e^(-bVal_t * ADC_t) * - Save S_t in targetSignalVector */ // outImageIterator set S_t // ** /* int vecLength; // initialize output image typename OutputImageType::Pointer 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( vecLength ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); MITK_INFO << "...done"; */ } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp index 20902c1114..66cd2d5e87 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp @@ -1,169 +1,205 @@ /*=================================================================== 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. ===================================================================*/ #include "mitkDiffusionFunctionCollection.h" #include #include "mitkVector.h" // for Windows #ifndef M_PI #define M_PI 3.14159265358979323846 #endif // Namespace ::SH #include #include #include // Namespace ::Gradients #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *cart) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = M_PI/2; phi = M_PI/2; } else { th = acos(z/rad); phi = atan2(y, x); } cart[0] = phi; cart[1] = th; cart[2] = rad; } double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i mitk::gradients::GetAllUniqueDirections(const std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; BValueMap::const_iterator mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); IndiciesVector::iterator containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } bool mitk::gradients::CheckForDifferingShellDirections(const std::map > & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { BValueMap::const_iterator mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ BValueMap::const_iterator mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } template double mitk::gradients::dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ) { double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm()); return result ; } + +vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, GradientDirectionContainerType::Pointer refGradientsContainer) +{ + + vnl_matrix Q(3, refShell.size()); + + for(unsigned int i = 0; i < refShell.size(); i++) + { + double x = refGradientsContainer->ElementAt(refShell[i]).normalize().get(0); + double y = refGradientsContainer->ElementAt(refShell[i]).normalize().get(1); + double z = refGradientsContainer->ElementAt(refShell[i]).normalize().get(2); + double cart[3]; + mitk::sh::Cart2Sph(x,y,z,cart); + Q(0,i) = cart[0]; + Q(1,i) = cart[1]; + Q(2,i) = cart[2]; + } + return Q; +} + +vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) +{ + vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); + for(unsigned int i=0; i< SHBasisOutput.rows(); i++) + for(int k = 0; k <= LOrder; k += 2) + for(int m =- k; m <= k; m++) + { + int j = ( k * k + k + 2 ) / 2 + m - 1; + double phi = QBallReference(0,i); + double th = QBallReference(1,i); + SHBasisOutput(i,j) = mitk::sh::Yj(m,k,th,phi); + } + return SHBasisOutput; +} diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h index 66f91ffdd1..8d71f85cb9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h @@ -1,61 +1,62 @@ /*=================================================================== 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 __mitkDiffusionFunctionCollection_h_ #define __mitkDiffusionFunctionCollection_h_ #include "DiffusionCoreExports.h" #include "vector.h" #include "map.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_vector_fixed.h" #include "itkVectorContainer.h" namespace mitk{ class DiffusionCore_EXPORT sh { public: static double factorial(int number); static void Cart2Sph(double x, double y, double z, double* cart); static double legendre0(int l); static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); static double Yj(int m, int k, double theta, double phi); - }; class DiffusionCore_EXPORT gradients { private: typedef std::vector IndiciesVector; typedef std::map BValueMap; typedef itk::VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; public: static std::vector GetAllUniqueDirections(const std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ); static bool CheckForDifferingShellDirections(const std::map > & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); + static vnl_matrix ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder); + static vnl_matrix ComputeSphericalFromCartesian(const IndiciesVector & refShell, GradientDirectionContainerType::Pointer refGradientsContainer); template static double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ); }; } #endif //__mitkDiffusionFunctionCollection_h_