otbLocalRxDetectorFilter.h 4.77 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20 21


22 23
#ifndef otbLocalRxDetectorFilter_h
#define otbLocalRxDetectorFilter_h
24 25 26 27

#include "itkImageToImageFilter.h"

#include "itkListSample.h"
28
#include "itkCovarianceSampleFilter.h"
29 30
#include "itkConstNeighborhoodIterator.h"
#include "otbVectorImage.h"
31 32 33 34

namespace otb
{

35 36 37 38 39 40 41 42 43 44 45 46
/** \class LocalRxDetectionFunctor
 * \brief This functor computes a local Rx score on an input neighborhood. Pixel of the neighborhood
 * inside the internal radius are not considered during the computation of local statistics.
 *
 * \ingroup ImageFilters
 *
 * \ingroup OTBAnomalyDetection
 */
namespace Functor
{
template<typename TInput, typename TOutput =TInput> 
class LocalRxDetectionFunctor
47 48 49 50
{
public:

  /** typedef */
51 52
  typedef typename itk::Neighborhood<itk::VariableLengthVector<TInput>>::OffsetType OffsetType;
  typedef typename itk::VariableLengthVector<TInput>                         VectorMeasurementType;
53 54 55 56 57 58
  typedef itk::Statistics::ListSample<VectorMeasurementType>            ListSampleType;
  typedef itk::Statistics::CovarianceSampleFilter<ListSampleType>       CovarianceCalculatorType;
  typedef typename CovarianceCalculatorType::MeasurementVectorRealType  MeasurementVectorRealType;
  typedef typename CovarianceCalculatorType::MatrixType                 MatrixType;

private:
59
  // Internal radius along the X axis
60
  unsigned int m_InternalRadiusX;
61

62
  // Internal radius along the Y axis
63
  unsigned int m_InternalRadiusY;
64 65

public:
66
  LocalRxDetectionFunctor() : m_InternalRadiusX(1), m_InternalRadiusY(1) {};
67

68
  void SetInternalRadius(const unsigned int internalRadiusX, const unsigned int internalRadiusY)
69
  {
70 71
    m_InternalRadiusX = internalRadiusX;
    m_InternalRadiusY = internalRadiusY;
72 73
  };

74
  int GetInternalRadiusX()
75
  {
76
    return m_InternalRadiusX;
77 78
  };

79 80 81 82
  int GetInternalRadiusY()
  {
    return m_InternalRadiusY;
  };
83

84
  auto operator()(const itk::ConstNeighborhoodIterator<otb::VectorImage<TInput>> & in) const
85 86 87 88 89 90
  {
    // Create a list sample with the pixels of the neighborhood located between
    // the two radius.
    typename ListSampleType::Pointer listSample = ListSampleType::New();

    // The pixel on whih we will compute the Rx score, we load it now to get the input vector size.
91
    auto centerPixel = in.GetCenterPixel();
92 93 94
    listSample->SetMeasurementVectorSize(centerPixel.Size());

    OffsetType off;
95 96 97 98 99

    // Cache radiuses attributes for threading performances
    const int internalRadiusX = m_InternalRadiusX;
    const int internalRadiusY = m_InternalRadiusY;

100
    auto externalRadius = in.GetRadius();
101
    for (int y = -externalRadius[1]; y <= static_cast<int>(externalRadius[1]); y++)
102 103
      {
      off[1] = y;
104
      for (int x = -externalRadius[0]; x <= static_cast<int>(externalRadius[0]); x++)
105 106
        {
        off[0] = x;
107
        if ((abs(x) > internalRadiusX) || (abs(y) > internalRadiusY))
108
          {
109
            listSample->PushBack(in.GetPixel(off));
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
          }
        }
      }

    // Compute mean & inverse covariance matrix
    typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
    covarianceCalculator->SetInput(listSample);
    covarianceCalculator->Update();

    MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
    
    VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
    for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
      {
      meanVec.SetElement(i, meanVector.GetElement(i));
      }

    MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
    typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();

    typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
    
    for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
      {
      centeredTestPixMat.put(i, 0, (centerPixel.GetElement(i) - meanVector.GetElement(i)));
      }

    // Rx score computation
    typename MatrixType::InternalMatrixType rxValue 
      = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;

141
    return static_cast<TOutput> (rxValue.get(0, 0));
142 143 144 145
  }

};

146
} // end namespace functor
147 148 149
} // end namespace otb

#endif