diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index 022776dee5..69af5a8ce3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,241 +1,275 @@ /*=================================================================== 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()); + // [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(); // 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 = allDirectionsIndicies.size(); - std::vector SHMaxOrderVector(m_BValueMap.size()-1); - std::vector WeightsVector(m_BValueMap.size()-1); - std::vector > ShellInterpolationMatrixVector(m_BValueMap.size()-1); + 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 - IndicesVector currentShell = (*it).second; + const 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 sphericalCoordinates; + sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(m_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); + 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(int i = 0 ; i < WeightsVector.size(); i++) - WeightsVector.at(i) /= maxShellSize; - + for(;it != m_BValueMap.end(); it++) + m_WeightsVector.push_back(it->second.size() / maxShellSize); } - // calculate average b-Value for target b-Value [bVal_t] - IndicesVector BZeroIndices = m_BValueMap.at(0.0); + // 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) { - for(int i = 0 ; i < numberOfBZeroImages; i++) + 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); } - // calculate target bZero-Value [b0_t] - double BZeroAverage = 0.0; - for(int i = 0 ; i < numberOfBZeroImages; i++) - BZeroAverage += m_OriginalGradientDirections->ElementAt(BZeroIndices.at(i)); - BZeroAverage /= numberOfBZeroImages; + + // 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( m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer + outImage->Allocate(); 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 + // 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 // create nx1 targetSignalVector + + // TO DO + // ** 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 - */ + while(!iit.IsAtEnd()) + { + InputPixelType b = iit.Get(); + + for(int i = 0 ; i < BZeroIndices.size(); i++) + BZeroAverage += b[BZeroIndices[i]]; + BZeroAverage /= BZeroIndices.size(); + + OutputPixelType out; + out.AllocateElements(m_allDirectionsSize + 1); + out.Fill(0.0); + out.SetElement(0,BZeroAverage); + + /* 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/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h index 249a298e63..6893cca1bd 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,106 +1,115 @@ /*=================================================================== 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 ); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions BValueMap m_BValueMap; float m_BValue; float 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