diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h index 3a6b34f27d..5e2fa898a1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h @@ -1,143 +1,155 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.h $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_TensorImageToDiffusionImageFilter_h_ #define _itk_TensorImageToDiffusionImageFilter_h_ #include "itkImageToImageFilter.h" #include namespace itk { template class TensorImageToDiffusionImageFilter : public ImageToImageFilter,3>, itk::VectorImage > { public: typedef TInputScalarType InputScalarType; typedef itk::DiffusionTensor3D InputPixelType; typedef itk::Image InputImageType; typedef typename InputImageType::RegionType InputImageRegionType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef typename BaselineImageType::RegionType BaselineImageRegionType; typedef TensorImageToDiffusionImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; itkTypeMacro (TensorImageToDiffusionImageFilter, ImageToImageFilter); itkStaticConstMacro (ImageDimension, unsigned int, OutputImageType::ImageDimension); itkNewMacro (Self); typedef Vector GradientType; typedef std::vector GradientListType; /** Manually Set/Get a list of gradients */ void SetGradientList(const GradientListType list) { m_GradientList = list; this->Modified(); } GradientListType GetGradientList(void) const {return m_GradientList;} void SetBValue( const double& bval) { m_BValue = bval; } + + /** + * @brief Set an external baseline image for signal generation (optional) + * + * An option to enforce a specific baseline image. If none provided (default) the filter uses + * the itk::TensorToL2NormImageFilter to generate the modelled baseline image. + */ + void SetExternalBaselineImage( typename BaselineImageType::Pointer bimage) + { + m_BaselineImage = bimage; + } + itkSetMacro(Min, OutputScalarType); itkSetMacro(Max, OutputScalarType); protected: TensorImageToDiffusionImageFilter() { m_BValue = 1.0; m_BaselineImage = 0; m_Min = 0.0; m_Max = 10000.0; - }; - ~TensorImageToDiffusionImageFilter(){}; + } + + ~TensorImageToDiffusionImageFilter(){} void PrintSelf (std::ostream& os, Indent indent) const { Superclass::PrintSelf (os, indent); } void BeforeThreadedGenerateData( void ); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); //void GenerateData(); - private: TensorImageToDiffusionImageFilter (const Self&); void operator=(const Self&); GradientListType m_GradientList; double m_BValue; typename BaselineImageType::Pointer m_BaselineImage; OutputScalarType m_Min; OutputScalarType m_Max; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkTensorImageToDiffusionImageFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx index f8427db3aa..4029ca4a0b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx @@ -1,229 +1,211 @@ /*=================================================================== 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. ===================================================================*/ /*========================================================================= 2 3 Program: Tensor ToolKit - TTK 4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $ 5 Language: C++ 6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ 7 Version: $Revision: 68 $ 8 9 Copyright (c) INRIA 2010. All rights reserved. 10 See LICENSE.txt for details. 11 12 This software is distributed WITHOUT ANY WARRANTY; without even 13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 14 PURPOSE. See the above copyright notices for more information. 15 16 =========================================================================*/ #ifndef _itk_TensorImageToDiffusionImageFilter_txx_ #define _itk_TensorImageToDiffusionImageFilter_txx_ #endif #include "itkTensorImageToDiffusionImageFilter.h" #include "itkTensorToL2NormImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include #include namespace itk { - //template - //void - // TensorImageToDiffusionImageFilter - // ::GenerateData() - //{ - // // Call a method that can be overriden by a subclass to allocate - // // memory for the filter's outputs - // this->AllocateOutputs(); - - // // Call a method that can be overridden by a subclass to perform - // // some calculations prior to splitting the main computations into - // // separate threads - // this->BeforeThreadedGenerateData(); - - // // Set up the multithreaded processing - // ThreadStruct str; - // str.Filter = this; - - // this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads()); - // this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str); - - // // multithread the execution - // this->GetMultiThreader()->SingleMethodExecute(); - - // // Call a method that can be overridden by a subclass to perform - // // some calculations after all the threads have completed - // this->AfterThreadedGenerateData(); - - //} - - template - void - TensorImageToDiffusionImageFilter - ::BeforeThreadedGenerateData() - { +template +void +TensorImageToDiffusionImageFilter +::BeforeThreadedGenerateData() +{ - if( m_GradientList.size()==0 ) - { - throw itk::ExceptionObject (__FILE__,__LINE__,"Error: gradient list is empty, cannot generate DWI."); - } + if( m_GradientList->Size()==0 ) + { + throw itk::ExceptionObject (__FILE__,__LINE__,"Error: gradient list is empty, cannot generate DWI."); + } + if( m_BaselineImage.IsNull() ) + { // create a B0 image by taking the norm of the tensor field * scale: typedef itk::TensorToL2NormImageFilter > - TensorToL2NormFilterType; + TensorToL2NormFilterType; typename TensorToL2NormFilterType::Pointer myFilter1 = TensorToL2NormFilterType::New(); myFilter1->SetInput (this->GetInput()); try { myFilter1->Update(); } catch (itk::ExceptionObject &e) { std::cerr << e; return; } typename itk::RescaleIntensityImageFilter< itk::Image, BaselineImageType>::Pointer rescaler= - itk::RescaleIntensityImageFilter, BaselineImageType>::New(); + itk::RescaleIntensityImageFilter, BaselineImageType>::New(); rescaler->SetOutputMinimum ( m_Min ); rescaler->SetOutputMaximum ( m_Max ); rescaler->SetInput ( myFilter1->GetOutput() ); try { rescaler->Update(); } catch (itk::ExceptionObject &e) { std::cerr << e; return; } m_BaselineImage = rescaler->GetOutput(); + } - typename OutputImageType::Pointer outImage = OutputImageType::New(); - outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing - outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin - outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction - outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); - outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); - outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); - outImage->SetVectorLength(m_GradientList.size()); - outImage->Allocate(); + typename OutputImageType::Pointer outImage = OutputImageType::New(); + outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing + outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin + outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction + outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); + outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); + outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); + outImage->SetVectorLength(m_GradientList->Size()); + outImage->Allocate(); + + this->SetNumberOfRequiredOutputs (1); + this->SetNthOutput (0, outImage); + +} + +template +void +TensorImageToDiffusionImageFilter +::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId ) +{ + typedef ImageRegionIterator IteratorOutputType; + typedef ImageRegionConstIterator IteratorInputType; + typedef ImageRegionConstIterator IteratorBaselineType; - this->SetNumberOfRequiredOutputs (1); - this->SetNthOutput (0, outImage); + unsigned long numPixels = outputRegionForThread.GetNumberOfPixels(); + unsigned long step = numPixels/100; + unsigned long progress = 0; - } + IteratorOutputType itOut (this->GetOutput(), outputRegionForThread); + IteratorInputType itIn (this->GetInput(), outputRegionForThread); + IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread); - template - void - TensorImageToDiffusionImageFilter - ::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId ) + if( threadId==0 ) { + this->UpdateProgress (0.0); + } - typedef ImageRegionIterator IteratorOutputType; - typedef ImageRegionConstIterator IteratorInputType; - typedef ImageRegionConstIterator IteratorBaselineType; - - unsigned long numPixels = outputRegionForThread.GetNumberOfPixels(); - unsigned long step = numPixels/100; - unsigned long progress = 0; - IteratorOutputType itOut (this->GetOutput(), outputRegionForThread); - IteratorInputType itIn (this->GetInput(), outputRegionForThread); - IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread); + while(!itIn.IsAtEnd()) + { - if( threadId==0 ) + if( this->GetAbortGenerateData() ) { - this->UpdateProgress (0.0); + throw itk::ProcessAborted(__FILE__,__LINE__); } + InputPixelType T = itIn.Get(); - while(!itIn.IsAtEnd()) - { - - if( this->GetAbortGenerateData() ) - { - throw itk::ProcessAborted(__FILE__,__LINE__); - } + BaselinePixelType b0 = itB0.Get(); - InputPixelType T = itIn.Get(); + OutputPixelType out; + out.SetSize(m_GradientList->Size()); + out.Fill(0); - BaselinePixelType b0 = itB0.Get(); - - OutputPixelType out; - out.SetSize(m_GradientList.size()); - out.Fill(0); - - if( b0 > 0) + if( b0 > 0) + { + for( unsigned int i=0; iSize()-1; i++) { - for( unsigned int i=0; iat(i+1); - InputPixelType S; - S[0] = g[0]*g[0]; - S[1] = g[1]*g[0]; - S[2] = g[2]*g[0]; - S[3] = g[1]*g[1]; - S[4] = g[2]*g[1]; - S[5] = g[2]*g[2]; + // normalize vector so the following computations work + const double twonorm = g.two_norm(); + GradientType gn = g.normalize(); - double res = - T[0]*S[0] + T[1]*S[1] + T[2]*S[2] + - T[1]*S[1] + T[3]*S[3] + T[4]*S[4] + - T[2]*S[2] + T[4]*S[4] + T[5]*S[5]; + InputPixelType S; + S[0] = gn[0]*gn[0]; + S[1] = gn[1]*gn[0]; + S[2] = gn[2]*gn[0]; + S[3] = gn[1]*gn[1]; + S[4] = gn[2]*gn[1]; + S[5] = gn[2]*gn[2]; - // check for corrupted tensor - if (res>=0) - out[i] = static_cast( 1.0*b0*exp ( -1.0 * m_BValue * res ) ); + const double res = + T[0]*S[0] + + 2 * T[1]*S[1] + T[3]*S[3] + + 2 * T[2]*S[2] + 2 * T[4]*S[4] + T[5]*S[5]; + + // check for corrupted tensor + if (res>=0) + { + // estimate the bvalue from the base value and the norm of the gradient + // - because of this estimation the vector have to be normalized beforehand + // otherwise the modelled signal is wrong ( i.e. not scaled properly ) + const double bval = m_BValue * twonorm; + out[i+1] = static_cast( 1.0 * b0 * exp ( -1.0 * bval * res ) ); } } + } - out[m_GradientList.size()-1] = b0; + out[0] = b0; - itOut.Set(out); + itOut.Set(out); - if( threadId==0 && step>0) + if( threadId==0 && step>0) + { + if( (progress%step)==0 ) { - if( (progress%step)==0 ) - { - this->UpdateProgress ( double(progress)/double(numPixels) ); - } + this->UpdateProgress ( double(progress)/double(numPixels) ); } + } - ++progress; - ++itB0; - ++itIn; - ++itOut; + ++progress; + ++itB0; + ++itIn; + ++itOut; - } + } - if( threadId==0 ) - { - this->UpdateProgress (1.0); - } + if( threadId==0 ) + { + this->UpdateProgress (1.0); } +} } // end of namespace