Orfeo Toolbox  4.0
otbLocalHistogramImageFunction.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbLocalHistogramImageFunction_txx
19 #define __otbLocalHistogramImageFunction_txx
20 
23 #include "itkNumericTraits.h"
24 #include "itkMacro.h"
25 #include "itkExtractImageFilter.h"
26 #include "otbMath.h"
27 
28 namespace otb
29 {
30 
34 template <class TInputImage, class TCoordRep>
37  m_NeighborhoodRadius(1), m_NumberOfHistogramBins(128), m_HistogramMin(0), m_HistogramMax(1), m_GaussianSmoothing(true)
38 {
39 }
40 
41 template <class TInputImage, class TCoordRep>
42 void
44 ::PrintSelf(std::ostream& os, itk::Indent indent) const
45 {
46  this->Superclass::PrintSelf(os, indent);
47  os << indent << " Neighborhood radius value : " << this->GetNeighborhoodRadius() << std::endl;
48  os << indent << " Number Of Histogram Bins : " << this->GetNumberOfHistogramBins() << std::endl;
49  os << indent << " Histogram Minimum : " << this->GetHistogramMin() << std::endl;
50  os << indent << " Histogram Maximum : " << this->GetHistogramMax() << std::endl;
51 }
52 
53 template <class TInputImage, class TCoordRep>
56 ::EvaluateAtIndex(const IndexType& index) const
57 {
58  typename HistogramType::Pointer histogram = HistogramType::New();
59 
60  typename HistogramType::SizeType size(this->GetInputImage()->GetNumberOfComponentsPerPixel());
61  size.Fill(this->GetNumberOfHistogramBins());
62 
63  typename HistogramType::MeasurementVectorType lowerBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
64  typename HistogramType::MeasurementVectorType upperBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
65 
66  lowerBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMin()) );
67  upperBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMax()) );
68 
69  histogram->SetMeasurementVectorSize(this->GetInputImage()->GetNumberOfComponentsPerPixel());
70  histogram->Initialize(size, lowerBound, upperBound );
71  histogram->SetToZero();
72 
73  // Check for input image
74  if( !this->GetInputImage() )
75  {
76  return histogram;
77  }
78 
79  // Check for out of buffer
80  if ( !this->IsInsideBuffer( index ) )
81  {
82  return histogram;
83  }
84 
85  // Create an N-d neighborhood kernel, using a zeroflux boundary condition
86  typename InputImageType::SizeType kernelSize;
87  kernelSize.Fill( m_NeighborhoodRadius );
88 
90  it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
91 
92  // Set the iterator at the desired location
93  it.SetLocation(index);
94 
95  // Define a gaussian kernel around the center location
96  double squaredRadius = m_NeighborhoodRadius * m_NeighborhoodRadius;
97  double squaredSigma = 0.25 * squaredRadius;
98 
99  // Offset to be used in the loops
100  typename InputImageType::OffsetType offset;
101 
102  // Fill the histogram
103  for(int i = -(int)m_NeighborhoodRadius; i< (int)m_NeighborhoodRadius; ++i)
104  {
105  for(int j = -(int)m_NeighborhoodRadius; j< (int)m_NeighborhoodRadius; ++j)
106  {
107  // Check if the current pixel lies within a disc of radius m_NeighborhoodRadius
108  double currentSquaredRadius = i*i+j*j;
109  if(currentSquaredRadius < squaredRadius)
110  {
111  // If so, compute the gaussian weighting (this could be
112  // computed once for all for the sake of optimisation) if necessary
113  double gWeight = 1.;
114  if(m_GaussianSmoothing)
115  {
116  gWeight = (1/vcl_sqrt(otb::CONST_2PI*squaredSigma)) * vcl_exp(- currentSquaredRadius/(2*squaredSigma));
117  }
118 
119  // Compute pixel location
120  offset[0]=i;
121  offset[1]=j;
122 
123  // Get the current value
124  typename HistogramType::MeasurementVectorType sample(this->GetInputImage()->GetNumberOfComponentsPerPixel());
125  sample[0] = it.GetPixel(offset);
126 
127  // Populate histogram
128  histogram->IncreaseFrequencyOfMeasurement(sample, gWeight);
129  }
130  }
131  }
132  return histogram;
133 }
134 
135 } // namespace otb
136 
137 #endif

Generated at Sat Mar 8 2014 16:06:01 for Orfeo Toolbox with doxygen 1.8.3.1