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 Program: Tensor ToolKit - TTK
19 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $
21 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
22 Version: $Revision: 68 $
24 Copyright (c) INRIA 2010. All rights reserved.
25 See LICENSE.txt for details.
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.
31 =========================================================================*/
32 #ifndef _itk_ResidualImageFilter_txx_
33 #define _itk_ResidualImageFilter_txx_
36 #include "itkResidualImageFilter.h"
37 #include <mitkCommon.h>
46 template <class TInputScalarType, class TOutputScalarType>
48 ResidualImageFilter<TInputScalarType, TOutputScalarType>
51 typename InputImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
52 typename InputImageType::SizeType size2 = m_SecondDiffusionImage->GetLargestPossibleRegion().GetSize();
56 MITK_ERROR << "Sizes do not match";
60 // Initialize output image
61 typename OutputImageType::Pointer outputImage =
62 static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
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);
75 std::vector< std::vector<double> > residuals;
77 // per slice, per volume
78 std::vector< std::vector <std::vector<double> > > residualsPerSlice;
80 // Detrmine number of B0 images
82 for(unsigned int i=0; i<m_Gradients->Size(); i++)
84 GradientDirectionType grad = m_Gradients->ElementAt(i);
86 if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001)
92 residuals.resize(this->GetInput()->GetVectorLength()-numberB0);
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++)
97 std::vector< std::vector<double> > sliceResiduals; // residuals per volume for this slice
98 sliceResiduals.resize(this->GetInput()->GetVectorLength()-numberB0);
100 for(unsigned int y=0; y<size[1]; y++)
102 for(unsigned int x=0; x<size[0]; x++)
105 // Check if b0 exceeds threshold
112 typename InputImageType::PixelType p1 = this->GetInput()->GetPixel(ix);
113 typename InputImageType::PixelType p2 = m_SecondDiffusionImage->GetPixel(ix);
115 if(p1.GetSize() != p2.GetSize())
117 MITK_ERROR << "Vector sizes do not match";
122 if(p1.GetElement(m_B0Index) <= m_B0Threshold)
130 int shift = 0; // correction for the skipped B0 images
132 for(unsigned int i = 0; i<p1.GetSize(); i++)
137 GradientDirectionType grad = m_Gradients->ElementAt(i);
138 if(!(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001))
140 double val1 = (double)p1.GetElement(i);
141 double val2 = (double)p2.GetElement(i);
143 res += abs(val1-val2);
145 residuals[i-shift].push_back(val1-val2);
146 sliceResiduals[i-shift].push_back(val1-val2);
159 res = res/p1.GetSize();
161 outputImage->SetPixel(ix, res);
166 residualsPerSlice.push_back(sliceResiduals);
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
177 double q1,q3, median;
178 std::vector< std::vector<double> >::iterator it = residuals.begin();
179 while(it != residuals.end())
181 std::vector<double> res = *it;
184 std::sort(res.begin(), res.end());
186 q1 = res[0.25*res.size()];
188 q3 = res[0.75*res.size()];
193 median = res[0.5*res.size()];
195 double outlierThreshold = median + 1.5*iqr;
196 double numberOfOutliers = 0.0;
198 std::vector<double>::iterator resIt = res.begin();
200 while(resIt != res.end())
203 if(f>outlierThreshold)
211 double percOfOutliers = 100 * numberOfOutliers / res.size();
212 m_PercentagesOfOutliers.push_back(percOfOutliers);
215 m_Means.push_back(mean);
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();
225 while(sliceIt != residualsPerSlice.end())
227 std::vector< std::vector<double> > currentSlice = *sliceIt;
228 std::vector<double> percentages;
230 std::vector< std::vector<double> >::iterator volIt = currentSlice.begin();
231 while(volIt != currentSlice.end())
235 std::vector<double> currentVolume = *volIt;
238 std::sort(currentVolume.begin(), currentVolume.end());
242 q1 = currentVolume[0.25*currentVolume.size()];
243 q3 = currentVolume[0.75*currentVolume.size()];
244 median = currentVolume[0.5*currentVolume.size()];
246 double outlierThreshold = median + 1.5*iqr;
247 double numberOfOutliers = 0.0;
249 std::vector<double>::iterator resIt = currentVolume.begin();
251 while(resIt != currentVolume.end())
254 if(f>outlierThreshold)
262 double percOfOutliers = 100 * numberOfOutliers / currentVolume.size();
263 percentages.push_back(percOfOutliers);
268 m_OutliersPerSlice.push_back(percentages);
285 } // end of namespace