diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h index 549bc1549d..0bb074566a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h @@ -1,89 +1,114 @@ /*=================================================================== 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 ) */ 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; + /** 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; + } + 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 dc55d58f6c..e311e6ad5b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx @@ -1,54 +1,159 @@ /*=================================================================== 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 >::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 = outputIt.GoToBegin(); + inputIt = 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