Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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