Orfeo Toolbox  3.16
itkBayesianClassifierImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkBayesianClassifierImageFilter.txx,v $
5  Language: C++
6  Date: $Date: 2010-01-30 21:05:25 $
7  Version: $Revision: 1.12 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  Portions of this code are covered under the VTK copyright.
13  See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details.
14 
15  This software is distributed WITHOUT ANY WARRANTY; without even
16  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
17  PURPOSE. See the above copyright notices for more information.
18 
19 =========================================================================*/
20 #ifndef __itkBayesianClassifierImageFilter_txx
21 #define __itkBayesianClassifierImageFilter_txx
22 
23 #include "itkBayesianClassifierImageFilter.h"
24 #include "itkImageRegionConstIterator.h"
25 
26 namespace itk
27 {
28 
32 template < class TInputVectorImage, class TLabelsType,
33  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
34 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
35  TPosteriorsPrecisionType, TPriorsPrecisionType >
36 ::BayesianClassifierImageFilter()
37 {
38  m_UserProvidedPriors = false;
39  m_UserProvidedSmoothingFilter = false;
40  this->SetNumberOfRequiredOutputs( 2 );
41  m_NumberOfSmoothingIterations = 0;
42  m_SmoothingFilter = NULL;
43  PosteriorsImagePointer p =
44  static_cast<PosteriorsImageType*>(this->MakeOutput(1).GetPointer());
45  this->SetNthOutput( 1 , p.GetPointer() );
46 }
47 
51 template < class TInputVectorImage, class TLabelsType,
52  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
53 void
54 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
55  TPosteriorsPrecisionType, TPriorsPrecisionType >
56 ::PrintSelf(std::ostream& os, Indent indent) const
57 {
58  Superclass::PrintSelf(os,indent);
59 
60  os << indent << "User provided priors = " << m_UserProvidedPriors << std::endl;
61  os << indent << "User provided smooting filter = " << m_UserProvidedSmoothingFilter << std::endl;
62  os << indent << "Smoothing filter pointer = " << m_SmoothingFilter.GetPointer() << std::endl;
63  os << indent << "Number of smoothing iterations = " << m_NumberOfSmoothingIterations << std::endl;
64 
65 }
66 
70 template < class TInputVectorImage, class TLabelsType,
71  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
72 void
73 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
74  TPosteriorsPrecisionType, TPriorsPrecisionType >
75 ::GenerateData()
76 {
77  // Setup input image
78  const InputImageType * membershipImage = this->GetInput();
79 
80  // Setup general parameters
81  const unsigned int numberOfClasses = membershipImage->GetVectorLength();
82 
83  if( numberOfClasses == 0 )
84  {
85  itkExceptionMacro("The number of components in the input Membership image is Zero !");
86  return;
87  }
88 
89  this->AllocateOutputs();
90 
91 
92  this->ComputeBayesRule();
93 
94  if( m_UserProvidedSmoothingFilter )
95  {
96  this->NormalizeAndSmoothPosteriors();
97  }
98 
99  this->ClassifyBasedOnPosteriors();
100 
101 }
102 
103 
104 template < class TInputVectorImage, class TLabelsType,
105  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
106 typename BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
107  TPosteriorsPrecisionType, TPriorsPrecisionType >
108 ::PosteriorsImageType *
109 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
110  TPosteriorsPrecisionType, TPriorsPrecisionType >
111 ::GetPosteriorImage()
112 {
113  PosteriorsImageType * ptr = dynamic_cast< PosteriorsImageType * >(
114  this->ProcessObject::GetOutput(1) );
115  return ptr;
116 }
117 
118 template < class TInputVectorImage, class TLabelsType,
119  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
120 typename BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
121  TPosteriorsPrecisionType, TPriorsPrecisionType >
122 ::DataObjectPointer
123 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
124  TPosteriorsPrecisionType, TPriorsPrecisionType >
125 ::MakeOutput(unsigned int idx)
126 {
127  if (idx == 1)
128  {
129  return static_cast<DataObject*>(PosteriorsImageType::New().GetPointer());
130  }
131  return Superclass::MakeOutput(idx);
132 }
133 
134 template < class TInputVectorImage, class TLabelsType,
135  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
136 void
137 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
138  TPosteriorsPrecisionType, TPriorsPrecisionType >
139 ::GenerateOutputInformation(void)
140 {
141  Superclass::GenerateOutputInformation();
142 
143  if ( !this->GetPosteriorImage() )
144  {
145  return;
146  }
147 
148  // the vector length is part of the output information that must be
149  // updated here
150  this->GetPosteriorImage()->SetVectorLength( this->GetInput()->GetVectorLength() );
151 }
152 
156 template < class TInputVectorImage, class TLabelsType,
157  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
158  void
159 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
160  TPosteriorsPrecisionType, TPriorsPrecisionType >
161 ::ComputeBayesRule()
162 {
163  itkDebugMacro( << "Computing Bayes Rule" );
164  const InputImageType * membershipImage = this->GetInput();
165 
166  ImageRegionType imageRegion = membershipImage->GetBufferedRegion();
167 
168  if( m_UserProvidedPriors )
169  {
170  const PriorsImageType * priorsImage =
171  dynamic_cast< const PriorsImageType * >( this->GetInput(1) );
172 
173  if( priorsImage == NULL )
174  {
175  itkExceptionMacro("Second input type does not correspond to expected Priors Image Type");
176  }
177 
178  PosteriorsImageType * posteriorsImage =
179  dynamic_cast< PosteriorsImageType * >( this->GetPosteriorImage() );
180 
181  if( posteriorsImage == NULL )
182  {
183  itkExceptionMacro("Second output type does not correspond to expected Posteriors Image Type");
184  }
185 
186 
187  InputImageIteratorType itrMembershipImage( membershipImage, imageRegion );
188  PriorsImageIteratorType itrPriorsImage( priorsImage, imageRegion );
189  PosteriorsImageIteratorType itrPosteriorsImage( posteriorsImage, imageRegion );
190 
191  itrMembershipImage.GoToBegin();
192  itrPriorsImage.GoToBegin();
193 
194  const unsigned int numberOfClasses = membershipImage->GetVectorLength();
195 
196  itkDebugMacro( << "Computing Bayes Rule nclasses in membershipImage: " << numberOfClasses );
197 
198  while( !itrMembershipImage.IsAtEnd() )
199  {
200  PosteriorsPixelType posteriors(numberOfClasses);
201  const PriorsPixelType priors = itrPriorsImage.Get();
202  const MembershipPixelType memberships = itrMembershipImage.Get();
203  for( unsigned int i=0; i<numberOfClasses; i++)
204  {
205  posteriors[i] =
206  static_cast< TPosteriorsPrecisionType >( memberships[i] * priors[i] );
207  }
208  itrPosteriorsImage.Set( posteriors );
209  ++itrMembershipImage;
210  ++itrPriorsImage;
211  ++itrPosteriorsImage;
212  }
213  }
214  else
215  {
216  PosteriorsImageType * posteriorsImage =
217  dynamic_cast< PosteriorsImageType * >( this->GetPosteriorImage() );
218 
219  if( posteriorsImage == NULL )
220  {
221  itkExceptionMacro("Second output type does not correspond to expected Posteriors Image Type");
222  }
223 
224  InputImageIteratorType itrMembershipImage( membershipImage, imageRegion );
225  PosteriorsImageIteratorType itrPosteriorsImage( posteriorsImage, imageRegion );
226 
227  itrMembershipImage.GoToBegin();
228  itrPosteriorsImage.GoToBegin();
229 
230  while( !itrMembershipImage.IsAtEnd() )
231  {
232  itrPosteriorsImage.Set( itrMembershipImage.Get() );
233  ++itrMembershipImage;
234  ++itrPosteriorsImage;
235  }
236 
237  }
238 }
239 
240 
241 template < class TInputVectorImage, class TLabelsType,
242  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
243 void
244 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
245  TPosteriorsPrecisionType, TPriorsPrecisionType >
246 ::SetSmoothingFilter( SmoothingFilterType * smoothingFilter )
247 {
248  this->m_SmoothingFilter = smoothingFilter;
249  this->m_UserProvidedSmoothingFilter = true;
250  this->Modified();
251 }
252 
253 
257 template < class TInputVectorImage, class TLabelsType,
258  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
259 void
260 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
261  TPosteriorsPrecisionType, TPriorsPrecisionType >
262 ::SetPriors( const PriorsImageType * priors )
263 {
264  this->ProcessObject::SetNthInput(1, const_cast< PriorsImageType * >( priors ) );
265  this->m_UserProvidedPriors = true;
266  this->Modified();
267 }
268 
269 
273 template < class TInputVectorImage, class TLabelsType,
274  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
275 void
276 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
277  TPosteriorsPrecisionType, TPriorsPrecisionType >
278 ::NormalizeAndSmoothPosteriors()
279 {
280 
281  PosteriorsImageIteratorType itrPosteriorImage(
282  this->GetPosteriorImage(), this->GetPosteriorImage()->GetBufferedRegion() );
283 
284  PosteriorsPixelType p;
285  const unsigned int numberOfClasses = this->GetPosteriorImage()->GetVectorLength();
286 
287  for( unsigned int iter=0; iter< m_NumberOfSmoothingIterations; iter++)
288  {
289  itrPosteriorImage.GoToBegin();
290  while( !itrPosteriorImage.IsAtEnd() )
291  {
292  p = itrPosteriorImage.Get();
293 
294  // Normalize P so the probablity across components sums to 1
295  TPosteriorsPrecisionType probability=0;
296  for( unsigned int i=0; i< numberOfClasses; i++ )
297  {
298  probability += p[i];
299  }
300  p /= probability;
301  itrPosteriorImage.Set( p );
302  ++itrPosteriorImage;
303  }
304 
305 
306  for( unsigned int componentToExtract=0; componentToExtract < numberOfClasses; componentToExtract++)
307  {
308  // Create an auxillary image to store one component of the vector image.
309  // Smoothing filters typically can't handle multi-component images, so we
310  // will extract each component and smooth it.
311  typename ExtractedComponentImageType::Pointer extractedComponentImage =
312  ExtractedComponentImageType::New();
313  extractedComponentImage->CopyInformation( this->GetPosteriorImage());
314  extractedComponentImage->SetBufferedRegion(
315  this->GetPosteriorImage()->GetBufferedRegion() );
316  extractedComponentImage->SetRequestedRegion(
317  this->GetPosteriorImage()->GetRequestedRegion() );
318  extractedComponentImage->Allocate();
319  typedef itk::ImageRegionIterator< ExtractedComponentImageType > IteratorType;
320 
321  itrPosteriorImage.GoToBegin();
322  IteratorType it( extractedComponentImage,
323  extractedComponentImage->GetBufferedRegion() );
324 
325  it.GoToBegin();
326  while( !itrPosteriorImage.IsAtEnd() )
327  {
328  it.Set(itrPosteriorImage.Get()[componentToExtract]);
329  ++it;
330  ++itrPosteriorImage;
331  }
332 
333  m_SmoothingFilter->SetInput( extractedComponentImage );
334  m_SmoothingFilter->Modified(); // Force an update
335  m_SmoothingFilter->Update();
336 
337  itrPosteriorImage.GoToBegin();
338 
339  IteratorType sit( m_SmoothingFilter->GetOutput(),
340  m_SmoothingFilter->GetOutput()->GetBufferedRegion() );
341  sit.GoToBegin();
342  while( !itrPosteriorImage.IsAtEnd() )
343  {
344  PosteriorsPixelType posteriorPixel = itrPosteriorImage.Get();
345  posteriorPixel[componentToExtract] = sit.Get();
346  itrPosteriorImage.Set(posteriorPixel);
347  ++sit;
348  ++itrPosteriorImage;
349  }
350  }
351  }
352 }
353 
354 
358 template < class TInputVectorImage, class TLabelsType,
359  class TPosteriorsPrecisionType, class TPriorsPrecisionType >
360 void
361 BayesianClassifierImageFilter<TInputVectorImage, TLabelsType,
362  TPosteriorsPrecisionType, TPriorsPrecisionType >
363 ::ClassifyBasedOnPosteriors()
364 {
365  OutputImagePointer labels = this->GetOutput();
366 
367  ImageRegionType imageRegion = labels->GetBufferedRegion();
368 
369  PosteriorsImageType * posteriorsImage =
370  dynamic_cast< PosteriorsImageType * >( this->GetPosteriorImage() );
371 
372  if( posteriorsImage == NULL )
373  {
374  itkExceptionMacro("Second output type does not correspond to expected Posteriors Image Type");
375  }
376 
377  OutputImageIteratorType itrLabelsImage( labels, imageRegion );
378  PosteriorsImageIteratorType itrPosteriorsImage( posteriorsImage,imageRegion );
379 
380  DecisionRulePointer decisionRule = DecisionRuleType::New();
381 
382  itrLabelsImage.GoToBegin();
383  itrPosteriorsImage.GoToBegin();
384 
385  while ( !itrLabelsImage.IsAtEnd() )
386  {
387  itrLabelsImage.Set( static_cast< TLabelsType >(
388  decisionRule->Evaluate( itrPosteriorsImage.Get())) );
389  ++itrLabelsImage;
390  ++itrPosteriorsImage;
391  }
392 }
393 
394 } // end namespace itk
395 
396 #endif
Generated at Sat Feb 2 2013 23:25:29 for Orfeo Toolbox with doxygen 1.8.1.1