Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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