1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
16 /*=========================================================================
18 3 Program: Tensor ToolKit - TTK
19 4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkMergeDiffusionImagesFilter.txx $
21 6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
22 7 Version: $Revision: 68 $
24 9 Copyright (c) INRIA 2010. All rights reserved.
25 10 See LICENSE.txt for details.
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.
31 16 =========================================================================*/
32 #ifndef _itk_MergeDiffusionImagesFilter_txx_
33 #define _itk_MergeDiffusionImagesFilter_txx_
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>
46 template <class TScalarType>
47 MergeDiffusionImagesFilter<TScalarType>::MergeDiffusionImagesFilter()
52 template <class TScalarType>
53 MergeDiffusionImagesFilter<TScalarType>::~MergeDiffusionImagesFilter()
58 template <class TScalarType>
59 void MergeDiffusionImagesFilter<TScalarType>
60 ::SetImageVolumes(DwiImageContainerType cont)
65 template <class TScalarType>
66 void MergeDiffusionImagesFilter<TScalarType>
67 ::SetGradientLists(GradientListContainerType cont)
72 template <class TScalarType>
73 void MergeDiffusionImagesFilter<TScalarType>
74 ::SetBValues(std::vector< double > bvals)
79 template <class TScalarType>
80 MergeDiffusionImagesFilter<TScalarType>::GradientListType::Pointer MergeDiffusionImagesFilter<TScalarType>
81 ::GetOutputGradients()
83 return m_OutputGradients;
86 template <class TScalarType>
87 double MergeDiffusionImagesFilter<TScalarType>
93 template <class TScalarType>
95 MergeDiffusionImagesFilter<TScalarType>
99 if( m_ImageVolumes.size()<2 )
100 throw itk::ExceptionObject (__FILE__,__LINE__,"Error: cannot combine less than two DWIs.");
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.");
105 typename DwiImageType::Pointer img = m_ImageVolumes.at(0);
108 for (unsigned int i=0; i<m_GradientLists.size(); i++)
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.");
116 m_BValue = m_BValues.at(0);
117 m_OutputGradients = GradientListType::New();
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();
129 this->SetNumberOfRequiredOutputs(1);
130 this->SetNthOutput (0, outImage);
132 typedef ImageRegionIterator<DwiImageType> IteratorOutputType;
133 IteratorOutputType itOut (this->GetOutput(), this->GetOutput()->GetLargestPossibleRegion());
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())
142 out.SetSize(m_NumGradients);
146 for (unsigned int i=0; i<m_GradientLists.size(); i++)
148 GradientListType::Pointer gradients = m_GradientLists.at(i);
149 typename DwiImageType::Pointer img = m_ImageVolumes.at(i);
151 for (unsigned int j=0; j<gradients->Size(); j++)
153 GradientType g = gradients->GetElement(j);
154 double mag = g.two_norm();
158 double frac = m_BValues.at(i)*mag*mag/m_BValue;
165 m_OutputGradients->InsertElement(c, g);
166 out[c] = static_cast<TScalarType>(img->GetPixel(itOut.GetIndex())[j]);
176 } // end of namespace