diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index 94e57efd1c..e5d20a10a9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,349 +1,344 @@ /*=================================================================== 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 ); + //this->SetNumberOfThreads(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()); // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions m_allDirectionsIndicies = mitk::gradients::GetAllUniqueDirections(m_BValueMap, m_OriginalGradientDirections); // [sizeAllDirections] size of GradientContainer cointaining all unique directions m_allDirectionsSize = m_allDirectionsIndicies.size(); + //MITK_INFO << m_allDirectionsSize; + //exit(0); + + m_TargetGradientDirections = mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(m_BValueMap,m_OriginalGradientDirections); + //MITK_INFO << m_allDirectionsSize; + // 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"); } } m_ShellInterpolationMatrixVector.reserve(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 const IndicesVector currentShell = (*it).second; unsigned int SHMaxOrder = 12; while( ((SHMaxOrder+1)*(SHMaxOrder+2)/2) > currentShell.size()) SHMaxOrder -= 2 ; - //- save shell size - - if(currentShell.size() > maxShellSize) - maxShellSize = currentShell.size(); - //- 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(); m_ShellInterpolationMatrixVector.push_back(shellInterpolationMatrix); + + + //- save interpolationSHBasis } - m_WeightsVector.reserve(m_BValueMap.size()-1); - it = m_BValueMap.begin(); - it++; // skip ReferenceImages - //- calculate Weights [Weigthing = shell_size / max_shell_size] - for(;it != m_BValueMap.end(); it++) - m_WeightsVector.push_back(it->second.size() / maxShellSize); + + } + m_WeightsVector.reserve(m_BValueMap.size()-1); + BValueMap::const_iterator itt = m_BValueMap.begin(); + itt++; // skip ReferenceImages + //- calculate Weights [Weigthing = shell_size / max_shell_size] + unsigned int maxShellSize = 0; + for(;itt != m_BValueMap.end(); itt++){ + if(itt->second.size() > maxShellSize) + maxShellSize = itt->second.size(); + } + itt = m_BValueMap.begin(); + itt++; // skip ReferenceImages + for(;itt != m_BValueMap.end(); itt++){ + m_WeightsVector.push_back(itt->second.size() / (double)maxShellSize); } // calculate average b-Value for target b-Value [bVal_t] - const IndicesVector BZeroIndices = m_BValueMap.at(0.0); - const unsigned int numberOfBZeroImages = BZeroIndices.size(); - if(numberOfBZeroImages % (m_BValueMap.size()-1) == 0) - { - MITK_INFO << "referenceImage avareg for each shell"; - m_bZeroIndicesSplitVectors.reserve(m_BValueMap.size()-1); - const int stepWidth = BZeroIndices.size() / (m_BValueMap.size()-1); - IndicesVector::const_iterator it1 = BZeroIndices.begin(); - IndicesVector::const_iterator it2 = BZeroIndices.begin() + stepWidth; - for(int i = 0 ; i < m_BValueMap.size()-1; i++) - { - m_bZeroIndicesSplitVectors.push_back(IndicesVector(it1,it2)); - it1 += stepWidth; - it2 += stepWidth; - } - }else - { - MITK_INFO << "overall referenceImage avareg"; - m_bZeroIndicesSplitVectors.reserve(1); - m_bZeroIndicesSplitVectors.push_back(BZeroIndices); - } + // MITK_INFO << "overall referenceImage avareg"; + // initialize output image - typename OutputImageType::Pointer outImage = OutputImageType::New(); + 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( m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer + outImage->SetVectorLength( 1+m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); + BValueMap::iterator it = m_BValueMap.begin(); + it++; // skip bZeroImages corresponding to 0-bValue + m_TargetBValue = 0; + while(it!=m_BValueMap.end()) + { + m_TargetBValue += it->first; + it++; + } + m_TargetBValue /= (double)(m_BValueMap.size()-1); + - 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; + MITK_INFO << "Input:" << std::endl << 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 << std::endl; + MITK_INFO << "Output:" << std::endl << std::endl + << " OutImageVectorLength: " << outImage->GetVectorLength() << std::endl + << " TargetDirections: " << m_allDirectionsSize << std::endl + << " TargetBValue: " << m_TargetBValue << std::endl + << std::endl; } 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); - // calculate target bZero-Value [b0_t] const IndicesVector BZeroIndices = m_BValueMap[0.0]; double BZeroAverage = 0.0; // 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, m_BValueMap.size()); + vnl_matrix SignalMatrix(m_allDirectionsSize, m_BValueMap.size()-1); // create nx1 targetSignalVector vnl_vector SignalVector(m_allDirectionsSize); // ** walking over each Voxel while(!iit.IsAtEnd()) { InputPixelType b = iit.Get(); - for(int i = 0 ; i < BZeroIndices.size(); i++) + for(unsigned int i = 0 ; i < BZeroIndices.size(); i++) BZeroAverage += b[BZeroIndices[i]]; - BZeroAverage /= BZeroIndices.size(); + BZeroAverage /= (double)BZeroIndices.size(); OutputPixelType out = oit.Get(); - //out.AllocateElements(m_allDirectionsSize + 1); out.Fill(0.0); out.SetElement(0,BZeroAverage); BValueMap::const_iterator shellIterator = m_BValueMap.begin(); shellIterator++; //skip bZeroImages unsigned int shellIndex = 0; - bool multipleBValues = m_bZeroIndicesSplitVectors.size() > 1; - while(shellIterator != m_BValueMap.end()) { // reset Data SignalVector.fill(0.0); // - get the RawSignal const IndicesVector currentShell = shellIterator->second; - for(int i = 0 ; i < currentShell.size(); i++) + for(unsigned int i = 0 ; i < currentShell.size(); i++) SignalVector.put(i,b[currentShell[i]]); + //MITK_INFO <<"RawSignal: "<< SignalVector; + MITK_INFO << "möp" <first); + //MITK_INFO <<"ADC: "<< SignalVector; //- weight the signal if(m_Interpolation){ const double shellWeight = m_WeightsVector.at(shellIndex); SignalVector *= shellWeight; } //- save the (interpolated) ShellSignalVector as the ith column in the SignalMatrix SignalMatrix.set_column(shellIndex, SignalVector); shellIterator++; shellIndex++; + //MITK_INFO << SignalMatrix; } - - for(int i = 0 ; i < SignalMatrix.rows(); i++) + // MITK_INFO <<"finalSignalMatrix: " << SignalMatrix; + // ADC averaging Sum(ADC)/n + for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++) SignalVector.put(i,SignalMatrix.get_row(i).sum()); + SignalVector/= (double)(m_BValueMap.size()-1); + // MITK_INFO << "AveragedADC: " << SignalVector; + // recalculate signal from ADC + calculateSignalFromAdc(SignalVector, m_TargetBValue, BZeroAverage); - SignalVector/= SignalMatrix.cols(); - - calculateSignalFromAdc(SignalVector,,BZeroAverage) - - for(int i = 1 ; i < out.Size(); i ++) + for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); - + //MITK_INFO << "ADCAveragedSignal: " << out; + //MITK_INFO << "-------------"; oit.Set(out); ++oit; ++iit; } - - - - - // 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; this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); MITK_INFO << "...done"; */ } template void MultiShellAdcAverageReconstructionImageFilter ::S_S0Normalization( vnl_vector & vec, const double & S0 ) { for(unsigned int i = 0; i < vec.size(); i++) vec[i] /= S0; } template void MultiShellAdcAverageReconstructionImageFilter ::calculateAdcFromSignal( vnl_vector & vec, const double & bValue) { for(unsigned int i = 0; i < vec.size(); i++) vec[i] = log(vec[i])/-bValue; } template void MultiShellAdcAverageReconstructionImageFilter -::calculateSignalFromAdc( vnl_vector & vec, const double & bValue, const double & bZero) +::calculateSignalFromAdc(vnl_vector & vec, const double & bValue, const double & referenceSignal) { - for(unsigned int i = 0 ; i < vec.size(); i++) - vec[i] = bZero * exp(-bValue * vec[i]) + for(unsigned int i = 0 ; i < vec.size(); i++){ + //MITK_INFO << vec[i]; + //MITK_INFO << bValue; + //MITK_INFO << referenceSignal; + vec[i] = referenceSignal * exp((-bValue) * vec[i]); + } } - } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h index 8a0eb15dbb..39441a4f95 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,119 +1,119 @@ /*=================================================================== 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_MultiShellAdcAverageReconstructionImageFilter_h_ #define _itk_MultiShellAdcAverageReconstructionImageFilter_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 MultiShellAdcAverageReconstructionImageFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef MultiShellAdcAverageReconstructionImageFilter 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; itkGetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(TargetGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(TargetBValue, float) itkGetMacro(BValue,float) itkSetMacro(BValue,float) inline void SetOriginalBValueMap(BValueMap inp){m_BValueMap = inp;} protected: MultiShellAdcAverageReconstructionImageFilter(); ~MultiShellAdcAverageReconstructionImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads ); void S_S0Normalization( vnl_vector & vec, const double & S0 ); void calculateAdcFromSignal( vnl_vector & vec, const double & bValue); - void calculateSignalFromAdc( vnl_vector & vec, const double & bValue); + void calculateSignalFromAdc( vnl_vector & vec, 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_BValueMap; - float m_BValue; + double m_BValue; - float m_TargetBValue; + double m_TargetBValue; bool m_Interpolation; std::vector m_WeightsVector; std::vector > m_ShellInterpolationMatrixVector; std::vector m_bZeroIndicesSplitVectors; IndicesVector m_allDirectionsIndicies; unsigned int m_allDirectionsSize; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellAdcAverageReconstructionImageFilter.cpp" #endif #endif