Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkMergeDiffusionImagesFilter.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  2
18  3 Program: Tensor ToolKit - TTK
19  4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkMergeDiffusionImagesFilter.txx $
20  5 Language: C++
21  6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
22  7 Version: $Revision: 68 $
23  8
24  9 Copyright (c) INRIA 2010. All rights reserved.
25  10 See LICENSE.txt for details.
26  11
27  12 This software is distributed WITHOUT ANY WARRANTY; without even
28  13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
29  14 PURPOSE. See the above copyright notices for more information.
30  15
31  16 =========================================================================*/
32 #ifndef _itk_MergeDiffusionImagesFilter_txx_
33 #define _itk_MergeDiffusionImagesFilter_txx_
34 #endif
35 
36 #include "itkMergeDiffusionImagesFilter.h"
37 #include "itkTensorToL2NormImageFilter.h"
38 #include "itkRescaleIntensityImageFilter.h"
39 #include <itkImageRegionIterator.h>
40 #include <itkImageRegionConstIterator.h>
41 #include <boost/progress.hpp>
42 
43 namespace itk
44 {
45 
46 template <class TScalarType>
47 MergeDiffusionImagesFilter<TScalarType>::MergeDiffusionImagesFilter()
48 {
49 
50 }
51 
52 template <class TScalarType>
53 MergeDiffusionImagesFilter<TScalarType>::~MergeDiffusionImagesFilter()
54 {
55 
56 }
57 
58 template <class TScalarType>
59 void MergeDiffusionImagesFilter<TScalarType>
60 ::SetImageVolumes(DwiImageContainerType cont)
61 {
62  m_ImageVolumes=cont;
63 }
64 
65 template <class TScalarType>
66 void MergeDiffusionImagesFilter<TScalarType>
67 ::SetGradientLists(GradientListContainerType cont)
68 {
69  m_GradientLists=cont;
70 }
71 
72 template <class TScalarType>
73 void MergeDiffusionImagesFilter<TScalarType>
74 ::SetBValues(std::vector< double > bvals)
75 {
76  m_BValues=bvals;
77 }
78 
79 template <class TScalarType>
80 MergeDiffusionImagesFilter<TScalarType>::GradientListType::Pointer MergeDiffusionImagesFilter<TScalarType>
81 ::GetOutputGradients()
82 {
83  return m_OutputGradients;
84 }
85 
86 template <class TScalarType>
87 double MergeDiffusionImagesFilter<TScalarType>
88 ::GetB_Value()
89 {
90  return m_BValue;
91 }
92 
93 template <class TScalarType>
94 void
95 MergeDiffusionImagesFilter<TScalarType>
96 ::GenerateData ()
97 {
98 
99  if( m_ImageVolumes.size()<2 )
100  throw itk::ExceptionObject (__FILE__,__LINE__,"Error: cannot combine less than two DWIs.");
101 
102  if( m_GradientLists.size()!=m_ImageVolumes.size() || m_ImageVolumes.size()!=m_BValues.size() || m_BValues.size()!=m_GradientLists.size() )
103  throw itk::ExceptionObject (__FILE__,__LINE__,"Error: need same number of b-values, image volumes and gradient containers.");
104 
105  typename DwiImageType::Pointer img = m_ImageVolumes.at(0);
106 
107  m_NumGradients = 0;
108  for (unsigned int i=0; i<m_GradientLists.size(); i++)
109  {
110  m_NumGradients += m_GradientLists.at(i)->Size();
111  typename DwiImageType::Pointer tmp = m_ImageVolumes.at(i);
112  if ( img->GetLargestPossibleRegion()!=tmp->GetLargestPossibleRegion() )
113  throw itk::ExceptionObject (__FILE__,__LINE__,"Error: images are not of same size.");
114  }
115 
116  m_BValue = m_BValues.at(0);
117  m_OutputGradients = GradientListType::New();
118 
119  typename DwiImageType::Pointer outImage = DwiImageType::New();
120  outImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
121  outImage->SetOrigin( img->GetOrigin() ); // Set the image origin
122  outImage->SetDirection( img->GetDirection() ); // Set the image direction
123  outImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
124  outImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
125  outImage->SetRequestedRegion( img->GetLargestPossibleRegion() );
126  outImage->SetVectorLength(m_NumGradients);
127  outImage->Allocate();
128 
129  this->SetNumberOfRequiredOutputs(1);
130  this->SetNthOutput (0, outImage);
131 
132  typedef ImageRegionIterator<DwiImageType> IteratorOutputType;
133  IteratorOutputType itOut (this->GetOutput(), this->GetOutput()->GetLargestPossibleRegion());
134 
135  MITK_INFO << "MergeDiffusionImagesFilter: merging images";
136  GradientType zeroG; zeroG.fill(0.0);
137  boost::progress_display disp(this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels());
138  while(!itOut.IsAtEnd())
139  {
140  ++disp;
141  DwiPixelType out;
142  out.SetSize(m_NumGradients);
143  out.Fill(0);
144 
145  int c=0;
146  for (unsigned int i=0; i<m_GradientLists.size(); i++)
147  {
148  GradientListType::Pointer gradients = m_GradientLists.at(i);
149  typename DwiImageType::Pointer img = m_ImageVolumes.at(i);
150 
151  for (unsigned int j=0; j<gradients->Size(); j++)
152  {
153  GradientType g = gradients->GetElement(j);
154  double mag = g.two_norm();
155 
156  if (mag>0.0001)
157  {
158  double frac = m_BValues.at(i)*mag*mag/m_BValue;
159  g.normalize();
160  g *= sqrt(frac);
161  }
162  else
163  g = zeroG;
164 
165  m_OutputGradients->InsertElement(c, g);
166  out[c] = static_cast<TScalarType>(img->GetPixel(itOut.GetIndex())[j]);
167  c++;
168  }
169  }
170 
171  itOut.Set(out);
172  ++itOut;
173  }
174 }
175 
176 } // end of namespace