OTB  5.0.0
Orfeo Toolbox
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 "itkExtractImageFilter.h"
24 #include "otbMath.h"
25 
26 namespace otb
27 {
28 
32 template <class TInputImage, class TCoordRep>
35  m_NeighborhoodRadius(1), m_NumberOfHistogramBins(128), m_HistogramMin(0), m_HistogramMax(1), m_GaussianSmoothing(true)
36 {
37 }
38 
39 template <class TInputImage, class TCoordRep>
40 void
42 ::PrintSelf(std::ostream& os, itk::Indent indent) const
43 {
44  this->Superclass::PrintSelf(os, indent);
45  os << indent << " Neighborhood radius value : " << this->GetNeighborhoodRadius() << std::endl;
46  os << indent << " Number Of Histogram Bins : " << this->GetNumberOfHistogramBins() << std::endl;
47  os << indent << " Histogram Minimum : " << this->GetHistogramMin() << std::endl;
48  os << indent << " Histogram Maximum : " << this->GetHistogramMax() << std::endl;
49 }
50 
51 template <class TInputImage, class TCoordRep>
54 ::EvaluateAtIndex(const IndexType& index) const
55 {
56  typename HistogramType::Pointer histogram = HistogramType::New();
57 
58  typename HistogramType::SizeType size(this->GetInputImage()->GetNumberOfComponentsPerPixel());
59  size.Fill(this->GetNumberOfHistogramBins());
60 
61  typename HistogramType::MeasurementVectorType lowerBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
62  typename HistogramType::MeasurementVectorType upperBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
63 
64  lowerBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMin()) );
65  upperBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMax()) );
66 
67  histogram->SetMeasurementVectorSize(this->GetInputImage()->GetNumberOfComponentsPerPixel());
68  histogram->Initialize(size, lowerBound, upperBound );
69  histogram->SetToZero();
70 
71  // Check for input image
72  if( !this->GetInputImage() )
73  {
74  return histogram;
75  }
76 
77  // Check for out of buffer
78  if ( !this->IsInsideBuffer( index ) )
79  {
80  return histogram;
81  }
82 
83  // Create an N-d neighborhood kernel, using a zeroflux boundary condition
84  typename InputImageType::SizeType kernelSize;
85  kernelSize.Fill( m_NeighborhoodRadius );
86 
88  it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
89 
90  // Set the iterator at the desired location
91  it.SetLocation(index);
92 
93  // Define a gaussian kernel around the center location
94  double squaredRadius = m_NeighborhoodRadius * m_NeighborhoodRadius;
95  double squaredSigma = 0.25 * squaredRadius;
96 
97  // Offset to be used in the loops
98  typename InputImageType::OffsetType offset;
99 
100  // Fill the histogram
101  for(int i = -(int)m_NeighborhoodRadius; i< (int)m_NeighborhoodRadius; ++i)
102  {
103  for(int j = -(int)m_NeighborhoodRadius; j< (int)m_NeighborhoodRadius; ++j)
104  {
105  // Check if the current pixel lies within a disc of radius m_NeighborhoodRadius
106  double currentSquaredRadius = i*i+j*j;
107  if(currentSquaredRadius < squaredRadius)
108  {
109  // If so, compute the gaussian weighting (this could be
110  // computed once for all for the sake of optimisation) if necessary
111  double gWeight = 1.;
112  if(m_GaussianSmoothing)
113  {
114  gWeight = (1/vcl_sqrt(otb::CONST_2PI*squaredSigma)) * vcl_exp(- currentSquaredRadius/(2*squaredSigma));
115  }
116 
117  // Compute pixel location
118  offset[0]=i;
119  offset[1]=j;
120 
121  // Get the current value
122  typename HistogramType::MeasurementVectorType sample(this->GetInputImage()->GetNumberOfComponentsPerPixel());
123  sample[0] = it.GetPixel(offset);
124 
125  // Populate histogram
126  histogram->IncreaseFrequencyOfMeasurement(sample, gWeight);
127  }
128  }
129  }
130  return histogram;
131 }
132 
133 } // namespace otb
134 
135 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const
const double CONST_2PI
Definition: otbMath.h:51
void Fill(TValue const &v)
void SetLocation(const IndexType &position)
Superclass::MeasurementVectorType MeasurementVectorType
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
virtual OutputType EvaluateAtIndex(const IndexType &index) const