diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index 24c760e790..fe9d79e4a2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,356 +1,324 @@ /*=================================================================== 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) +::MultiShellAdcAverageReconstructionImageFilter() { this->SetNumberOfRequiredInputs( 1 ); - this->SetNumberOfThreads(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 - // 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; + const IndicesVector currentShell = it->second; unsigned int SHMaxOrder = 12; while( ((SHMaxOrder+1)*(SHMaxOrder+2)/2) > currentShell.size()) + { 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(); + //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(); + //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); - - - + //MITK_INFO << "shellInterpolationMatrix " << shellInterpolationMatrix.rows() << " x " << shellInterpolationMatrix.cols(); //- save interpolationSHBasis + m_ShellInterpolationMatrixVector.push_back(shellInterpolationMatrix); } - } 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); + MITK_INFO << m_WeightsVector.back(); } - // calculate average b-Value for target b-Value [bVal_t] - // MITK_INFO << "overall referenceImage avareg"; - - - // 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 it = m_BValueMap.begin(); - it++; // skip bZeroImages corresponding to 0-bValue + BValueMap::iterator ittt = m_BValueMap.begin(); + ittt++; // skip bZeroImages corresponding to 0-bValue m_TargetBValue = 0; - while(it!=m_BValueMap.end()) + while(ittt!=m_BValueMap.end()) { - m_TargetBValue += it->first; - it++; + m_TargetBValue += ittt->first; + ittt++; } m_TargetBValue /= (double)(m_BValueMap.size()-1); 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; + << " ReferenceImages: " << m_BValueMap.at(0.0).size() << 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(); // 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()-1); // create nx1 targetSignalVector vnl_vector SignalVector(m_allDirectionsSize); // ** walking over each Voxel while(!iit.IsAtEnd()) { InputPixelType b = iit.Get(); - - for(unsigned int i = 0 ; i < BZeroIndices.size(); i++) + 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]]); - MITK_INFO <<"RawSignal: "<< InterpVector; - + //MITK_INFO << "InputVector(RAWSIGNAL:"<< shellIndex <<")"; + //MITK_INFO << InterpVector ; //- normalization of the raw Signal S_S0Normalization(InterpVector, BZeroAverage); - MITK_INFO <<"Normalized: "<< InterpVector; - + //MITK_INFO << "InputVector(NORMALIZED:"<< shellIndex <<")"; + //MITK_INFO << InterpVector ; //- interpolate the Signal if necessary using corresponding interpolationSHBasis - if(m_Interpolation) - SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; - else - SignalVector = InterpVector; - MITK_INFO <<"Interpolated: "<< SignalVector; + + SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; + + //MITK_INFO << "SignalVector(NORMALIZED:"<< shellIndex <<")"; + //MITK_INFO << SignalVector ; // - ADC calculation for the signalVector calculateAdcFromSignal(SignalVector, shellIterator->first); - MITK_INFO << "ADCVector: " << SignalVector; - + //MITK_INFO << "SignalVector(SIGNaL->ADC:"<< shellIndex <<")"; + //MITK_INFO << SignalVector ; //- weight the signal - //if(m_Interpolation){ - // const double shellWeight = m_WeightsVector.at(shellIndex); - // SignalVector *= shellWeight; - //} + SignalVector *= m_WeightsVector.at(shellIndex); //- save the (interpolated) ShellSignalVector as the ith column in the SignalMatrix SignalMatrix.set_column(shellIndex, SignalVector); shellIterator++; shellIndex++; - //MITK_INFO << SignalMatrix; + } - exit(0); - // 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; + //MITK_INFO << "SignalVector(AVERAGEDADC)"; + //MITK_INFO << SignalVector ; + // recalculate signal from ADC calculateSignalFromAdc(SignalVector, m_TargetBValue, BZeroAverage); + //MITK_INFO << "SignalVector(ADC->SIGNAL)"; + //MITK_INFO << SignalVector ; for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); - //MITK_INFO << "ADCAveragedSignal: " << out; + //MITK_INFO << "NewSignal: " << out; //MITK_INFO << "-------------"; oit.Set(out); ++oit; ++iit; } - - - // 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++) + for(unsigned int i = 0; i < vec.size(); i++){ vec[i] /= S0; + if(vec[i]>1.0) vec[i] = 1; + if(vec[i]<0.0) vec[i] = 0; + } } 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 & referenceSignal) { 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 39441a4f95..4b828e7e05 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,119 +1,117 @@ /*=================================================================== 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, const double & referenceSignal); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions BValueMap m_BValueMap; double m_BValue; 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 diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx index afb1972fad..4778589dae 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,430 +1,430 @@ /*=================================================================== 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 "itkImageRegionIterator.h" #include "itkImageRegionConstIterator.h" #include "mitkImageCast.h" template mitk::DiffusionImage::DiffusionImage() : m_VectorImage(0), m_Directions(0), m_OriginalDirections(0), m_B_Value(-1.0), m_VectorImageAdaptor(0) { MeasurementFrameType mf; for(int i=0; i<3; i++) for(int j=0; j<3; j++) mf[i][j] = 0; for(int i=0; i<3; i++) mf[i][i] = 1; m_MeasurementFrame = mf; } template mitk::DiffusionImage::~DiffusionImage() { // Remove Observer for m_Directions RemoveDirectionsContainerObserver(); } template void mitk::DiffusionImage ::InitializeFromVectorImage() { if(!m_VectorImage || !m_Directions || m_B_Value==-1.0) { MITK_INFO << "DiffusionImage could not be initialized. Set all members first!" << std::endl; return; } // find bzero index int firstZeroIndex = -1; for(GradientDirectionContainerType::ConstIterator it = m_Directions->Begin(); it != m_Directions->End(); ++it) { firstZeroIndex++; GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); img->SetSpacing( m_VectorImage->GetSpacing() ); // Set the image spacing img->SetOrigin( m_VectorImage->GetOrigin() ); // Set the image origin img->SetDirection( m_VectorImage->GetDirection() ); // Set the image direction img->SetLargestPossibleRegion( m_VectorImage->GetLargestPossibleRegion()); img->SetBufferedRegion( m_VectorImage->GetLargestPossibleRegion() ); img->Allocate(); int vecLength = m_VectorImage->GetVectorLength(); InitializeByItk( img.GetPointer(), 1, vecLength ); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw = itw.Begin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr = itr.Begin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(firstZeroIndex)); ++itr; ++itw; } // init SetImportVolume(img->GetBufferPointer()); m_DisplayIndex = firstZeroIndex; MITK_INFO << "Diffusion-Image successfully initialized."; } template void mitk::DiffusionImage ::SetDisplayIndexForRendering(int displayIndex) { int index = displayIndex; int vecLength = m_VectorImage->GetVectorLength(); index = index > vecLength-1 ? vecLength-1 : index; if( m_DisplayIndex != index ) { typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); CastToItkImage(this, img); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw = itw.Begin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr = itr.Begin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(index)); ++itr; ++itw; } } m_DisplayIndex = index; } template bool mitk::DiffusionImage::AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision) { GradientDirectionType diff = g1 - g2; GradientDirectionType diff2 = g1 + g2; return diff.two_norm() < precision || diff2.two_norm() < precision; } template void mitk::DiffusionImage::CorrectDKFZBrokenGradientScheme(double precision) { GradientDirectionContainerType::Pointer directionSet = CalcAveragedDirectionSet(precision, m_Directions); if(directionSet->size() < 7) { MITK_INFO << "Too few directions, assuming and correcting DKFZ-bogus sequence details."; double v [7][3] = {{ 0, 0, 0 }, {-0.707057, 0, 0.707057 }, { 0.707057, 0, 0.707057 }, { 0, 0.707057, 0.707057 }, { 0, 0.707057, -0.707057 }, {-0.707057, 0.707057, 0 }, { 0.707057, 0.707057, 0 } }; int i=0; for(GradientDirectionContainerType::Iterator it = m_OriginalDirections->Begin(); it != m_OriginalDirections->End(); ++it) { it.Value().set(v[i++%7]); } ApplyMeasurementFrame(); } } template mitk::DiffusionImage::GradientDirectionContainerType::Pointer mitk::DiffusionImage::CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions) { // save old and construct new direction container GradientDirectionContainerType::Pointer newDirections = GradientDirectionContainerType::New(); // fill new direction container for(GradientDirectionContainerType::ConstIterator gdcitOld = directions->Begin(); gdcitOld != directions->End(); ++gdcitOld) { // already exists? bool found = false; for(GradientDirectionContainerType::ConstIterator gdcitNew = newDirections->Begin(); gdcitNew != newDirections->End(); ++gdcitNew) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { found = true; break; } } // if not found, add it to new container if(!found) { newDirections->push_back(gdcitOld.Value()); } } return newDirections; } template void mitk::DiffusionImage::AverageRedundantGradients(double precision) { GradientDirectionContainerType::Pointer newDirs = CalcAveragedDirectionSet(precision, m_Directions); GradientDirectionContainerType::Pointer newOriginalDirs = CalcAveragedDirectionSet(precision, m_OriginalDirections); // if sizes equal, we do not need to do anything in this function if(m_Directions->size() == newDirs->size() || m_OriginalDirections->size() == newOriginalDirs->size()) return; GradientDirectionContainerType::Pointer oldDirections = m_OriginalDirections; m_Directions = newDirs; m_OriginalDirections = newOriginalDirs; // new image typename ImageType::Pointer oldImage = m_VectorImage; m_VectorImage = ImageType::New(); m_VectorImage->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing m_VectorImage->SetOrigin( oldImage->GetOrigin() ); // Set the image origin m_VectorImage->SetDirection( oldImage->GetDirection() ); // Set the image direction m_VectorImage->SetLargestPossibleRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->SetVectorLength( m_Directions->size() ); m_VectorImage->SetBufferedRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->Allocate(); // average image data that corresponds to identical directions itk::ImageRegionIterator< ImageType > newIt(m_VectorImage, m_VectorImage->GetLargestPossibleRegion()); newIt.GoToBegin(); itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel typename ImageType::PixelType newVec; newVec.SetSize(m_Directions->size()); newVec.AllocateElements(m_Directions->size()); std::vector > dirIndices; for(GradientDirectionContainerType::ConstIterator gdcitNew = m_Directions->Begin(); gdcitNew != m_Directions->End(); ++gdcitNew) { dirIndices.push_back(std::vector(0)); for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin(); gdcitOld != oldDirections->End(); ++gdcitOld) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { //MITK_INFO << gdcitNew.Value() << " " << gdcitOld.Value(); dirIndices[gdcitNew.Index()].push_back(gdcitOld.Index()); } } } //int ind1 = -1; while(!newIt.IsAtEnd()) { // progress //typename ImageType::IndexType ind = newIt.GetIndex(); //ind1 = ind.m_Index[2]; // init new vector with zeros newVec.Fill(0.0); // the old voxel value with duplicates typename ImageType::PixelType oldVec = oldIt.Get(); for(unsigned int i=0; i void mitk::DiffusionImage::ApplyMeasurementFrame() { RemoveDirectionsContainerObserver(); m_Directions = GradientDirectionContainerType::New(); int c = 0; for(GradientDirectionContainerType::ConstIterator gdcit = m_OriginalDirections->Begin(); gdcit != m_OriginalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); vec = vec.pre_multiply(m_MeasurementFrame); m_Directions->InsertElement(c, vec); c++; } UpdateBValueList(); AddDirectionsContainerObserver(); } // returns number of gradients template int mitk::DiffusionImage::GetNumDirections() { int gradients = m_OriginalDirections->Size(); for (int i=0; iSize(); i++) if (GetB_Value(i)<=0) { gradients--; } return gradients; } // returns number of not diffusion weighted images template int mitk::DiffusionImage::GetNumB0() { int b0 = 0; for (int i=0; iSize(); i++) if (GetB_Value(i)<=0) { b0++; } return b0; } // returns a list of indices belonging to the not diffusion weighted images template typename mitk::DiffusionImage::IndicesVector mitk::DiffusionImage::GetB0Indices() { IndicesVector indices; for (int i=0; iSize(); i++) if (GetB_Value(i)<=0) { indices.push_back(i); } return indices; } template bool mitk::DiffusionImage::IsMultiBval() { int gradients = m_OriginalDirections->Size(); for (int i=0; i0 && std::fabs(m_B_Value-GetB_Value(i))>50) return true; return false; } template void mitk::DiffusionImage::UpdateBValueList() { m_B_ValueMap.clear(); GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = this->m_Directions->Begin(); gdcit != this->m_Directions->End(); ++gdcit) { float currentBvalue = std::floor(GetB_Value(gdcit.Index())); - double rounded = int((currentBvalue+7.5)/10)*10; + double rounded = int((currentBvalue+7.5)/100)*100; m_B_ValueMap[rounded].push_back(gdcit.Index()); } /* BValueMap::iterator it = m_B_ValueMap.begin(); for(;it != m_B_ValueMap.end(); it++) { MITK_INFO << it->first << " : " << it->second.size(); } */ } template float mitk::DiffusionImage::GetB_Value(int i) { if(i > m_Directions->Size()-1) return -1; if(m_Directions->ElementAt(i).one_norm() <= 0.0) { return 0; } else { double twonorm = m_Directions->ElementAt(i).two_norm(); return m_B_Value*twonorm*twonorm ; } } template void mitk::DiffusionImage::SetDirections(const std::vector > directions) { m_OriginalDirections = GradientDirectionContainerType::New(); for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); } this->ApplyMeasurementFrame(); } template void mitk::DiffusionImage::AddDirectionsContainerObserver() { // Add Modified Observer to invoke an UpdateBValueList by modifieng the DirectionsContainer (m_Directions) typedef DiffusionImage< TPixelType > Self; typedef itk::SimpleMemberCommand< Self > DCCommand ; typename DCCommand::Pointer command = DCCommand::New(); command->SetCallbackFunction(this, &Self::UpdateBValueList); } template void mitk::DiffusionImage::RemoveDirectionsContainerObserver() { if(m_Directions){ m_Directions->RemoveAllObservers(); } }