diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h index f014bc9215..7ccfb4b56c 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h @@ -1,116 +1,140 @@ /*=================================================================== 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 +#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; + 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; + // 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); + mitkClassMacro( DiffusionImage, Image ); + itkNewMacro(Self); - void AverageRedundantGradients(double precision); + void AverageRedundantGradients(double precision); - GradientDirectionContainerType::Pointer CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions); + GradientDirectionContainerType::Pointer CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions); - void CorrectDKFZBrokenGradientScheme(double precision); + 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); + 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); + GradientDirectionContainerType::Pointer GetDirections() { return m_Directions; } + void SetDirections( GradientDirectionContainerType::Pointer directions ) { this->m_Directions = directions; } + void SetDirections(const std::vector > directions); - GradientDirectionContainerType::Pointer GetOriginalDirections() { return m_OriginalDirections; } - void SetOriginalDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; this->ApplyMeasurementFrame(); } - void SetOriginalDirections(const std::vector > directions); + GradientDirectionContainerType::Pointer GetOriginalDirections() { return m_OriginalDirections; } + void SetOriginalDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; this->ApplyMeasurementFrame(); } + void SetOriginalDirections(const std::vector > directions); - MeasurementFrameType GetMeasurementFrame() { return m_MeasurementFrame; } - void SetMeasurementFrame( MeasurementFrameType mFrame ) { this->m_MeasurementFrame = mFrame; this->ApplyMeasurementFrame(); } + 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(); + bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); + int GetNumDirections(); + int GetNumB0(); - float GetB_Value(int i); - bool IsMultiBval(); - void UpdateBValueList(); + float GetB_Value(int i); + bool IsMultiBval(); + void UpdateBValueList(); - IndicesVector GetB0Indices(); + void AddObserver() { + // Add Observer for m_Directions + typedef DiffusionImage< TPixelType > Self; - itkGetMacro(B_Value, float); - itkSetMacro(B_Value, float); + typedef itk::SimpleMemberCommand< Self > DCCommand ; + + typename DCCommand::Pointer command = DCCommand::New(); + + command->SetCallbackFunction(this, &Self::UpdateBValueList); + + m_DirectionsObserverTag = m_Directions->AddObserver(itk::ModifiedEvent(), command); + } + + void RemoveObserver() { + // Remove Observer for m_Directions + m_Directions->RemoveObserver(m_DirectionsObserverTag ); + } + + IndicesVector GetB0Indices(); + + itkGetMacro(B_Value, float); + itkSetMacro(B_Value, float); + + BValueMap GetB_ValueMap(){ return m_B_ValueMap; } - 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; + 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; + + + + unsigned long m_DirectionsObserverTag; }; } // 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 1ccd9b185a..c882c09ce8 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,473 +1,477 @@ /*=================================================================== 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; + + AddObserver(); } template mitk::DiffusionImage::~DiffusionImage() { - + RemoveObserver(); } 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 ); //for(int i=0; i 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());//, 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 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 } }; 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_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) { dirIndices.push_back(std::vector(0)); for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin(); 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) { 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 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() { + + MITK_INFO << " TEST123| MODIFIED OOOWWEEEE"; 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() ); } }