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