diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h index 06f9cee942..a7883718d9 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h @@ -1,163 +1,124 @@ /*=================================================================== 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; - - 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; } - - 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(); } - - itkGetMacro(B_Value, float); - itkSetMacro(B_Value, float); - - 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 ; - } - } - - bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); - - int GetNumDirections(); - int GetNumB0(); - std::vector GetB0Indices(); - bool IsMultiBval(); - - 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"; - } - } + 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 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; } + + 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 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(); } + + bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); + int GetNumDirections(); + int GetNumB0(); + + float GetB_Value(int i); + bool IsMultiBval(); + void UpdateBValueList(); + + IndicesVector GetB0Indices(); + + itkGetMacro(B_Value, float); + itkSetMacro(B_Value, float); + + BValueMap GetB_ValueMap(){ return m_B_ValueMap; } + + void AddDirectionsContainerObserver(); + void RemoveDirectionsContainerObserver(); 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; + 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 53eec36178..fdccc5023f 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,417 +1,512 @@ /*=================================================================== 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() { - + // Remove Observer for m_Directions + m_Directions->RemoveObserver(m_DirectionsObserverTag ); } 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() { + RemoveDirectionsContainerObserver(); + 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++; } + + 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 -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() +{ + 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()); + } + + /* + 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::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) +{ + + RemoveDirectionsContainerObserver(); + + m_Directions = GradientDirectionContainerType::New(); + for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); + } + + UpdateBValueList(); + AddDirectionsContainerObserver(); +} + +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(); + } +} + diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp index f4fd5c7c4c..ab313319e9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp @@ -1,441 +1,459 @@ /*=================================================================== 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. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkPreprocessingView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include // itk includes #include "itkTimeProbe.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "itkBrainMaskExtractionImageFilter.h" #include "itkCastImageFilter.h" #include "itkVectorContainer.h" #include // mitk includes #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include #include const std::string QmitkPreprocessingView::VIEW_ID = -"org.mitk.views.preprocessing"; + "org.mitk.views.preprocessing"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; QmitkPreprocessingView::QmitkPreprocessingView() -: QmitkFunctionality(), -m_Controls(NULL), -m_MultiWidget(NULL), -m_DiffusionImage(NULL) + : QmitkFunctionality(), + m_Controls(NULL), + m_MultiWidget(NULL), + m_DiffusionImage(NULL) { } QmitkPreprocessingView::QmitkPreprocessingView(const QmitkPreprocessingView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } QmitkPreprocessingView::~QmitkPreprocessingView() { } void QmitkPreprocessingView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkPreprocessingViewControls; m_Controls->setupUi(parent); this->CreateConnections(); m_Controls->m_MeasurementFrameTable->horizontalHeader()->setResizeMode(QHeaderView::Stretch); m_Controls->m_MeasurementFrameTable->verticalHeader()->setResizeMode(QHeaderView::Stretch); } } void QmitkPreprocessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkPreprocessingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkPreprocessingView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonAverageGradients), SIGNAL(clicked()), this, SLOT(AverageGradients()) ); connect( (QObject*)(m_Controls->m_ButtonExtractB0), SIGNAL(clicked()), this, SLOT(ExtractB0()) ); connect( (QObject*)(m_Controls->m_ButtonBrainMask), SIGNAL(clicked()), this, SLOT(BrainMask()) ); connect( (QObject*)(m_Controls->m_ModifyMeasurementFrame), SIGNAL(clicked()), this, SLOT(DoApplyMesurementFrame()) ); connect( (QObject*)(m_Controls->m_ReduceGradientsButton), SIGNAL(clicked()), this, SLOT(DoReduceGradientDirections()) ); connect( (QObject*)(m_Controls->m_ShowGradientsButton), SIGNAL(clicked()), this, SLOT(DoShowGradientDirections()) ); connect( (QObject*)(m_Controls->m_MirrorGradientToHalfSphereButton), SIGNAL(clicked()), this, SLOT(DoHalfSphereGradientDirections()) ); } } void QmitkPreprocessingView::OnSelectionChanged( std::vector nodes ) { bool foundDwiVolume = false; m_DiffusionImage = NULL; m_SelectedDiffusionNodes = mitk::DataStorage::SetOfObjects::New(); // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast*>(node->GetData()) ) { foundDwiVolume = true; m_DiffusionImage = dynamic_cast*>(node->GetData()); m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_SelectedDiffusionNodes->push_back(node); } } m_Controls->m_ButtonBrainMask->setEnabled(foundDwiVolume); m_Controls->m_ButtonAverageGradients->setEnabled(foundDwiVolume); m_Controls->m_ButtonExtractB0->setEnabled(foundDwiVolume); m_Controls->m_CheckExtractAll->setEnabled(foundDwiVolume); m_Controls->m_ModifyMeasurementFrame->setEnabled(foundDwiVolume); m_Controls->m_MeasurementFrameTable->setEnabled(foundDwiVolume); m_Controls->m_ReduceGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_ShowGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_MirrorGradientToHalfSphereButton->setEnabled(foundDwiVolume); if (foundDwiVolume) { vnl_matrix_fixed< double, 3, 3 > mf = m_DiffusionImage->GetMeasurementFrame(); for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); item->setTextAlignment(Qt::AlignCenter | Qt::AlignVCenter); item->setText(QString::number(mf.get(r,c))); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } - m_Controls->m_GradientsLabel->setText(QString::number(m_DiffusionImage->GetNumDirections())); - if (m_DiffusionImage->IsMultiBval()) - m_Controls->m_BvalLabel->setText("Acquisition with multiple b-values!"); - else - m_Controls->m_BvalLabel->setText(QString::number(m_DiffusionImage->GetB_Value())); + typedef mitk::DiffusionImage::BValueMap BValueMap; + typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; + + BValueMap bValMap = m_DiffusionImage->GetB_ValueMap(); + BValueMapIterator it = bValMap.begin(); + m_Controls->m_BvalueTable->clear(); + m_Controls->m_BvalueTable->setRowCount(bValMap.size() ); + QStringList headerList; + headerList << "b-Value" << "Number of gradients"; + m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); + int i = 0 ; + for(;it != bValMap.end(); it++) + { + m_Controls->m_BvalueTable->setItem(i,0,new QTableWidgetItem(QString::number(it->first))); + m_Controls->m_BvalueTable->setItem(i,1,new QTableWidgetItem(QString::number(it->second.size()))); + i++; + } + } else { for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } - m_Controls->m_GradientsLabel->setText("-"); - m_Controls->m_BvalLabel->setText("-"); + m_Controls->m_BvalueTable->clear(); + m_Controls->m_BvalueTable->setRowCount(1); + QStringList headerList; + headerList << "b-Value" << "Number of gradients"; + m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); + m_Controls->m_BvalueTable->setItem(0,0,new QTableWidgetItem("-")); + m_Controls->m_BvalueTable->setItem(0,1,new QTableWidgetItem("-")); m_Controls->m_DiffusionImageLabel->setText("-"); } } void QmitkPreprocessingView::Activated() { QmitkFunctionality::Activated(); } void QmitkPreprocessingView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkPreprocessingView::DoHalfSphereGradientDirections() { if (m_DiffusionImage.IsNull()) return; GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetOriginalDirections(); for (int j=0; jSize(); j++) if (gradientContainer->at(j)[0]<0) gradientContainer->at(j) = -gradientContainer->at(j); } void QmitkPreprocessingView::DoApplyMesurementFrame() { if (m_DiffusionImage.IsNull()) return; vnl_matrix_fixed< double, 3, 3 > mf; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); if (!item) return; mf[r][c] = item->text().toDouble(); } m_DiffusionImage->SetMeasurementFrame(mf); } void QmitkPreprocessingView::DoShowGradientDirections() { if (m_DiffusionImage.IsNull()) return; GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetOriginalDirections(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); for (int j=0; jSize(); j++) { mitk::Point3D p; vnl_vector_fixed< double, 3 > v = gradientContainer->at(j); if (fabs(v[0])>0.001 || fabs(v[1])>0.001 || fabs(v[2])>0.001) { p[0] = v[0]; p[1] = v[1]; p[2] = v[2]; pointset->InsertPoint(j, p); } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(pointset); node->SetName("gradient directions"); node->SetProperty("pointsize", mitk::FloatProperty::New(0.05)); node->SetProperty("color", mitk::ColorProperty::New(1,0,0)); GetDefaultDataStorage()->Add(node); } void QmitkPreprocessingView::DoReduceGradientDirections() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage DiffusionImageType; typedef itk::ReduceDirectionGradientsFilter FilterType; GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetOriginalDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(m_DiffusionImage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(m_Controls->m_ReduceGradientsBox->value()); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(m_DiffusionImage->GetB_Value()); image->SetDirections(filter->GetGradientDirections()); image->SetOriginalDirections(filter->GetGradientDirections()); image->SetMeasurementFrame(m_DiffusionImage->GetMeasurementFrame()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); imageNode->SetName("reduced_image"); GetDefaultDataStorage()->Add(imageNode); } void QmitkPreprocessingView::ExtractB0() { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; int nrFiles = m_SelectedDiffusionNodes->size(); if (!nrFiles) return; // call the extraction withou averaging if the check-box is checked if( this->m_Controls->m_CheckExtractAll->isChecked() ) { DoExtractBOWithoutAveraging(); return; } mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = - static_cast( - (*itemiter)->GetData()); + static_cast( + (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::DoExtractBOWithoutAveraging() { // typedefs typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; typedef itk::B0ImageExtractionToSeparateImageFilter< short, short> FilterType; // check number of selected objects, return if empty int nrFiles = m_SelectedDiffusionNodes->size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes->end() ); std::vector< mitk::DataNode::Pointer > nodes; while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = - static_cast( - (*itemiter)->GetData()); + static_cast( + (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0_ALL")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::AverageGradients() { int nrFiles = m_SelectedDiffusionNodes->size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = - static_cast*>( - (*itemiter)->GetData()); + static_cast*>( + (*itemiter)->GetData()); vols->AverageRedundantGradients(m_Controls->m_Blur->value()); ++itemiter; } } void QmitkPreprocessingView::BrainMask() { int nrFiles = m_SelectedDiffusionNodes->size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes->end() ); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = - static_cast*>( - (*itemiter)->GetData()); + static_cast*>( + (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); typedef itk::CastImageFilter, itk::Image > CastFilterType; CastFilterType::Pointer castfilter = CastFilterType::New(); castfilter->SetInput(filter->GetOutput()); typedef itk::BrainMaskExtractionImageFilter MaskFilterType; MaskFilterType::Pointer maskfilter = MaskFilterType::New(); maskfilter->SetInput(castfilter->GetOutput()); maskfilter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( maskfilter->GetOutput() ); mitkImage->SetVolume( maskfilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_Mask")); GetDefaultDataStorage()->Add(node); ++itemiter; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui index c67d78669e..7fe02fe94a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui @@ -1,529 +1,555 @@ QmitkPreprocessingViewControls 0 0 389 - 913 + 984 0 0 - true + false QmitkPreprocessingViewControls + + true + Data Diffusion Image: - + + + 0 + 0 + + Info - - - QFormLayout::AllNonFixedFieldsGrow - + - - - Number of Gradients: - - - - - - - - - - - Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - - - - - - b-Value: + + + + 0 + 0 + - - - - - - - + + Qt::ScrollBarAsNeeded - - Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + Qt::ScrollBarAlwaysOff + + true + + + 100 + + + true + + + false + + + true + + + + b-Value + + + + + Number of gradients + + - + false Generate pointset displaying the gradient vectors. Show gradients + + + + Qt::Horizontal + + + + 40 + 20 + + + + Reduce size 0 70 Multiple acquistions of one gradient direction can be averaged. Due to rounding errors, similar gradients often differ in the last decimal positions. The Merge radius allows to average them anyway by taking into account all directions within a certain radius. true QFrame::NoFrame QFrame::Raised 0 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. 6 2.000000000000000 0.000100000000000 0.001000000000000 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Merge radius false Average redundant gradients false Mirror all gradients around one axis. Mirror gradients to half sphere Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 0 0 New number of gradients: 1 30 false Retain only the specified number of gradient directions and according image volumes. The retained directions are spread equally over the half sphere. Reduce number of gradients Non diffusion weighted image 0 30 Average and extract all images that were acquired without diffusion weighting. true false Extract B0 Create a 3D+t data set containing all b0 images as timesteps Extract all B0 without averaging Brain mask false Estimate binary brain mask 0 0 Measurment frame Qt::Horizontal 40 20 false 0 0 10 10 IBeamCursor true Qt::ScrollBarAlwaysOff Qt::ScrollBarAlwaysOff true false false true true false true true New Row New Row New Row New Column New Column New Column false Apply new mesurement frame Qt::Vertical 20 40 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp index 9c33c46767..96d51602d2 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp @@ -1,1419 +1,1419 @@ /*=================================================================== 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 "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include #include #include #include #include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.h" //#include "itkFreeWaterEliminationFilter.h" #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkDiffusionImageMapper.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include #include const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; using namespace berry; struct TrSelListener : ISelectionListener { berryObjectMacro(TrSelListener); TrSelListener(QmitkTensorReconstructionView* view) { m_View = view; } void DoSelectionChanged(ISelection::ConstPointer selection) { /* // if(!m_View->IsVisible()) // return; // save current selection in member variable m_View->m_CurrentSelection = selection.Cast(); // do something with the selected items if(m_View->m_CurrentSelection) { bool foundDwiVolume = false; bool foundTensorVolume = false; // iterate selection for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin(); i != m_View->m_CurrentSelection->End(); ++i) { // extract datatree node if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); // only look at interesting types if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { foundDwiVolume = true; } // only look at interesting types if(QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0) { foundTensorVolume = true; } } } m_View->m_Controls->m_ItkReconstruction->setEnabled(foundDwiVolume); m_View->m_Controls->m_TeemReconstruction->setEnabled(foundDwiVolume); m_View->m_Controls->m_ReconstructionWithCorrection->setEnabled(foundDwiVolume); m_View->m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_View->m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); m_View->m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_View->m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); m_View->m_Controls->m_PerSliceView->setEnabled(foundDwiVolume && foundTensorVolume); } */ } void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection) { // check, if selection comes from datamanager if (part) { QString partname(part->GetPartName().c_str()); if(partname.compare("Datamanager")==0) { // apply selection DoSelectionChanged(selection); } } } QmitkTensorReconstructionView* m_View; }; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); } QmitkTensorReconstructionView::QmitkTensorReconstructionView(const QmitkTensorReconstructionView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); QStringList items; items << "LLS (Linear Least Squares)" << "MLE (Maximum Likelihood)" << "NLS (Nonlinear Least Squares)" << "WLS (Weighted Least Squares)"; m_Controls->m_TensorEstimationTeemEstimationMethodCombo->addItems(items); m_Controls->m_TensorEstimationTeemEstimationMethodCombo->setCurrentIndex(0); m_Controls->m_TensorEstimationManualThreashold->setChecked(false); m_Controls->m_TensorEstimationTeemSigmaEdit->setText("NaN"); m_Controls->m_TensorEstimationTeemNumItsSpin->setValue(1); m_Controls->m_TensorEstimationTeemFuzzyEdit->setText("0.0"); m_Controls->m_TensorEstimationTeemMinValEdit->setText("1.0"); m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(true); m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(true); m_Controls->m_TensorsToDWIBValueEdit->setText("1000"); Advanced1CheckboxClicked(); Advanced2CheckboxClicked(); Advanced3CheckboxClicked(); TeemCheckboxClicked(); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_TeemToggle->setVisible(false); #endif // define data type for combobox //m_Controls->m_ImageSelector->SetDataStorage( this->GetDefaultDataStorage() ); //m_Controls->m_ImageSelector->SetPredicate( mitk::NodePredicateDataType::New("DiffusionImage") ); } } void QmitkTensorReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkTensorReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_TeemToggle), SIGNAL(clicked()), this, SLOT(TeemCheckboxClicked()) ); connect( (QObject*)(m_Controls->m_ItkReconstruction), SIGNAL(clicked()), this, SLOT(ItkReconstruction()) ); connect( (QObject*)(m_Controls->m_TeemReconstruction), SIGNAL(clicked()), this, SLOT(TeemReconstruction()) ); connect( (QObject*)(m_Controls->m_ReconstructionWithCorrection), SIGNAL(clicked()), this, SLOT(ReconstructionWithCorrection()) ); connect( (QObject*)(m_Controls->m_TensorEstimationTeemEstimationMethodCombo), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_Advanced2), SIGNAL(clicked()), this, SLOT(Advanced2CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_Advanced3), SIGNAL(clicked()), this, SLOT(Advanced3CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorEstimationManualThreashold), SIGNAL(clicked()), this, SLOT(ManualThresholdClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); connect( (QObject*)(m_Controls->m_PerSliceView), SIGNAL(pointSelected(int, int)), this, SLOT(ResidualClicked(int, int)) ); } } void QmitkTensorReconstructionView::ResidualClicked(int slice, int volume) { // Use image coord to reset crosshair // Find currently selected diffusion image // Update Label // to do: This position should be modified in order to skip B0 volumes that are not taken into account // when calculating residuals // Find the diffusion image mitk::DiffusionImage* diffImage; mitk::DataNode::Pointer correctNode; mitk::Geometry3D* geometry; if (m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); geometry = diffImage->GetGeometry(); // Remember the node whose display index must be updated correctNode = mitk::DataNode::New(); correctNode = m_DiffusionImage; } if(diffImage != NULL) { typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; GradientDirectionContainerType::Pointer dirs = diffImage->GetDirections(); for(int i=0; iSize() && i<=volume; i++) { GradientDirectionType grad = dirs->ElementAt(i); // check if image is b0 weighted if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { volume++; } } QString pos = "Volume: "; pos.append(QString::number(volume)); pos.append(", Slice: "); pos.append(QString::number(slice)); m_Controls->m_PositionLabel->setText(pos); if(correctNode) { int oldDisplayVal; correctNode->GetIntProperty("DisplayChannel", oldDisplayVal); std::string oldVal = QString::number(oldDisplayVal).toStdString(); std::string newVal = QString::number(volume).toStdString(); correctNode->SetIntProperty("DisplayChannel",volume); correctNode->SetSelected(true); this->FirePropertyChanged("DisplayChannel", oldVal, newVal); correctNode->UpdateOutputInformation(); mitk::Point3D p3 = m_MultiWidget->GetCrossPosition(); itk::Index<3> ix; geometry->WorldToIndex(p3, ix); // ix[2] = slice; mitk::Vector3D vec; vec[0] = ix[0]; vec[1] = ix[1]; vec[2] = slice; mitk::Vector3D v3New; geometry->IndexToWorld(vec, v3New); mitk::Point3D origin = geometry->GetOrigin(); mitk::Point3D p3New; p3New[0] = v3New[0] + origin[0]; p3New[1] = v3New[1] + origin[1]; p3New[2] = v3New[2] + origin[2]; m_MultiWidget->MoveCrossToPosition(p3New); m_MultiWidget->RequestUpdate(); } } } void QmitkTensorReconstructionView::TeemCheckboxClicked() { m_Controls->groupBox_3->setVisible(m_Controls-> m_TeemToggle->isChecked()); } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Advanced2CheckboxClicked() { bool check = m_Controls-> m_Advanced2->isChecked(); m_Controls->frame_2->setVisible(check); } void QmitkTensorReconstructionView::Advanced3CheckboxClicked() { bool check = m_Controls-> m_Advanced3->isChecked(); m_Controls->frame_6->setVisible(check); } void QmitkTensorReconstructionView::ManualThresholdClicked() { m_Controls->m_TensorReconstructionThreasholdEdit_2->setEnabled( m_Controls->m_TensorEstimationManualThreashold->isChecked()); } void QmitkTensorReconstructionView::Activated() { QmitkFunctionality::Activated(); } void QmitkTensorReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkTensorReconstructionView::MethodChoosen(int method) { m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(method==3); m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(method==3); } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::DiffusionImage::Pointer diffImage = mitk::DiffusionImage::New(); TensorImageType::Pointer tensorImage; std::string nodename; if(m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); } else return; if(m_TensorImage.IsNotNull()) { mitk::TensorImage* mitkVol; mitkVol = static_cast(m_TensorImage->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); m_TensorImage->GetStringProperty("name", nodename); } else return; typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListType gradientList; mitk::DiffusionImage::GradientDirectionContainerType* gradients = diffImage->GetDirections(); // Copy gradients vectors from gradients to gradientList for(int i=0; iSize(); i++) { mitk::DiffusionImage::GradientDirectionType vec = gradients->at(i); itk::Vector grad; grad[0] = vec[0]; grad[1] = vec[1]; grad[2] = vec[2]; gradientList.push_back(grad); } // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue(diffImage->GetB_Value()); filter->SetGradientList(gradientList); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(500); filter->Update(); // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(diffImage->GetB_Value()); image->SetDirections(gradientList); image->SetOriginalDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); - std::vector b0Indices = image->GetB0Indices(); + std::vector b0Indices = image->GetB0Indices(); typedef itk::ResidualImageFilter ResidualImageFilterType; ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput(diffImage->GetVectorImage()); residualFilter->SetSecondDiffusionImage(image->GetVectorImage()); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); int size = lookupTable->GetTable()->GetSize(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residual Image"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDefaultDataStorage()->Add(resNode); m_MultiWidget->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } // Draw Graph for volumes per slice in the QGraphicsView std::vector< std::vector > outliersPerSlice = residualFilter->GetOutliersPerSlice(); int xSize = outliersPerSlice.size(); if(xSize == 0) { return; } int ySize = outliersPerSlice[0].size(); // Find maximum in outliersPerSlice double maxOutlier= 0.0; for(int i=0; imaxOutlier) { maxOutlier = outliersPerSlice[i][j]; } } } // Create some QImage QImage qImage(xSize, ySize, QImage::Format_RGB32); QImage legend(1, 256, QImage::Format_RGB32); QRgb value; vtkSmartPointer lookup = vtkSmartPointer::New(); lookup->SetTableRange(0.0, maxOutlier); lookup->Build(); reversedlookupTable->SetTableRange(0, maxOutlier); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookup->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } // Fill qImage for(int i=0; iMapValue(out); int r, g, b; r = _rgba[0]; g = _rgba[1]; b = _rgba[2]; value = qRgb(r, g, b); qImage.setPixel(i,j,value); } } for(int i=0; i<256; i++) { double* rgba = lookup->GetTableValue(i); int r, g, b; r = rgba[0]*255; g = rgba[1]*255; b = rgba[2]*255; value = qRgb(r, g, b); legend.setPixel(0,255-i,value); } QString upper = QString::number(maxOutlier, 'g', 3); upper.append(" %"); QString lower = QString::number(0.0); lower.append(" %"); m_Controls->m_UpperLabel->setText(upper); m_Controls->m_LowerLabel->setText(lower); QGraphicsScene* scene = new QGraphicsScene; QGraphicsScene* scene2 = new QGraphicsScene; QPixmap pixmap(QPixmap::fromImage(qImage)); QGraphicsPixmapItem *item = new QGraphicsPixmapItem( pixmap, 0, scene); item->scale(10.0, 3.0); QPixmap pixmap2(QPixmap::fromImage(legend)); QGraphicsPixmapItem *item2 = new QGraphicsPixmapItem( pixmap2, 0, scene2); item2->scale(20.0, 1.0); m_Controls->m_PerSliceView->SetResidualPixmapItem(item); m_Controls->m_PerSliceView->setScene(scene); m_Controls->m_LegendView->setScene(scene2); m_Controls->m_PerSliceView->show(); m_Controls->m_PerSliceView->repaint(); m_Controls->m_LegendView->setHorizontalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->setVerticalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->show(); m_Controls->m_LegendView->repaint(); } void QmitkTensorReconstructionView::ItkReconstruction() { Reconstruct(0); } void QmitkTensorReconstructionView::TeemReconstruction() { Reconstruct(1); } void QmitkTensorReconstructionView::ReconstructionWithCorrection() { Reconstruct(2); } void QmitkTensorReconstructionView::Reconstruct(int method) { if(method == 0) ItkTensorReconstruction(m_DiffusionImages); if(method == 1) TeemTensorReconstruction(m_DiffusionImages); if(method == 2) { TensorReconstructionWithCorr(m_DiffusionImages); } } void QmitkTensorReconstructionView::TensorReconstructionWithCorr (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { typedef mitk::DiffusionImage DiffusionImageType; DiffusionImageType* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MBI_INFO << "Tensor reconstruction with correction for negative eigenvalues"; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::TensorReconstructionWithEigenvalueCorrectionFilter< DiffusionPixelType, TTensorPixelType > ReconstructionFilter; float b0Threshold = m_Controls->m_ReconstructionThreshold->text().toFloat(); ReconstructionFilter::Pointer reconFilter = ReconstructionFilter::New(); reconFilter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); reconFilter->SetBValue(vols->GetB_Value()); reconFilter->SetB0Threshold(b0Threshold); reconFilter->Update(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer outputTensorImg = reconFilter->GetOutput(); typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(outputTensorImg, outputTensorImg->GetRequestedRegion()); tensorIt.GoToBegin(); int negatives = 0; while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iInitializeByItk( outputTensorImg.GetPointer() ); image->SetVolume( outputTensorImg->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dti"); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); // Corrected diffusion image typedef itk::VectorImage ImageType; ImageType::Pointer correctedVols = reconFilter->GetVectorImage(); DiffusionImageType::Pointer correctedDiffusion = DiffusionImageType::New(); correctedDiffusion->SetVectorImage(correctedVols); correctedDiffusion->SetDirections(vols->GetDirections()); correctedDiffusion->SetB_Value(vols->GetB_Value()); correctedDiffusion->SetOriginalDirections(vols->GetDirections()); correctedDiffusion->InitializeFromVectorImage(); mitk::DataNode::Pointer diffNode=mitk::DataNode::New(); diffNode->SetData( correctedDiffusion ); QString diffname; diffname = diffname.append(nodename.c_str()); diffname = diffname.append("corrDiff"); SetDefaultNodeProperties(diffNode, diffname.toStdString()); nodes.push_back(diffNode); // B0 mask as used in tensorreconstructionwithcorrectionfilter typedef itk::Image MaskImageType; MaskImageType::Pointer mask = reconFilter->GetMask(); mitk::Image::Pointer mitkMask; mitk::CastToMitkImage(mask, mitkMask); mitk::DataNode::Pointer maskNode=mitk::DataNode::New(); maskNode->SetData( mitkMask ); QString maskname; maskname = maskname.append(nodename.c_str()); maskname = maskname.append("_mask"); SetDefaultNodeProperties(maskNode, maskname.toStdString()); nodes.push_back(maskNode); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkTensorReconstructionView::ItkTensorReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MBI_INFO << "Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); tensorReconstructionFilter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); tensorReconstructionFilter->SetBValue(vols->GetB_Value()); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreasholdEdit->text().toFloat() ); tensorReconstructionFilter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s."; // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); // TensorType2 tensor2; /* for(int i=0; i SymEigenSystemType; SymEigenSystemType eig (tensor2.GetVnlMatrix()); for(unsigned int i=0; iInitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dti"); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkTensorReconstructionView::TeemTensorReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MBI_INFO << "Teem Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Teem Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef mitk::TeemDiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); tensorReconstructionFilter->SetInput( vols ); if(!m_Controls->m_TensorEstimationTeemSigmaEdit->text().contains(QString("NaN"))) tensorReconstructionFilter->SetSigma( m_Controls->m_TensorEstimationTeemSigmaEdit->text().toFloat() ); switch(m_Controls->m_TensorEstimationTeemEstimationMethodCombo->currentIndex()) { // items << "LLS (Linear Least Squares)" //<< "MLE (Maximum Likelihood)" //<< "NLS (Nonlinear Least Squares)" //<< "WLS (Weighted Least Squares)"; case 0: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS); break; case 1: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsMLE); break; case 2: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsNLS); break; case 3: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsWLS); break; default: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS); } tensorReconstructionFilter->SetNumIterations( m_Controls->m_TensorEstimationTeemNumItsSpin->value() ); if(m_Controls->m_TensorEstimationManualThreashold->isChecked()) tensorReconstructionFilter->SetConfidenceThreshold( m_Controls->m_TensorReconstructionThreasholdEdit_2->text().toDouble() ); tensorReconstructionFilter->SetConfidenceFuzzyness( m_Controls->m_TensorEstimationTeemFuzzyEdit->text().toFloat() ); tensorReconstructionFilter->SetMinPlausibleValue( m_Controls->m_TensorEstimationTeemMinValEdit->text().toDouble() ); tensorReconstructionFilter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ; // TENSORS TO DATATREE mitk::DataNode::Pointer node2=mitk::DataNode::New(); node2->SetData( tensorReconstructionFilter->GetOutputItk() ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dtix"); SetDefaultNodeProperties(node2, newname.toStdString()); nodes.push_back(node2); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); //node->SetProperty( "opacity", mitk::FloatProperty::New(1.0f) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } //node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) ); //node->SetProperty( "use color", mitk::BoolProperty::New( true ) ); //node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) ); //node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ); //node->SetProperty( "layer", mitk::IntProperty::New(0)); //node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ); //node->SetOpacity(1.0f); //node->SetColor(1.0,1.0,1.0); //node->SetVisibility(true); //node->SetProperty( "IsTensorVolume", mitk::BoolProperty::New( true ) ); //mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); //mitk::LevelWindow levelwindow; //// levelwindow.SetAuto( image ); //levWinProp->SetLevelWindow( levelwindow ); //node->GetPropertyList()->SetProperty( "levelwindow", levWinProp ); //// add a default rainbow lookup table for color mapping //if(!node->GetProperty("LookupTable")) //{ // mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New(); // vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable(); // vtkLut->SetHueRange(0.6667, 0.0); // vtkLut->SetTableRange(0.0, 20.0); // vtkLut->Build(); // mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New(); // mitkLutProp->SetLookupTable(mitkLut); // node->SetProperty( "LookupTable", mitkLutProp ); //} //if(!node->GetProperty("binary")) // node->SetProperty( "binary", mitk::BoolProperty::New( false ) ); //// add a default transfer function //mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New(); //node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) ); //// set foldername as string property //mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name ); //node->SetProperty( "name", nameProp ); void QmitkTensorReconstructionView::TensorsToDWI() { DoTensorsToDWI(m_TensorImages); } void QmitkTensorReconstructionView::TensorsToQbi() { for (int i=0; isize(); i++) { mitk::DataNode::Pointer tensorImageNode = m_TensorImages->at(i); MITK_INFO << "starting Q-Ball estimation"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::QBallImage::Pointer image = mitk::QBallImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(tensorImageNode->GetName().c_str()); newname = newname.append("_qbi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); } } void QmitkTensorReconstructionView::OnSelectionChanged( std::vector nodes ) { if ( !this->IsVisible() ) return; m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); bool foundDwiVolume = false; bool foundTensorVolume = false; m_Controls->m_DiffusionImageLabel->setText("-"); m_Controls->m_TensorImageLabel->setText("-"); m_DiffusionImage = NULL; m_TensorImage = NULL; // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNull()) continue; // only look at interesting types if(dynamic_cast*>(node->GetData())) { foundDwiVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_DiffusionImages->push_back(node); m_DiffusionImage = node; } else if(dynamic_cast(node->GetData())) { foundTensorVolume = true; m_Controls->m_TensorImageLabel->setText(node->GetName().c_str()); m_TensorImages->push_back(node); m_TensorImage = node; } } m_Controls->m_ItkReconstruction->setEnabled(foundDwiVolume); m_Controls->m_TeemReconstruction->setEnabled(foundDwiVolume); m_Controls->m_ReconstructionWithCorrection->setEnabled(foundDwiVolume); m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PerSliceView->setEnabled(foundDwiVolume && foundTensorVolume); } template std::vector > QmitkTensorReconstructionView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } // Add 0 vector for B0 itk::Vector v; v.Fill(0.0); retval.push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { std::string nodename; (*itemiter)->GetStringProperty("name", nodename); mitk::TensorImage* vol = static_cast((*itemiter)->GetData()); ++itemiter; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListType gradientList; switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); // DWI ESTIMATION clock.Start(); MBI_INFO << "DWI Estimation "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "DWI Estimation for %s", nodename.c_str()).toAscii()); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); //filter->SetNumberOfThreads(1); filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s."; // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(bVal); image->SetDirections(gradientList); image->SetOriginalDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } }