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
itkResidualImageFilter.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 Program: Tensor ToolKit - TTK
19 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $
20 Language: C++
21 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
22 Version: $Revision: 68 $
23 
24 Copyright (c) INRIA 2010. All rights reserved.
25 See LICENSE.txt for details.
26 
27 This software is distributed WITHOUT ANY WARRANTY; without even
28 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
29 PURPOSE. See the above copyright notices for more information.
30 
31 =========================================================================*/
32 #ifndef _itk_ResidualImageFilter_txx_
33 #define _itk_ResidualImageFilter_txx_
34 #endif
35 
36 #include "itkResidualImageFilter.h"
37 #include <mitkCommon.h>
38 #include <cmath>
39 
40 namespace itk
41 {
42 
43 
44 
45 
46  template <class TInputScalarType, class TOutputScalarType>
47  void
48  ResidualImageFilter<TInputScalarType, TOutputScalarType>
49  ::GenerateData()
50  {
51  typename InputImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
52  typename InputImageType::SizeType size2 = m_SecondDiffusionImage->GetLargestPossibleRegion().GetSize();
53 
54  if(size != size2)
55  {
56  MITK_ERROR << "Sizes do not match";
57  return;
58  }
59 
60  // Initialize output image
61  typename OutputImageType::Pointer outputImage =
62  static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
63 
64  outputImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
65  outputImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
66  outputImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
67  outputImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
68  outputImage->Allocate();
69  outputImage->FillBuffer(0.0);
70 
71 
72 
73 
74 
75  std::vector< std::vector<double> > residuals;
76 
77  // per slice, per volume
78  std::vector< std::vector <std::vector<double> > > residualsPerSlice;
79 
80  // Detrmine number of B0 images
81  int numberB0=0;
82  for(unsigned int i=0; i<m_Gradients->Size(); i++)
83  {
84  GradientDirectionType grad = m_Gradients->ElementAt(i);
85 
86  if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001)
87  {
88  numberB0++;
89  }
90  }
91 
92  residuals.resize(this->GetInput()->GetVectorLength()-numberB0);
93 
94  // Calculate the standard residual image and for each volume put all residuals in a vector
95  for(unsigned int z=0; z<size[2]; z++)
96  {
97  std::vector< std::vector<double> > sliceResiduals; // residuals per volume for this slice
98  sliceResiduals.resize(this->GetInput()->GetVectorLength()-numberB0);
99 
100  for(unsigned int y=0; y<size[1]; y++)
101  {
102  for(unsigned int x=0; x<size[0]; x++)
103  {
104 
105  // Check if b0 exceeds threshold
106  itk::Index<3> ix;
107  ix[0] = x;
108  ix[1] = y;
109  ix[2] = z;
110 
111 
112  typename InputImageType::PixelType p1 = this->GetInput()->GetPixel(ix);
113  typename InputImageType::PixelType p2 = m_SecondDiffusionImage->GetPixel(ix);
114 
115  if(p1.GetSize() != p2.GetSize())
116  {
117  MITK_ERROR << "Vector sizes do not match";
118  return;
119  }
120 
121 
122  if(p1.GetElement(m_B0Index) <= m_B0Threshold)
123  {
124  continue;
125  }
126 
127 
128 
129  double res = 0;
130  int shift = 0; // correction for the skipped B0 images
131 
132  for(unsigned int i = 0; i<p1.GetSize(); i++)
133  {
134 
135 
136 
137  GradientDirectionType grad = m_Gradients->ElementAt(i);
138  if(!(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001))
139  {
140  double val1 = (double)p1.GetElement(i);
141  double val2 = (double)p2.GetElement(i);
142 
143  res += abs(val1-val2);
144 
145  residuals[i-shift].push_back(val1-val2);
146  sliceResiduals[i-shift].push_back(val1-val2);
147 
148  }
149  else
150  {
151  shift++;
152  }
153 
154 
155 
156  }
157 
158 
159  res = res/p1.GetSize();
160 
161  outputImage->SetPixel(ix, res);
162 
163  }
164  }
165 
166  residualsPerSlice.push_back(sliceResiduals);
167  }
168 
169 
170 
171 
172 
173  // for each dw volume: sort the the measured residuals (for each voxel) to enable determining Q1 and Q3; calculate means
174  // determine percentage of errors as described in QUALITY ASSESSMENT THROUGH ANALYSIS OF RESIDUALS OF DIFFUSION IMAGE FITTING
175  // Leemans et al 2008
176 
177  double q1,q3, median;
178  std::vector< std::vector<double> >::iterator it = residuals.begin();
179  while(it != residuals.end())
180  {
181  std::vector<double> res = *it;
182 
183  // sort
184  std::sort(res.begin(), res.end());
185 
186  q1 = res[0.25*res.size()];
187  m_Q1.push_back(q1);
188  q3 = res[0.75*res.size()];
189  m_Q3.push_back(q3);
190 
191 
192 
193  median = res[0.5*res.size()];
194  double iqr = q3-q1;
195  double outlierThreshold = median + 1.5*iqr;
196  double numberOfOutliers = 0.0;
197 
198  std::vector<double>::iterator resIt = res.begin();
199  double mean = 0;
200  while(resIt != res.end())
201  {
202  double f = *resIt;
203  if(f>outlierThreshold)
204  {
205  numberOfOutliers++;
206  }
207  mean += f;
208  ++resIt;
209  }
210 
211  double percOfOutliers = 100 * numberOfOutliers / res.size();
212  m_PercentagesOfOutliers.push_back(percOfOutliers);
213 
214  mean /= res.size();
215  m_Means.push_back(mean);
216 
217  ++it;
218  }
219 
220 
221 
222  // Calculate for each slice the number of outliers per volume(dw volume)
223  std::vector< std::vector <std::vector<double> > >::iterator sliceIt = residualsPerSlice.begin();
224 
225  while(sliceIt != residualsPerSlice.end())
226  {
227  std::vector< std::vector<double> > currentSlice = *sliceIt;
228  std::vector<double> percentages;
229 
230  std::vector< std::vector<double> >::iterator volIt = currentSlice.begin();
231  while(volIt != currentSlice.end())
232  {
233 
234 
235  std::vector<double> currentVolume = *volIt;
236 
237  //sort
238  std::sort(currentVolume.begin(), currentVolume.end());
239 
240 
241 
242  q1 = currentVolume[0.25*currentVolume.size()];
243  q3 = currentVolume[0.75*currentVolume.size()];
244  median = currentVolume[0.5*currentVolume.size()];
245  double iqr = q3-q1;
246  double outlierThreshold = median + 1.5*iqr;
247  double numberOfOutliers = 0.0;
248 
249  std::vector<double>::iterator resIt = currentVolume.begin();
250  double mean(0.0);
251  while(resIt != currentVolume.end())
252  {
253  double f = *resIt;
254  if(f>outlierThreshold)
255  {
256  numberOfOutliers++;
257  }
258  mean += f;
259  ++resIt;
260  }
261 
262  double percOfOutliers = 100 * numberOfOutliers / currentVolume.size();
263  percentages.push_back(percOfOutliers);
264 
265  ++volIt;
266  }
267 
268  m_OutliersPerSlice.push_back(percentages);
269 
270 
271 
272  ++sliceIt;
273  }
274 
275 
276 
277 
278  }
279 
280 
281 
282 
283 
284 
285 } // end of namespace