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
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