diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx index af04974c9d..006e59d3ba 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,424 +1,424 @@ /*=================================================================== 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++; } UpdateBValueMap(); 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::UpdateBValueMap() { m_B_ValueMap.clear(); GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = this->m_Directions->Begin(); gdcit != this->m_Directions->End(); ++gdcit) - m_B_ValueMap[std::floor((GetB_Value(gdcit.Index())+0.5))].push_back(gdcit.Index()); + m_B_ValueMap[GetB_Value(gdcit.Index())].push_back(gdcit.Index()); } 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(); double bval = m_B_Value*twonorm*twonorm; if (bval<0) bval = ceil(bval - 0.5); else bval = floor(bval + 0.5); return bval; } } 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::UpdateBValueMap); } template void mitk::DiffusionImage::RemoveDirectionsContainerObserver() { if(m_Directions){ m_Directions->RemoveAllObservers(); } } diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h index cdff3c6495..b8a88facb7 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h @@ -1,62 +1,63 @@ /*=================================================================== 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 "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; + typedef vnl_vector_fixed GradientDirectionType; public: static std::vector GetAllUniqueDirections(const BValueMap &bValueMap, GradientDirectionContainerType *refGradientsContainer ); static bool CheckForDifferingShellDirections(const BValueMap &bValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); static vnl_matrix ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder); static vnl_matrix ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer); static mitk::gradients::GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType * origninalGradentcontainer); template static double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ); }; } #endif //__mitkDiffusionFunctionCollection_h_