diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h index 566c4806ca..75430dc8b4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h @@ -1,145 +1,152 @@ /*=================================================================== 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 ITKSPLITDWIMAGEFILTER_H #define ITKSPLITDWIMAGEFILTER_H #include "itkImageToImageFilter.h" #include namespace itk { /** \class SplitDWImageFilter * * \brief Splits a DW-Image passed in as input into a 3D-t image where each volume coresponds to a * gradient image ( or the unweighted b0 ) * * Several applications require to get the gradient images as a separate volume, f.e. the registration for * head-motion correction. Also a reduction of the DW Image is possible when combined with its counterpart filter, * the \sa mitkImageToDiffusionImageSource, which can reinterpret a 3d+t (scalar) image into a diffusion weighted image. */ template< class TInputImagePixelType, class TOutputImagePixelType > class SplitDWImageFilter : public ImageToImageFilter< VectorImage, Image< TOutputImagePixelType, 4 > > { public: typedef SplitDWImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef TInputImagePixelType InputPixelType; typedef TOutputImagePixelType OutputPixelType; typedef ImageToImageFilter< VectorImage, Image< TOutputImagePixelType, 4 > > Superclass; /** typedefs from superclass */ typedef typename Superclass::InputImageType InputImageType; //typedef typename Superclass::OutputImageType OutputImageType; typedef Image< TOutputImagePixelType, 4 > OutputImageType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef std::vector< unsigned int > IndexListType; typedef std::map< unsigned int, IndexListType> BValueMapType; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(SplitDWImageFilter, ImageToImageFilter); /** * @brief Set the indices of the images to be extracted. * * */ void SetExtractIndices( IndexListType list) { m_IndexList = list; } /** * @brief Extract all images * * The same as setting SetExtractIndices( ) with [0,1,2,...,N] */ void SetExtractAll() { m_ExtractAllImages = true; } /** * @brief Selects only the weighted images with b-value above the given b_threshold to be extracted * * Setting b_threshold to 0 will do the same as \sa SetExtractAll. Please note that some images have no true * unweighted images as the minimal b-value is something like 5 so for extracting all. * * @note It will reorder the images! * * @param b_threshold the minimal b-value to be extracted * @param map the map with b-values to the corresponding image * * @sa GetIndexList */ void SetExtractAllAboveThreshold( double b_threshold, BValueMapType map); + /** + * @brief SetExtractSingleShell + * @param b_value b-value of the shell to be extracted + * @param tol the tolerance of the shell choice, i.e. all shells within [ b_value - tol, b_value + tol ] will be extracted + */ + void SetExtractSingleShell( double b_value, BValueMapType map, double tol ); + /** * @brief Returns the index list used for extraction * * The list is necessary for further processing, especially when a b-value threshold is used ( like in \sa SetExtractAllAboveThreshold ) * * @return The index list used during the extraction */ const IndexListType GetIndexList() const { return m_IndexList; } protected: SplitDWImageFilter(); virtual ~SplitDWImageFilter(){}; void GenerateData(); /** The dimension of the output does not match the dimension of the input hence we need to re-implement the CopyInformation method to avoid executing the default implementation which tries to copy the input information to the output */ virtual void CopyInformation( const DataObject *data); /** Override of the ProcessObject::GenerateOutputInformation() because of different dimensionality of the input and the output */ virtual void GenerateOutputInformation(); IndexListType m_IndexList; bool m_ExtractAllImages; }; } //end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkSplitDWImageFilter.txx" #endif #endif // ITKSPLITDWIMAGEFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx index 79a97ab912..0455285422 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx @@ -1,189 +1,212 @@ /*=================================================================== 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 __itkSplitDWImageFilter_txx #define __itkSplitDWImageFilter_txx #include "itkSplitDWImageFilter.h" #include #include template< class TInputImagePixelType, class TOutputImagePixelType> itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::SplitDWImageFilter() : m_IndexList(0), m_ExtractAllImages(true) { this->SetNumberOfRequiredInputs(1); } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::CopyInformation( const DataObject* /*data*/) { } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType >::GenerateOutputInformation() { } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::SetExtractAllAboveThreshold(double b_threshold, BValueMapType map) { m_ExtractAllImages = false; m_IndexList.clear(); // create the index list following the given threshold // iterate over the b-value map BValueMapType::const_iterator bvalueIt = map.begin(); while( bvalueIt != map.end() ) { // check threshold if (bvalueIt->first > b_threshold) { // the map contains an index container, this needs to be inserted into the // index list IndexListType::const_iterator listIt = bvalueIt->second.begin(); while( listIt != bvalueIt->second.end() ) { m_IndexList.push_back( *listIt ); ++listIt; } } ++bvalueIt; } } +template< class TInputImagePixelType, + class TOutputImagePixelType> +void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > +::SetExtractSingleShell(double b_value, BValueMapType map, double tol) +{ + m_ExtractAllImages = false; + m_IndexList.clear(); + + // create index list + BValueMapType::const_iterator bvalueIt = map.begin(); + while( bvalueIt != map.end() ) + { + IndexListType::const_iterator listIt = bvalueIt->second.begin(); + if( std::fabs( bvalueIt->first - b_value) < tol) + { + m_IndexList.push_back(*listIt); + ++listIt; + } + + ++bvalueIt; + } +} + template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType >::GenerateData() { if( m_IndexList.empty() && !m_ExtractAllImages ) { itkExceptionMacro(<<"The index list is empty and the option to extract all images is disabled. "); } // construct the index list (for each component) if( m_ExtractAllImages ) { m_IndexList.clear(); for( unsigned int i=0; i< this->GetInput()->GetNumberOfComponentsPerPixel(); i++) m_IndexList.push_back(i); } // declare the output image // this will have the b0 images stored as timesteps typename OutputImageType::Pointer outputImage = this->GetOutput(); //OutputImageType::New(); // allocate image with // - dimension: [DimX, DimY, DimZ, size(IndexList) ] // - spacing: old one, 1.0 time typename OutputImageType::SpacingType spacing; spacing.Fill(1); typename OutputImageType::PointType origin; origin.Fill(0); OutputImageRegionType outputRegion; // the spacing and origin corresponds to the respective values in the input image (3D) // the same for the region for (unsigned int i=0; i< 3; i++) { spacing[i] = (this->GetInput()->GetSpacing())[i]; origin[i] = (this->GetInput()->GetOrigin())[i]; outputRegion.SetSize(i, this->GetInput()->GetLargestPossibleRegion().GetSize()[i]); outputRegion.SetIndex(i, this->GetInput()->GetLargestPossibleRegion().GetIndex()[i]); } // number of timesteps = number of b0 images outputRegion.SetSize(3, m_IndexList.size()); outputRegion.SetIndex( 3, 0 ); // output image direction (4x4) typename OutputImageType::DirectionType outputDirection; //initialize to identity outputDirection.SetIdentity(); // get the input image direction ( 3x3 matrix ) typename InputImageType::DirectionType inputDirection = this->GetInput()->GetDirection(); for( unsigned int i=0; i< 3; i++) { outputDirection(0,i) = inputDirection(0,i); outputDirection(1,i) = inputDirection(1,i); outputDirection(2,i) = inputDirection(2,i); } outputImage->SetSpacing( spacing ); outputImage->SetOrigin( origin ); outputImage->SetDirection( outputDirection ); outputImage->SetRegions( outputRegion ); outputImage->Allocate(); // input iterator itk::ImageRegionConstIterator inputIt( this->GetInput(), this->GetInput()->GetLargestPossibleRegion() ); // we want to iterate separately over each 3D volume of the output image // so reset the regions last dimension outputRegion.SetSize(3,1); unsigned int currentTimeStep = 0; // for each index in the iterator list, extract the image and insert it as a new time step for(IndexListType::const_iterator indexIt = m_IndexList.begin(); indexIt != m_IndexList.end(); indexIt++) { // set the time step outputRegion.SetIndex(3, currentTimeStep); itk::ImageRegionIterator< OutputImageType> outputIt( outputImage.GetPointer(), outputRegion ); // iterate over the current b0 image and store it to corresponding place outputIt.GoToBegin(); inputIt.GoToBegin(); while( !outputIt.IsAtEnd() && !inputIt.IsAtEnd() ) { // the input vector typename InputImageType::PixelType vec = inputIt.Get(); outputIt.Set( vec[*indexIt]); ++outputIt; ++inputIt; } // increase time step currentTimeStep++; } } #endif // __itkSplitDWImageFilter_txx