OTB  6.1.0
Orfeo Toolbox
otbLocalHistogramImageFunction.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbLocalHistogramImageFunction_txx
22 #define otbLocalHistogramImageFunction_txx
23 
26 #include "itkExtractImageFilter.h"
27 #include "otbMath.h"
28 
29 namespace otb
30 {
31 
35 template <class TInputImage, class TCoordRep>
38  m_NeighborhoodRadius(1), m_NumberOfHistogramBins(128), m_HistogramMin(0), m_HistogramMax(1), m_GaussianSmoothing(true)
39 {
40 }
41 
42 template <class TInputImage, class TCoordRep>
43 void
45 ::PrintSelf(std::ostream& os, itk::Indent indent) const
46 {
47  this->Superclass::PrintSelf(os, indent);
48  os << indent << " Neighborhood radius value : " << this->GetNeighborhoodRadius() << std::endl;
49  os << indent << " Number Of Histogram Bins : " << this->GetNumberOfHistogramBins() << std::endl;
50  os << indent << " Histogram Minimum : " << this->GetHistogramMin() << std::endl;
51  os << indent << " Histogram Maximum : " << this->GetHistogramMax() << std::endl;
52 }
53 
54 template <class TInputImage, class TCoordRep>
57 ::EvaluateAtIndex(const IndexType& index) const
58 {
59  typename HistogramType::Pointer histogram = HistogramType::New();
60 
61  typename HistogramType::SizeType size(this->GetInputImage()->GetNumberOfComponentsPerPixel());
62  size.Fill(this->GetNumberOfHistogramBins());
63 
64  typename HistogramType::MeasurementVectorType lowerBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
65  typename HistogramType::MeasurementVectorType upperBound(this->GetInputImage()->GetNumberOfComponentsPerPixel());
66 
67  lowerBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMin()) );
68  upperBound.Fill( static_cast<typename HistogramType::MeasurementType>(this->GetHistogramMax()) );
69 
70  histogram->SetMeasurementVectorSize(this->GetInputImage()->GetNumberOfComponentsPerPixel());
71  histogram->Initialize(size, lowerBound, upperBound );
72  histogram->SetToZero();
73 
74  // Check for input image
75  if( !this->GetInputImage() )
76  {
77  return histogram;
78  }
79 
80  // Check for out of buffer
81  if ( !this->IsInsideBuffer( index ) )
82  {
83  return histogram;
84  }
85 
86  // Create an N-d neighborhood kernel, using a zeroflux boundary condition
87  typename InputImageType::SizeType kernelSize;
88  kernelSize.Fill( m_NeighborhoodRadius );
89 
91  it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
92 
93  // Set the iterator at the desired location
94  it.SetLocation(index);
95 
96  // Define a gaussian kernel around the center location
97  double squaredRadius = m_NeighborhoodRadius * m_NeighborhoodRadius;
98  double squaredSigma = 0.25 * squaredRadius;
99 
100  // Offset to be used in the loops
101  typename InputImageType::OffsetType offset;
102 
103  // Fill the histogram
104  for(int i = -(int)m_NeighborhoodRadius; i< (int)m_NeighborhoodRadius; ++i)
105  {
106  for(int j = -(int)m_NeighborhoodRadius; j< (int)m_NeighborhoodRadius; ++j)
107  {
108  // Check if the current pixel lies within a disc of radius m_NeighborhoodRadius
109  double currentSquaredRadius = i*i+j*j;
110  if(currentSquaredRadius < squaredRadius)
111  {
112  // If so, compute the gaussian weighting (this could be
113  // computed once for all for the sake of optimisation) if necessary
114  double gWeight = 1.;
115  if(m_GaussianSmoothing)
116  {
117  gWeight = (1/vcl_sqrt(otb::CONST_2PI*squaredSigma)) * vcl_exp(- currentSquaredRadius/(2*squaredSigma));
118  }
119 
120  // Compute pixel location
121  offset[0]=i;
122  offset[1]=j;
123 
124  // Get the current value
125  typename HistogramType::MeasurementVectorType sample(this->GetInputImage()->GetNumberOfComponentsPerPixel());
126  sample[0] = it.GetPixel(offset);
127 
128  // Populate histogram
129  histogram->IncreaseFrequencyOfMeasurement(sample, gWeight);
130  }
131  }
132  }
133  return histogram;
134 }
135 
136 } // namespace otb
137 
138 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:146
const double CONST_2PI
Definition: otbMath.h:54
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
void Fill(TValue const &v)
void SetLocation(const IndexType &position)
Superclass::MeasurementVectorType MeasurementVectorType
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
OutputType EvaluateAtIndex(const IndexType &index) const ITK_OVERRIDE
void Fill(SizeValueType value)