Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkB0ImageExtractionToSeparateImageFilter.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 #ifndef __itkB0ImageExtractionToSeparateImageFilter_txx
17 #define __itkB0ImageExtractionToSeparateImageFilter_txx
18 
19 #include "itkB0ImageExtractionToSeparateImageFilter.h"
20 
21 template< class TInputImagePixelType,
22  class TOutputImagePixelType>
23 itk::B0ImageExtractionToSeparateImageFilter< TInputImagePixelType, TOutputImagePixelType >
24 ::B0ImageExtractionToSeparateImageFilter()
25 {
26  this->SetNumberOfRequiredInputs( 1 );
27 }
28 
29 template< class TInputImagePixelType,
30  class TOutputImagePixelType>
31 void itk::B0ImageExtractionToSeparateImageFilter<
32 TInputImagePixelType, TOutputImagePixelType >::CopyInformation( const DataObject *)
33 {
34 
35 }
36 
37 template< class TInputImagePixelType,
38  class TOutputImagePixelType>
39 void itk::B0ImageExtractionToSeparateImageFilter<
40 TInputImagePixelType, TOutputImagePixelType >::GenerateOutputInformation()
41 {
42 
43 }
44 
45 template< class TInputImagePixelType,
46  class TOutputImagePixelType>
47 void itk::B0ImageExtractionToSeparateImageFilter<
48 TInputImagePixelType, TOutputImagePixelType >::GenerateData()
49 {
50 
51  GradContainerIteratorType begin = m_Directions->Begin();
52  GradContainerIteratorType end = m_Directions->End();
53 
54  // Find the indices of the b0 images in the diffusion image
55  std::vector<int> indices;
56  int index = 0;
57  while(begin!=end)
58  {
59  GradientDirectionType grad = begin->Value();
60 
61  if(grad[0] == 0 && grad[1] == 0 && grad[2] == 0)
62  {
63  indices.push_back(index);
64  }
65  ++index;
66  ++begin;
67  }
68 
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();
73 
74  // get the input region
75  typename Superclass::InputImageType::RegionType inputRegion =
76  this->GetInput()->GetLargestPossibleRegion();
77 
78  // allocate image with
79  // - dimension: [DimX, DimY, DimZ, NumOfb0 ]
80  // - spacing: old one, 1.0 time
81 
82  typename OutputImageType::SpacingType spacing;
83  spacing.Fill(1);
84 
85  typename OutputImageType::PointType origin;
86  origin.Fill(0);
87 
88  OutputImageRegionType outputRegion;
89 
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++)
93  {
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]);
98  }
99 
100  // number of timesteps = number of b0 images
101  outputRegion.SetSize(3, indices.size());
102  outputRegion.SetIndex( 3, 0 );
103 
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();
110 
111  for( unsigned int i=0; i< 3; i++)
112  {
113  outputDirection(0,i) = inputDirection(0,i);
114  outputDirection(1,i) = inputDirection(1,i);
115  outputDirection(2,i) = inputDirection(2,i);
116  }
117 
118  outputImage->SetSpacing( spacing );
119  outputImage->SetOrigin( origin );
120  outputImage->SetDirection( outputDirection );
121  outputImage->SetRegions( outputRegion );
122  outputImage->Allocate();
123 
124  // input iterator
125  itk::ImageRegionConstIterator<InputImageType> inputIt( this->GetInput(), this->GetInput()->GetLargestPossibleRegion() );
126 
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;
131 
132 
133  // extract b0 and insert them as a new time step
134  for(std::vector<int>::iterator indexIt = indices.begin();
135  indexIt != indices.end();
136  indexIt++)
137  {
138  // set the time step
139  outputRegion.SetIndex(3, currentTimeStep);
140  itk::ImageRegionIterator< OutputImageType> outputIt( outputImage.GetPointer(), outputRegion );
141 
142  // iterate over the current b0 image and store it to corresponding place
143  outputIt.GoToBegin();
144  inputIt.GoToBegin();
145  while( !outputIt.IsAtEnd() && !inputIt.IsAtEnd() )
146  {
147  // the input vector
148  typename InputImageType::PixelType vec = inputIt.Get();
149 
150  outputIt.Set( vec[*indexIt]);
151  ++outputIt;
152  ++inputIt;
153  }
154 
155  // increase time step
156  currentTimeStep++;
157  }
158 
159 
160 }
161 
162 #endif // ifndef __itkB0ImageExtractionToSeparateImageFilter_txx