ANIMA  4.0
animaComputeMahalanobisImagesFilter.hxx
Go to the documentation of this file.
1 #pragma once
2 
4 
5 namespace anima
6 {
7 
8 template <typename TInputImage, typename TMaskImage, typename TOutput>
10 {
11  this->SetNthInput(0, const_cast<TInputImage*>(image));
12 }
13 
14 template <typename TInputImage, typename TMaskImage, typename TOutput>
16 {
17  this->SetNthInput(1, const_cast<TInputImage*>(image));
18 }
19 
20 template <typename TInputImage, typename TMaskImage, typename TOutput>
22 {
23  this->SetNthInput(2, const_cast<TInputImage*>(image));
24 }
25 
26 template <typename TInputImage, typename TMaskImage, typename TOutput>
28 {
29  this->SetNthInput(3, const_cast<TMaskImage*>(mask));
30 }
31 
32 template <typename TInputImage, typename TMaskImage, typename TOutput>
34 {
35  return static_cast< const TInputImage * >
36  ( this->itk::ProcessObject::GetInput(0) );
37 }
38 
39 template <typename TInputImage, typename TMaskImage, typename TOutput>
41 {
42  return static_cast< const TInputImage * >
43  ( this->itk::ProcessObject::GetInput(1) );
44 }
45 
46 template <typename TInputImage, typename TMaskImage, typename TOutput>
48 {
49  return static_cast< const TInputImage * >
50  ( this->itk::ProcessObject::GetInput(2) );
51 }
52 
53 template <typename TInputImage, typename TMaskImage, typename TOutput>
55 {
56  return static_cast< const TMaskImage * >
57  ( this->itk::ProcessObject::GetInput(3) );
58 }
59 
60 
61 template <typename TInputImage, typename TMaskImage, typename TOutput>
63 {
64  itk::DataObject::Pointer output;
65 
66  switch ( idx )
67  {
68  case 0:
69  output = ( TOutput::New() ).GetPointer();
70  break;
71  case 1:
72  output = ( TOutput::New() ).GetPointer();
73  break;
74  case 2:
75  output = ( TOutput::New() ).GetPointer();
76  break;
77  case 3:
78  output = ( TOutput::New() ).GetPointer();
79  break;
80  case 4:
81  output = ( TOutput::New() ).GetPointer();
82  break;
83  default:
84  std::cerr << "No output " << idx << std::endl;
85  output = NULL;
86  break;
87  }
88  return output.GetPointer();
89 }
90 
91 template <typename TInputImage, typename TMaskImage, typename TOutput>
93 {
94  return dynamic_cast< OutputImageType * >( this->itk::ProcessObject::GetOutput(0) );
95 }
96 template <typename TInputImage, typename TMaskImage, typename TOutput>
98 {
99  return dynamic_cast< OutputImageType * >( this->itk::ProcessObject::GetOutput(1) );
100 }
101 template <typename TInputImage, typename TMaskImage, typename TOutput>
103 {
104  return dynamic_cast< OutputImageType * >( this->itk::ProcessObject::GetOutput(2) );
105 }
106 template <typename TInputImage, typename TMaskImage, typename TOutput>
108 {
109  return dynamic_cast< OutputImageType * >( this->itk::ProcessObject::GetOutput(3) );
110 }
111 template <typename TInputImage, typename TMaskImage, typename TOutput>
113 {
114  return dynamic_cast< OutputImageType * >( this->itk::ProcessObject::GetOutput(4) );
115 }
116 
117 template <typename TInputImage, typename TMaskImage, typename TOutput>
118 void
120 {
121  if( m_OutputMahaCSFFilename != "" )
122  {
123  std::cout << "Writing mahalanobis CSF image to: " << m_OutputMahaCSFFilename << std::endl;
124  anima::writeImage<TOutput>(m_OutputMahaCSFFilename, this->GetOutputMahaCSF());
125  }
126  if( m_OutputMahaGMFilename != "" )
127  {
128  std::cout << "Writing mahalanobis GM image to: " << m_OutputMahaGMFilename << std::endl;
129  anima::writeImage<TOutput>(m_OutputMahaGMFilename, this->GetOutputMahaGM());
130  }
131  if( m_OutputMahaWMFilename != "" )
132  {
133  std::cout << "Writing mahalanobis WM image to: " << m_OutputMahaWMFilename << std::endl;
134  anima::writeImage<TOutput>(m_OutputMahaWMFilename, this->GetOutputMahaWM());
135  }
136  if( m_OutputMahaMinimumFilename != "" )
137  {
138  std::cout << "Writing minimum mahalanobis image to: " << m_OutputMahaMinimumFilename << std::endl;
139  anima::writeImage<TOutput>(m_OutputMahaMinimumFilename, this->GetOutputMahaMinimum());
140  }
141  if( m_OutputMahaMaximumFilename != "" )
142  {
143  std::cout << "Writing maximum mahalanobis image to: " << m_OutputMahaMaximumFilename << std::endl;
144  anima::writeImage<TOutput>(m_OutputMahaMaximumFilename, this->GetOutputMahaMaximum());
145  }
146 }
147 
148 template <typename TInputImage, typename TMaskImage, typename TOutput>
149 void
151 ::BeforeThreadedGenerateData()
152 {
153  if(m_GaussianModel.size() != 3)
154  {
155  std::cout << "Error in mahalanobis filter: number of images do not correspond to solution size exiting..."<< std::endl;
156  exit(-1);
157  }
158 
159  m_InputImage_1_UC = ImageTypeUC::New();
160  m_InputImage_2_UC = ImageTypeUC::New();
161  m_InputImage_3_UC = ImageTypeUC::New();
162 
163  // Rescale images
164  double desiredMinimum=0,desiredMaximum=255;
165 
166  typename RescaleFilterType::Pointer rescaleFilter1 = RescaleFilterType::New();
167  rescaleFilter1->SetInput( this->GetInputImage1() );
168  rescaleFilter1->SetOutputMinimum( desiredMinimum );
169  rescaleFilter1->SetOutputMaximum( desiredMaximum );
170  rescaleFilter1->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
171  rescaleFilter1->SetCoordinateTolerance(m_Tol);
172  rescaleFilter1->SetDirectionTolerance(m_Tol);
173  rescaleFilter1->UpdateLargestPossibleRegion();
174  m_InputImage_1_UC = rescaleFilter1->GetOutput();
175 
176  typename RescaleFilterType::Pointer rescaleFilter2 = RescaleFilterType::New();
177  rescaleFilter2->SetInput( this->GetInputImage2() );
178  rescaleFilter2->SetOutputMinimum( desiredMinimum );
179  rescaleFilter2->SetOutputMaximum( desiredMaximum );
180  rescaleFilter2->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
181  rescaleFilter2->SetCoordinateTolerance(m_Tol);
182  rescaleFilter2->SetDirectionTolerance(m_Tol);
183  rescaleFilter2->UpdateLargestPossibleRegion();
184  m_InputImage_2_UC = rescaleFilter2->GetOutput();
185 
186  typename RescaleFilterType::Pointer rescaleFilter3 = RescaleFilterType::New();
187  rescaleFilter3->SetInput( this->GetInputImage3() );
188  rescaleFilter3->SetOutputMinimum( desiredMinimum );
189  rescaleFilter3->SetOutputMaximum( desiredMaximum );
190  rescaleFilter3->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
191  rescaleFilter3->SetCoordinateTolerance(m_Tol);
192  rescaleFilter3->SetDirectionTolerance(m_Tol);
193  rescaleFilter3->UpdateLargestPossibleRegion();
194  m_InputImage_3_UC = rescaleFilter3->GetOutput();
195 }
196 
197 
198 template <typename TInputImage, typename TMaskImage, typename TOutput>
199 void
201 {
202  const unsigned int Dimension = 3;
203  typedef double PixelComponentType;
204  typedef itk::RGBPixel<PixelComponentType> PixelType;
205  typedef itk::Image< PixelType, Dimension > ImageType;
206  typedef itk::MahalanobisDistanceThresholdImageFunction< ImageType > FunctionType;
207  typedef FunctionType::IndexType FunctionIndexType;
208 
209  OutputIteratorType mahaCSFIt(this->GetOutputMahaCSF(), outputRegionForThread );
210  OutputIteratorType mahaGMIt(this->GetOutputMahaGM(), outputRegionForThread );
211  OutputIteratorType mahaWMIt(this->GetOutputMahaWM(), outputRegionForThread );
212 
213  OutputIteratorType mahaMaxiIt(this->GetOutputMahaMaximum(), outputRegionForThread );
214  OutputIteratorType mahaMiniIt(this->GetOutputMahaMinimum(), outputRegionForThread );
215 
216  ImageType::Pointer image = ImageType::New();
217  image->SetRegions( this->GetMask()->GetLargestPossibleRegion() );
218  image->CopyInformation( this->GetMask() );
219  image->Allocate();
220  itk::ImageRegionIterator<ImageType> itRGB (image, outputRegionForThread);
221 
222  ImageConstIteratorTypeUC It1(m_InputImage_1_UC, outputRegionForThread );
223  ImageConstIteratorTypeUC It2(m_InputImage_2_UC, outputRegionForThread );
224  ImageConstIteratorTypeUC It3(m_InputImage_3_UC, outputRegionForThread );
225 
226  while(!itRGB.IsAtEnd())
227  {
228  PixelType pix;
229  pix[0] = It1.Get();
230  pix[1] = It2.Get();
231  pix[2] = It3.Get();
232  itRGB.Set(pix);
233  ++It1;
234  ++It2;
235  ++It3;
236  ++itRGB;
237  }
238 
239  FunctionType::Pointer function = FunctionType::New();
240  FunctionType::CovarianceMatrixType Covariance( m_NbModalities, m_NbModalities );
241  Covariance.fill( 0.0 );
242 
243  FunctionType::MeanVectorType Mean( m_NbModalities );
244  function->SetInputImage( image );
245 
246  MaskConstIteratorType MaskImageIt ( this->GetMask(), outputRegionForThread );
247 
248  itk::Statistics::ChiSquareDistribution::Pointer chi = itk::Statistics::ChiSquareDistribution::New();
249  itk::SizeValueType degree = (m_GaussianModel[0])->GetMean().Size();
250  chi->SetDegreesOfFreedom(degree);
251 
252  while(!MaskImageIt.IsAtEnd())
253  {
254  mahaCSFIt.Set(0);
255  mahaGMIt.Set(0);
256  mahaWMIt.Set(0);
257  mahaMaxiIt.Set(0);
258  mahaMiniIt.Set(0);
259 
260  if(MaskImageIt.Get()!=0)
261  {
262  FunctionIndexType ind;
263  ind[0] = MaskImageIt.GetIndex()[0];
264  ind[1] = MaskImageIt.GetIndex()[1];
265  ind[2] = MaskImageIt.GetIndex()[2];
266 
267  for(unsigned int i = 0; i < m_NbTissus; i++)
268  {
269  GaussianFunctionType::MeanVectorType mean = (m_GaussianModel[i])->GetMean();
270  Mean[0] = mean[0];
271  Mean[1] = mean[1];
272  Mean[2] = mean[2];
273 
274  GaussianFunctionType::CovarianceMatrixType cov = (m_GaussianModel[i])->GetCovariance();
275  Covariance[0][0] = cov[0][0]; Covariance[1][0] = cov[1][0]; Covariance[2][0] = cov[2][0];
276  Covariance[0][1] = cov[0][1]; Covariance[1][1] = cov[1][1]; Covariance[2][1] = cov[2][1];
277  Covariance[0][2] = cov[0][2]; Covariance[1][2] = cov[1][2]; Covariance[2][2] = cov[2][2];
278 
279  function->SetCovariance( Covariance );
280  function->SetMean( Mean );
281  double distanceD = function->EvaluateDistanceAtIndex( ind );
282  double chiRes = chi->EvaluateCDF(distanceD);
283  double distance = 1-chiRes;
284 
285  if(i==0)
286  {
287  mahaCSFIt.Set(static_cast<OutputPixelType>(distance));
288  }
289  if(i==1)
290  {
291  mahaGMIt.Set(static_cast<OutputPixelType>(distance));
292  }
293  if(i==2)
294  {
295  mahaWMIt.Set(static_cast<OutputPixelType>(distance));
296  }
297 
298 
299  if(mahaMaxiIt.Get() < mahaCSFIt.Get())
300  {
301  mahaMaxiIt.Set(mahaCSFIt.Get());
302  }
303  if(mahaMaxiIt.Get() < mahaGMIt.Get())
304  {
305  mahaMaxiIt.Set(mahaGMIt.Get());
306  }
307  if(mahaMaxiIt.Get() < mahaWMIt.Get())
308  {
309  mahaMaxiIt.Set(mahaWMIt.Get());
310  }
311  mahaMiniIt.Set(1.0 - mahaMaxiIt.Get());
312  }
313  }
314 
315  ++mahaCSFIt;
316  ++mahaGMIt;
317  ++mahaWMIt;
318  ++mahaMaxiIt;
319  ++mahaMiniIt;
320  ++MaskImageIt;
321  }
322 }
323 
324 } //end of namespace anima
itk::ImageRegionIterator< OutputImageType > OutputIteratorType
Compute the mahalanobis images from the NABT model.
itk::ImageRegionConstIterator< MaskImageType > MaskConstIteratorType
itk::ImageRegionConstIterator< ImageTypeUC > ImageConstIteratorTypeUC
Superclass::OutputImageRegionType OutputImageRegionType