Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkB0ImageExtractionImageFilter.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 
17 /*===================================================================
18 
19 This file is based heavily on a corresponding ITK filter.
20 
21 ===================================================================*/
22 #ifndef __itkB0ImageExtractionImageFilter_txx
23 #define __itkB0ImageExtractionImageFilter_txx
24 
25 #include "itkB0ImageExtractionImageFilter.h"
26 
27 #include "itkImageRegionIterator.h"
28 
29 namespace itk {
30 
31  template< class TInputImagePixelType,
32  class TOutputImagePixelType >
33  B0ImageExtractionImageFilter< TInputImagePixelType,
34  TOutputImagePixelType >
35  ::B0ImageExtractionImageFilter()
36  {
37  // At least 1 inputs is necessary for a vector image.
38  // For images added one at a time we need at least six
39  this->SetNumberOfRequiredInputs( 1 );
40  }
41 
42  template< class TInputImagePixelType,
43  class TOutputImagePixelType >
44  void B0ImageExtractionImageFilter< TInputImagePixelType,
45  TOutputImagePixelType >
46  ::GenerateData()
47  {
48 
49  typename GradientDirectionContainerType::Iterator begin = m_Directions->Begin();
50  typename GradientDirectionContainerType::Iterator end = m_Directions->End();
51 
52  // Find the index of the b0 image
53  std::vector<int> indices;
54  int index = 0;
55  while(begin!=end)
56  {
57  GradientDirectionType grad = begin->Value();
58 
59  if(grad[0] == 0 && grad[1] == 0 && grad[2] == 0)
60  {
61  indices.push_back(index);
62  }
63  ++index;
64  ++begin;
65  }
66 
67  typedef itk::Image<float,3> TempImageType;
68  TempImageType::Pointer tmp = TempImageType::New();
69  typename TempImageType::RegionType region = this->GetInput()->GetLargestPossibleRegion();
70 
71  tmp->SetSpacing(this->GetInput()->GetSpacing());
72  tmp->SetOrigin(this->GetInput()->GetOrigin());
73  tmp->SetDirection(this->GetInput()->GetDirection());
74  tmp->SetRegions(region);
75  tmp->Allocate();
76 
77  itk::ImageRegionIterator<TempImageType> it(tmp.GetPointer(), tmp->GetLargestPossibleRegion() );
78  itk::ImageRegionConstIterator<InputImageType> vectorIt(this->GetInput(), this->GetInput()->GetLargestPossibleRegion() );
79 
80  it.GoToBegin();
81  while(!it.IsAtEnd())
82  {
83  it.Set(0);
84  ++it;
85  }
86 
87  //Sum all images that have zero diffusion weighting (indices stored in vector index)
88  for(std::vector<int>::iterator indexIt = indices.begin();
89  indexIt != indices.end();
90  indexIt++)
91  {
92  it.GoToBegin();
93  vectorIt.GoToBegin();
94 
95  while(!it.IsAtEnd() && !vectorIt.IsAtEnd())
96  {
97  typename InputImageType::PixelType vec = vectorIt.Get();
98  it.Set((1.0 * it.Get()) + (1.0 * vec[*indexIt]) / (1.0 * indices.size()));
99  ++it;
100  ++vectorIt;
101  }
102  }
103 
104  typename OutputImageType::Pointer b0Image =
105  static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
106  typename OutputImageType::RegionType outregion = this->GetInput()->GetLargestPossibleRegion();
107  b0Image->SetSpacing(this->GetInput()->GetSpacing());
108  b0Image->SetOrigin(this->GetInput()->GetOrigin());
109  b0Image->SetDirection(this->GetInput()->GetDirection());
110  b0Image->SetRegions(outregion);
111  b0Image->Allocate();
112  itk::ImageRegionIterator<TempImageType> itIn(tmp, tmp->GetLargestPossibleRegion() );
113  itk::ImageRegionIterator<OutputImageType> itOut(b0Image, b0Image->GetLargestPossibleRegion() );
114  itIn.GoToBegin();
115  itOut.GoToBegin();
116  while(!itIn.IsAtEnd())
117  {
118  itOut.Set(itIn.Get());
119  ++itIn;
120  ++itOut;
121  }
122  }
123 
124 }
125 
126 #endif // __itkB0ImageExtractionImageFilter_txx