Orfeo Toolbox  3.16
itkHistogramMatchingImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkHistogramMatchingImageFilter.txx,v $
5  Language: C++
6  Date: $Date: 2009-05-05 17:00:00 $
7  Version: $Revision: 1.19 $
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  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkHistogramMatchingImageFilter_txx
18 #define __itkHistogramMatchingImageFilter_txx
19 
20 #include "itkHistogramMatchingImageFilter.h"
21 #include "itkImageRegionIterator.h"
22 #include "itkImageRegionConstIterator.h"
23 #include "itkNumericTraits.h"
24 #include <vector>
25 
26 namespace itk
27 {
28 
32 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
33 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
34 ::HistogramMatchingImageFilter()
35 {
36  this->SetNumberOfRequiredInputs( 2 );
37 
38  m_NumberOfHistogramLevels = 256;
39  m_NumberOfMatchPoints = 1;
40 
41  m_QuantileTable.set_size( 3, m_NumberOfMatchPoints + 2 );
42  m_QuantileTable.fill(0);
43  m_Gradients.set_size( m_NumberOfMatchPoints + 1 );
44  m_Gradients.fill(0);
45 
46  m_ThresholdAtMeanIntensity = true;
47  m_SourceIntensityThreshold = 0;
48  m_ReferenceIntensityThreshold = 0;
49  m_LowerGradient = 0.0;
50  m_UpperGradient = 0.0;
51 
52  // Create histograms.
53  m_SourceHistogram = HistogramType::New();
54  m_ReferenceHistogram = HistogramType::New();
55  m_OutputHistogram = HistogramType::New();
56 }
57 
58 
59 /*
60  *
61  */
62 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
63 void
64 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
65 ::PrintSelf(std::ostream& os, Indent indent) const
66 {
67  Superclass::PrintSelf(os,indent);
68 
69  os << indent << "NumberOfHistogramLevels: ";
70  os << m_NumberOfHistogramLevels << std::endl;
71  os << indent << "NumberOfMatchPoints: ";
72  os << m_NumberOfMatchPoints << std::endl;
73  os << indent << "ThresholdAtMeanIntensity: ";
74  os << m_ThresholdAtMeanIntensity << std::endl;
75 
76  os << indent << "SourceIntensityThreshold: ";
77  os << m_SourceIntensityThreshold << std::endl;
78  os << indent << "ReferenceIntensityThreshold: ";
79  os << m_ReferenceIntensityThreshold << std::endl;
80  os << indent << "OutputIntensityThreshold: ";
81  os << m_ReferenceIntensityThreshold << std::endl;
82  os << indent << "Source histogram: ";
83  os << m_SourceHistogram.GetPointer() << std::endl;
84  os << indent << "Reference histogram: ";
85  os << m_ReferenceHistogram.GetPointer() << std::endl;
86  os << indent << "Output histogram: ";
87  os << m_OutputHistogram.GetPointer() << std::endl;
88  os << indent << "QuantileTable: " << std::endl;
89  os << m_QuantileTable << std::endl;
90  os << indent << "Gradients: " << std::endl;
91  os << m_Gradients << std::endl;
92  os << indent << "LowerGradient: ";
93  os << m_LowerGradient << std::endl;
94  os << indent << "UpperGradient: ";
95  os << m_UpperGradient << std::endl;
96 }
97 
98 
99 /*
100  *
101  */
102 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
103 void
104 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
105 ::SetReferenceImage( const InputImageType * reference )
106 {
107  this->ProcessObject::SetNthInput(1,
108  const_cast< InputImageType * >( reference ) );
109 }
110 
111 
112 /*
113  *
114  */
115 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
116 const typename HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
117 ::InputImageType *
118 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
119 ::GetReferenceImage()
120 {
121  if ( this->GetNumberOfInputs() < 2 )
122  {
123  return NULL;
124  }
125 
126  return dynamic_cast<TInputImage*>(
127  this->ProcessObject::GetInput(1) );
128 }
129 
130 
131 /*
132  * This filter requires all of the input images to be
133  * in the buffer.
134  */
135 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
136 void
137 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
138 ::GenerateInputRequestedRegion()
139 {
140  this->Superclass::GenerateInputRequestedRegion();
141 
142  for ( unsigned int idx = 0; idx < this->GetNumberOfInputs(); ++idx )
143  {
144  if ( this->GetInput(idx) )
145  {
146  InputImagePointer image =
147  const_cast< InputImageType * >( this->GetInput(idx) );
148  image->SetRequestedRegionToLargestPossibleRegion();
149  }
150  }
151 }
152 
153 
157 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
158 void
159 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
160 ::BeforeThreadedGenerateData()
161 {
162  unsigned int j;
163 
164  InputImageConstPointer source = this->GetSourceImage();
165  InputImageConstPointer reference = this->GetReferenceImage();
166 
167  this->ComputeMinMaxMean( source, m_SourceMinValue,
168  m_SourceMaxValue, m_SourceMeanValue );
169  this->ComputeMinMaxMean( reference, m_ReferenceMinValue,
170  m_ReferenceMaxValue, m_ReferenceMeanValue );
171 
172  if ( m_ThresholdAtMeanIntensity )
173  {
174  m_SourceIntensityThreshold = static_cast<InputPixelType>(m_SourceMeanValue);
175  m_ReferenceIntensityThreshold = static_cast<InputPixelType>(m_ReferenceMeanValue);
176  }
177  else
178  {
179  m_SourceIntensityThreshold = static_cast<InputPixelType>(m_SourceMinValue);
180  m_ReferenceIntensityThreshold = static_cast<InputPixelType>(m_ReferenceMinValue);
181  }
182 
183  this->ConstructHistogram( source, m_SourceHistogram,
184  m_SourceIntensityThreshold, m_SourceMaxValue );
185  this->ConstructHistogram( reference, m_ReferenceHistogram,
186  m_ReferenceIntensityThreshold,
187  m_ReferenceMaxValue );
188 
189  // Fill in the quantile table.
190  m_QuantileTable.set_size( 3, m_NumberOfMatchPoints + 2 );
191  m_QuantileTable[0][0] = m_SourceIntensityThreshold;
192  m_QuantileTable[1][0] = m_ReferenceIntensityThreshold;
193 
194  m_QuantileTable[0][m_NumberOfMatchPoints + 1] = m_SourceMaxValue;
195  m_QuantileTable[1][m_NumberOfMatchPoints + 1] = m_ReferenceMaxValue;
196 
197  double delta = 1.0 / ( double(m_NumberOfMatchPoints) + 1.0 );
198 
199  for ( j = 1; j < m_NumberOfMatchPoints + 1; j++ )
200  {
201  m_QuantileTable[0][j] = m_SourceHistogram->Quantile(
202  0, double(j) * delta );
203  m_QuantileTable[1][j] = m_ReferenceHistogram->Quantile(
204  0, double(j) * delta );
205  }
206 
207  // Fill in the gradient array.
208  m_Gradients.set_size( m_NumberOfMatchPoints + 1 );
209  double denominator;
210  for ( j = 0; j < m_NumberOfMatchPoints + 1; j++ )
211  {
212  denominator = m_QuantileTable[0][j+1] -
213  m_QuantileTable[0][j];
214  if ( denominator != 0 )
215  {
216  m_Gradients[j] = m_QuantileTable[1][j+1] -
217  m_QuantileTable[1][j];
218  m_Gradients[j] /= denominator;
219  }
220  else
221  {
222  m_Gradients[j] = 0.0;
223  }
224  }
225 
226  denominator = m_QuantileTable[0][0] - m_SourceMinValue;
227  if ( denominator != 0 )
228  {
229  m_LowerGradient = m_QuantileTable[1][0] - m_ReferenceMinValue;
230  m_LowerGradient /= denominator;
231  }
232  else
233  {
234  m_LowerGradient = 0.0;
235  }
236 
237  denominator = m_QuantileTable[0][m_NumberOfMatchPoints+1] -
238  m_SourceMaxValue;
239  if ( denominator != 0 )
240  {
241  m_UpperGradient = m_QuantileTable[1][m_NumberOfMatchPoints+1] -
242  m_ReferenceMaxValue;
243  m_UpperGradient /= denominator;
244  }
245  else
246  {
247  m_UpperGradient = 0.0;
248  }
249 }
250 
251 
255 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
256 void
257 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
258 ::AfterThreadedGenerateData()
259 {
260 
261  OutputImagePointer output = this->GetOutput();
262 
263  this->ComputeMinMaxMean( output, m_OutputMinValue,
264  m_OutputMaxValue, m_OutputMeanValue );
265 
266  if ( m_ThresholdAtMeanIntensity )
267  {
268  m_OutputIntensityThreshold = static_cast<OutputPixelType>(m_OutputMeanValue);
269  }
270  else
271  {
272  m_OutputIntensityThreshold = static_cast<OutputPixelType>(m_OutputMinValue);
273  }
274 
275  this->ConstructHistogram( output, m_OutputHistogram,
276  m_OutputIntensityThreshold, m_OutputMaxValue );
277 
278  // Fill in the quantile table.
279  m_QuantileTable[2][0] = m_OutputIntensityThreshold;
280 
281  m_QuantileTable[2][m_NumberOfMatchPoints + 1] = m_OutputMaxValue;
282 
283  double delta = 1.0 / ( double(m_NumberOfMatchPoints) + 1.0 );
284 
285  for ( unsigned int j = 1; j < m_NumberOfMatchPoints + 1; j++ )
286  {
287  m_QuantileTable[2][j] = m_OutputHistogram->Quantile(
288  0, double(j) * delta );
289  }
290 }
291 
292 
296 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
297 void
298 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
299 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
300  int threadId )
301 {
302 
303  int i;
304  unsigned int j;
305 
306  // Get the input and output pointers;
307  InputImageConstPointer input = this->GetInput();
308  OutputImagePointer output = this->GetOutput();
309 
310 
311  // Transform the source image and write to output.
312  typedef ImageRegionConstIterator<InputImageType> InputConstIterator;
313  typedef ImageRegionIterator<OutputImageType> OutputIterator;
314 
315  InputConstIterator inIter( input, outputRegionForThread );
316  OutputIterator outIter( output, outputRegionForThread );
317 
318 
319  // support progress methods/callbacks
320  unsigned long updateVisits = 0;
321  unsigned long totalPixels = 0;
322  if ( threadId == 0 )
323  {
324  totalPixels = outputRegionForThread.GetNumberOfPixels();
325  updateVisits = totalPixels / 10;
326  if( updateVisits < 1 ) updateVisits = 1;
327  }
328 
329 
330  double srcValue, mappedValue;
331 
332  for ( i = 0; !outIter.IsAtEnd(); ++inIter, ++outIter, i++ )
333  {
334 
335  if ( threadId == 0 && !(i % updateVisits ) )
336  {
337  this->UpdateProgress((float)i / (float)totalPixels);
338  }
339 
340  srcValue = static_cast<double>( inIter.Get() );
341 
342  for ( j = 0; j < m_NumberOfMatchPoints + 2; j++ )
343  {
344  if ( srcValue < m_QuantileTable[0][j] )
345  {
346  break;
347  }
348  }
349 
350  if ( j == 0 )
351  {
352  // Linear interpolate from min to point[0]
353  mappedValue = m_ReferenceMinValue +
354  ( srcValue - m_SourceMinValue ) * m_LowerGradient;
355  }
356  else if ( j == m_NumberOfMatchPoints + 2 )
357  {
358  // Linear interpolate from point[m_NumberOfMatchPoints+1] to max
359  mappedValue = m_ReferenceMaxValue +
360  ( srcValue - m_SourceMaxValue ) * m_UpperGradient;
361  }
362  else
363  {
364  // Linear interpolate from point[j] and point[j+1].
365  mappedValue = m_QuantileTable[1][j-1] +
366  ( srcValue - m_QuantileTable[0][j-1] ) * m_Gradients[j-1];
367  }
368 
369  outIter.Set( static_cast<OutputPixelType>( mappedValue ) );
370 
371  }
372 
373 }
374 
375 
379 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
380 void
381 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
382 ::ComputeMinMaxMean(
383  const InputImageType * image,
384  THistogramMeasurement& minValue,
385  THistogramMeasurement& maxValue,
386  THistogramMeasurement& meanValue )
387 {
388  typedef ImageRegionConstIterator<InputImageType> ConstIterator;
389  ConstIterator iter( image, image->GetBufferedRegion() );
390 
391  double sum = 0.0;
392  long int count = 0;
393 
394  minValue = static_cast<THistogramMeasurement>( iter.Get() );
395  maxValue = minValue;
396 
397  while ( !iter.IsAtEnd() )
398  {
399  const THistogramMeasurement value = static_cast<THistogramMeasurement>( iter.Get() );
400  sum += static_cast<double>(value);
401 
402  if ( value < minValue ) { minValue = value; }
403  if ( value > maxValue ) { maxValue = value; }
404 
405  ++iter;
406  ++count;
407 
408  }
409 
410  meanValue = static_cast<THistogramMeasurement>( sum / static_cast<double>(count) );
411 
412 }
413 
417 template <class TInputImage, class TOutputImage, class THistogramMeasurement>
418 void
419 HistogramMatchingImageFilter<TInputImage,TOutputImage,THistogramMeasurement>
420 ::ConstructHistogram(
421  const InputImageType * image,
422  HistogramType * histogram,
423  const THistogramMeasurement minValue,
424  const THistogramMeasurement maxValue )
425 {
426  {
427  // allocate memory for the histogram
428  typename HistogramType::SizeType size;
429  typename HistogramType::MeasurementVectorType lowerBound;
430  typename HistogramType::MeasurementVectorType upperBound;
431 
432 #ifdef ITK_USE_REVIEW_STATISTICS
433  size.SetSize(1);
434  lowerBound.SetSize(1);
435  upperBound.SetSize(1);
436  histogram->SetMeasurementVectorSize(1);
437 #endif
438 
439  size[0] = m_NumberOfHistogramLevels;
440  lowerBound.Fill(minValue);
441  upperBound.Fill(maxValue);
442 
443  //Initialize with equally spaced bins.
444  histogram->Initialize( size , lowerBound, upperBound );
445  histogram->SetToZero();
446  }
447 
448  typename HistogramType::MeasurementVectorType measurement;
449 
450 #ifdef ITK_USE_REVIEW_STATISTICS
451  measurement.SetSize(1);
452 #endif
453 
454  typedef typename HistogramType::MeasurementType MeasurementType;
455  measurement[0] = NumericTraits<MeasurementType>::Zero;
456 
457  {
458  // put each image pixel into the histogram
459  typedef ImageRegionConstIterator<InputImageType> ConstIterator;
460  ConstIterator iter( image, image->GetBufferedRegion() );
461 
462  iter.GoToBegin();
463  while ( !iter.IsAtEnd() )
464  {
465  InputPixelType value = iter.Get();
466 
467  if ( static_cast<double>(value) >= minValue &&
468  static_cast<double>(value) <= maxValue )
469  {
470  // add sample to histogram
471  measurement[0] = value;
472  histogram->IncreaseFrequency( measurement, 1 );
473  }
474  ++iter;
475  }
476  }
477 }
478 
479 
480 } // end namespace itk
481 
482 #endif
Generated at Sat Feb 2 2013 23:41:39 for Orfeo Toolbox with doxygen 1.8.1.1