diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h index 06f9cee942..f014bc9215 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h @@ -1,163 +1,116 @@ /*=================================================================== 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 __mitkDiffusionImage__h #define __mitkDiffusionImage__h #include "mitkImage.h" #include "itkVectorImage.h" #include "itkVectorImageToImageAdaptor.h" #include namespace mitk { /** * \brief this class encapsulates diffusion volumes (vectorimages not * yet supported by mitkImage) */ template class DiffusionImage : public Image { public: typedef TPixelType PixelType; typedef typename itk::VectorImage ImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef itk::VectorImageToImageAdaptor< TPixelType, 3 > AdaptorType; typedef vnl_matrix_fixed< double, 3, 3 > MeasurementFrameType; + // BValue Map + // key := b-Value + // value := indicesVector (containing corresponding gradient directions for a b-Value-Shell + typedef std::vector< unsigned int > IndicesVector; + typedef std::map< double , IndicesVector > BValueMap; + mitkClassMacro( DiffusionImage, Image ); itkNewMacro(Self); - //void SetRequestedRegionToLargestPossibleRegion(); - //bool RequestedRegionIsOutsideOfTheBufferedRegion(); - //virtual bool VerifyRequestedRegion(); - //void SetRequestedRegion(itk::DataObject *data); void AverageRedundantGradients(double precision); GradientDirectionContainerType::Pointer CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions); void CorrectDKFZBrokenGradientScheme(double precision); - typename ImageType::Pointer GetVectorImage() - { return m_VectorImage; } - void SetVectorImage(typename ImageType::Pointer image ) - { this->m_VectorImage = image; } + typename ImageType::Pointer GetVectorImage() { return m_VectorImage; } + void SetVectorImage(typename ImageType::Pointer image ) { this->m_VectorImage = image; } void InitializeFromVectorImage(); void SetDisplayIndexForRendering(int displayIndex); - GradientDirectionContainerType::Pointer GetDirections() - { return m_Directions; } - void SetDirections( GradientDirectionContainerType::Pointer directions ) - { this->m_Directions = directions; } - void SetDirections(const std::vector > directions) - { - m_Directions = GradientDirectionContainerType::New(); - for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); - } - } - GradientDirectionContainerType::Pointer GetOriginalDirections() - { return m_OriginalDirections; } - void SetOriginalDirections( GradientDirectionContainerType::Pointer directions ) - { this->m_OriginalDirections = directions; this->ApplyMeasurementFrame(); } - void SetOriginalDirections(const std::vector > directions) - { - m_OriginalDirections = GradientDirectionContainerType::New(); - for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); - } - this->ApplyMeasurementFrame(); - } - - MeasurementFrameType GetMeasurementFrame() - { return m_MeasurementFrame; } - void SetMeasurementFrame( MeasurementFrameType mFrame ) - { this->m_MeasurementFrame = mFrame; this->ApplyMeasurementFrame(); } + GradientDirectionContainerType::Pointer GetDirections() { return m_Directions; } + void SetDirections( GradientDirectionContainerType::Pointer directions ) { this->m_Directions = directions; } + void SetDirections(const std::vector > directions); - itkGetMacro(B_Value, float); - itkSetMacro(B_Value, float); + GradientDirectionContainerType::Pointer GetOriginalDirections() { return m_OriginalDirections; } + void SetOriginalDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; this->ApplyMeasurementFrame(); } + void SetOriginalDirections(const std::vector > directions); - float 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 ; - } - } + MeasurementFrameType GetMeasurementFrame() { return m_MeasurementFrame; } + void SetMeasurementFrame( MeasurementFrameType mFrame ) { this->m_MeasurementFrame = mFrame; this->ApplyMeasurementFrame(); } bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); - int GetNumDirections(); int GetNumB0(); - std::vector GetB0Indices(); + + float GetB_Value(int i); bool IsMultiBval(); + void UpdateBValueList(); + + IndicesVector GetB0Indices(); - void GetBvalueList(std::map > * bValMap) - { - if(bValMap != 0){ - 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; - (*bValMap)[rounded].push_back(gdcit.Index()); - } - }else - { - MITK_ERROR << "Call method with 0 parameter"; - } - } + itkGetMacro(B_Value, float); + itkSetMacro(B_Value, float); + BValueMap GetB_ValueMap(){ return m_B_ValueMap; } protected: DiffusionImage(); virtual ~DiffusionImage(); void ApplyMeasurementFrame(); typename ImageType::Pointer m_VectorImage; GradientDirectionContainerType::Pointer m_Directions; GradientDirectionContainerType::Pointer m_OriginalDirections; float m_B_Value; typename AdaptorType::Pointer m_VectorImageAdaptor; int m_DisplayIndex; MeasurementFrameType m_MeasurementFrame; + BValueMap m_B_ValueMap; }; } // namespace mitk #include "mitkDiffusionImage.txx" #endif /* __mitkDiffusionImage__h */ diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx index 53eec36178..1ccd9b185a 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,417 +1,473 @@ /*=================================================================== 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) + : 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() { } 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) + 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 ); //for(int i=0; i itw (img, img->GetLargestPossibleRegion() ); - itw = itw.Begin(); + itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); + itw = itw.Begin(); - itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); - itr = itr.Begin(); + itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); + itr = itr.Begin(); - while(!itr.IsAtEnd()) - { - itw.Set(itr.Get().GetElement(firstZeroIndex)); - ++itr; - ++itw; - } + while(!itr.IsAtEnd()) + { + itw.Set(itr.Get().GetElement(firstZeroIndex)); + ++itr; + ++itw; + } - // init - SetImportVolume(img->GetBufferPointer());//, 0, 0, CopyMemory); - //SetVolume( img->GetBufferPointer(), i ); + // init + SetImportVolume(img->GetBufferPointer());//, 0, 0, CopyMemory); + //SetVolume( img->GetBufferPointer(), i ); //} 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::RequestedRegionIsOutsideOfTheBufferedRegion() //{ // return false; //} // //template //void mitk::DiffusionImage::SetRequestedRegion(itk::DataObject * /*data*/) //{ //} // //template //void mitk::DiffusionImage::SetRequestedRegionToLargestPossibleRegion() //{ //} // //template //bool mitk::DiffusionImage::VerifyRequestedRegion() //{ // return true; //} //template //void mitk::DiffusionImage::DuplicateIfSingleSlice() //{ // // new image // typename ImageType::Pointer oldImage = m_Image; // m_Image = ImageType::New(); // m_Image->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing // m_Image->SetOrigin( oldImage->GetOrigin() ); // Set the image origin // m_Image->SetDirection( oldImage->GetDirection() ); // Set the image direction // typename ImageType::RegionType region = oldImage->GetLargestPossibleRegion(); // if(region.GetSize(0) == 1) // region.SetSize(0,3); // if(region.GetSize(1) == 1) // region.SetSize(1,3); // if(region.GetSize(2) == 1) // region.SetSize(2,3); // m_Image->SetLargestPossibleRegion( region ); // m_Image->SetVectorLength( m_Directions->size() ); // m_Image->SetBufferedRegion( region ); // m_Image->Allocate(); // // // average image data that corresponds to identical directions // itk::ImageRegionIterator< ImageType > newIt(m_Image, region); // newIt.GoToBegin(); // itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); // oldIt.GoToBegin(); // // while(!newIt.IsAtEnd()) // { // newIt.Set(oldIt.Get()); // ++newIt; // ++oldIt; // if(oldIt.IsAtEnd()) // oldIt.GoToBegin(); // } // //} template bool mitk::DiffusionImage::AreAlike(GradientDirectionType g1, - GradientDirectionType g2, - double precision) + GradientDirectionType g2, + double precision) { GradientDirectionType diff = g1 - g2; return diff.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 } }; + {{ 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 != 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) + gdcitOld != directions->End(); ++gdcitOld) { // already exists? bool found = false; for(GradientDirectionContainerType::ConstIterator gdcitNew = newDirections->Begin(); - gdcitNew != newDirections->End(); ++gdcitNew) + 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); + 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_Directions; GradientDirectionContainerType::Pointer oldOriginalDirections = 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) + gdcitNew != m_Directions->End(); ++gdcitNew) { dirIndices.push_back(std::vector(0)); for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin(); - gdcitOld != oldDirections->End(); ++gdcitOld) + gdcitOld != oldDirections->End(); ++gdcitOld) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { 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() { m_Directions = GradientDirectionContainerType::New(); int c = 0; for(GradientDirectionContainerType::ConstIterator gdcit = m_OriginalDirections->Begin(); - gdcit != m_OriginalDirections->End(); ++gdcit) + gdcit != m_OriginalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); vec = vec.pre_multiply(m_MeasurementFrame); m_Directions->InsertElement(c, vec); c++; } } // 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 -std::vector mitk::DiffusionImage::GetB0Indices() +typename mitk::DiffusionImage::IndicesVector mitk::DiffusionImage::GetB0Indices() { - std::vector indices; + 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() +{ + if(m_B_ValueMap.size() != 0) + { + 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; + m_B_ValueMap[rounded].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(); + return m_B_Value*twonorm*twonorm ; + } +} + +template +void mitk::DiffusionImage::SetOriginalDirections(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::SetDirections(const std::vector > directions) +{ + m_Directions = GradientDirectionContainerType::New(); + for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); + } +}