1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
16 #ifndef __itkB0ImageExtractionToSeparateImageFilter_txx
17 #define __itkB0ImageExtractionToSeparateImageFilter_txx
19 #include "itkB0ImageExtractionToSeparateImageFilter.h"
21 template< class TInputImagePixelType,
22 class TOutputImagePixelType>
23 itk::B0ImageExtractionToSeparateImageFilter< TInputImagePixelType, TOutputImagePixelType >
24 ::B0ImageExtractionToSeparateImageFilter()
26 this->SetNumberOfRequiredInputs( 1 );
29 template< class TInputImagePixelType,
30 class TOutputImagePixelType>
31 void itk::B0ImageExtractionToSeparateImageFilter<
32 TInputImagePixelType, TOutputImagePixelType >::CopyInformation( const DataObject *)
37 template< class TInputImagePixelType,
38 class TOutputImagePixelType>
39 void itk::B0ImageExtractionToSeparateImageFilter<
40 TInputImagePixelType, TOutputImagePixelType >::GenerateOutputInformation()
45 template< class TInputImagePixelType,
46 class TOutputImagePixelType>
47 void itk::B0ImageExtractionToSeparateImageFilter<
48 TInputImagePixelType, TOutputImagePixelType >::GenerateData()
51 GradContainerIteratorType begin = m_Directions->Begin();
52 GradContainerIteratorType end = m_Directions->End();
54 // Find the indices of the b0 images in the diffusion image
55 std::vector<int> indices;
59 GradientDirectionType grad = begin->Value();
61 if(grad[0] == 0 && grad[1] == 0 && grad[2] == 0)
63 indices.push_back(index);
69 // declare the output image
70 // this will have the b0 images stored as timesteps
71 typename OutputImageType::Pointer outputImage = this->GetOutput();
72 //OutputImageType::New();
74 // get the input region
75 typename Superclass::InputImageType::RegionType inputRegion =
76 this->GetInput()->GetLargestPossibleRegion();
78 // allocate image with
79 // - dimension: [DimX, DimY, DimZ, NumOfb0 ]
80 // - spacing: old one, 1.0 time
82 typename OutputImageType::SpacingType spacing;
85 typename OutputImageType::PointType origin;
88 OutputImageRegionType outputRegion;
90 // the spacing and origin corresponds to the respective values in the input image (3D)
91 // the same for the region
92 for (unsigned int i=0; i< 3; i++)
94 spacing[i] = (this->GetInput()->GetSpacing())[i];
95 origin[i] = (this->GetInput()->GetOrigin())[i];
96 outputRegion.SetSize(i, this->GetInput()->GetLargestPossibleRegion().GetSize()[i]);
97 outputRegion.SetIndex(i, this->GetInput()->GetLargestPossibleRegion().GetIndex()[i]);
100 // number of timesteps = number of b0 images
101 outputRegion.SetSize(3, indices.size());
102 outputRegion.SetIndex( 3, 0 );
104 // output image direction (4x4)
105 typename OutputImageType::DirectionType outputDirection;
106 //initialize to identity
107 outputDirection.SetIdentity();
108 // get the input image direction ( 3x3 matrix )
109 typename InputImageType::DirectionType inputDirection = this->GetInput()->GetDirection();
111 for( unsigned int i=0; i< 3; i++)
113 outputDirection(0,i) = inputDirection(0,i);
114 outputDirection(1,i) = inputDirection(1,i);
115 outputDirection(2,i) = inputDirection(2,i);
118 outputImage->SetSpacing( spacing );
119 outputImage->SetOrigin( origin );
120 outputImage->SetDirection( outputDirection );
121 outputImage->SetRegions( outputRegion );
122 outputImage->Allocate();
125 itk::ImageRegionConstIterator<InputImageType> inputIt( this->GetInput(), this->GetInput()->GetLargestPossibleRegion() );
127 // we want to iterate separately over each 3D volume of the output image
128 // so reset the regions last dimension
129 outputRegion.SetSize(3,1);
130 unsigned int currentTimeStep = 0;
133 // extract b0 and insert them as a new time step
134 for(std::vector<int>::iterator indexIt = indices.begin();
135 indexIt != indices.end();
139 outputRegion.SetIndex(3, currentTimeStep);
140 itk::ImageRegionIterator< OutputImageType> outputIt( outputImage.GetPointer(), outputRegion );
142 // iterate over the current b0 image and store it to corresponding place
143 outputIt.GoToBegin();
145 while( !outputIt.IsAtEnd() && !inputIt.IsAtEnd() )
148 typename InputImageType::PixelType vec = inputIt.Get();
150 outputIt.Set( vec[*indexIt]);
155 // increase time step
162 #endif // ifndef __itkB0ImageExtractionToSeparateImageFilter_txx