KullbackLeiblerDistanceChDet.cxx 6.95 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.
 */
Emmanuel Christophe's avatar
Emmanuel Christophe committed
20 21


22 23 24 25
/* Example usage:
./KullbackLeiblerDistanceChDet Input/GomaAvant.png Input/GomaApres.png Output/KLdistanceChDet.png 35
*/

Emmanuel Christophe's avatar
Emmanuel Christophe committed
26 27 28 29 30 31 32 33 34 35

// This example illustrates the class
// \doxygen{otb}{KullbackLeiblerDistanceImageFilter} for detecting changes
// between pairs of images. This filter computes the Kullback-Leibler
// distance between probability density functions (pdfs).
// In fact, the Kullback-Leibler distance is itself approximated through
// a cumulant-based expansion, since the pdfs are approximated through an
// Edgeworth series.
// The Kullback-Leibler distance is evaluated by:
// \begin{multline}\label{eqKLapprox1D}
OTB Bot's avatar
STYLE  
OTB Bot committed
36 37 38 39
//    K_{\text{Edgeworth}}(X_1 | X_2) = \frac{1}{12} \frac{\kappa_{X_1; 3}^2}{\kappa_{X_1; 2}^2}
//        + \frac{1}{2} \left( \log \frac{\kappa_{X_2; 2}}{\kappa_{X_1; 2}}
//                             -1+\frac{1}{\kappa_{X_2; 2}}
//                             \left( \kappa_{X_1; 1} - \kappa_{X_2; 1} +  \kappa_{X_1; 2}^{1/2} \right)^2
Rashad Kanavath's avatar
Rashad Kanavath committed
40
//                        \right)
OTB Bot's avatar
STYLE  
OTB Bot committed
41 42 43
//        - \left( \kappa_{X_2; 3} \frac{a_1}{6} + \kappa_{X_2; 4} \frac{a_2}{24}
//            + \kappa_{X_2; 3}^2 \frac{a_3}{72} \right)
//        - \frac{1}{2} \frac{ \kappa_{X_2; 3}^2}{36}
Emmanuel Christophe's avatar
Emmanuel Christophe committed
44
//            \left(
OTB Bot's avatar
STYLE  
OTB Bot committed
45
//                c_6 - 6 \frac{c_4}{\kappa_{X_1; 2}} + 9 \frac{c_2}{\kappa_{X_2; 2}^2}
Rashad Kanavath's avatar
Rashad Kanavath committed
46
//            \right)
OTB Bot's avatar
STYLE  
OTB Bot committed
47 48 49
//        - 10 \frac{\kappa_{X_1; 3} \kappa_{X_2; 3}
//                        \left( \kappa_{X_1; 1} - \kappa_{X_2; 1} \right)
//                        \left( \kappa_{X_1; 2} - \kappa_{X_2; 2} \right)}{\kappa_{X_2; 2}^6} \qquad
Emmanuel Christophe's avatar
Emmanuel Christophe committed
50 51 52
// \end{multline}
// where
// \begin{align*}
Rashad Kanavath's avatar
Rashad Kanavath committed
53 54
//    a_1 &= c_3 - 3 \frac{\alpha}{\kappa_{X_2; 2}}
//    a_2 &= c_4 - 6 \frac{c_2}{\kappa_{X_2; 2}} + \frac{3}{\kappa_{X_2; 2}^2}
OTB Bot's avatar
STYLE  
OTB Bot committed
55
//    a_3 &= c_6 - 15\frac{c_4}{\kappa_{X_2; 2}} + 45\frac{c_2}{\kappa_{X_2; 2}^2} -
Rashad Kanavath's avatar
Rashad Kanavath committed
56 57 58 59 60 61
//            \frac{15}{\kappa_{X_2; 2}^3}
//    c_2 &= \alpha^2 + \beta^2
//    c_3 &= \alpha^3 + 3 \alpha \beta^2
//    c_4 &= \alpha^4 + 6 \alpha^2 \beta^2 + 3 \beta^4
//    c_6 &= \alpha^6 + 15\alpha^4 \beta^2 + 45 \alpha^2 \beta^4 + 15 \beta^6
//    \alpha &= \frac{\kappa_{X_1; 1} - \kappa_{X_2; 1}}{\kappa_{X_2; 2}}
OTB Bot's avatar
STYLE  
OTB Bot committed
62
//    \beta &= \frac{ \kappa_{X_1; 2}^{1/2} }{\kappa_{X_2; 2}}.
Emmanuel Christophe's avatar
Emmanuel Christophe committed
63
// \end{align*}
OTB Bot's avatar
STYLE  
OTB Bot committed
64 65
// $\kappa_{X_i; 1}$, $\kappa_{X_i; 2}$, $\kappa_{X_i; 3}$ and $\kappa_{X_i; 4}$
// are the cumulants up to order 4 of the random variable $X_i$ ($i=1, 2$).
Emmanuel Christophe's avatar
Emmanuel Christophe committed
66 67 68 69 70 71 72 73 74 75
// This example will use the images shown in
// figure~\ref{fig:RATCHDETINIM}. These correspond to 2 Radarsat fine
// mode acquisitions before and after a lava flow resulting from a
// volcanic eruption.
//
// The program itself is very similar to the ratio of means detector,
// implemented in \doxygen{otb}{MeanRatioImageFilter},
// in section~\ref{sec:RatioOfMeans}. Nevertheless
// the corresponding header file has to be used instead.

76
#include "itkMacro.h"
Emmanuel Christophe's avatar
Emmanuel Christophe committed
77 78 79
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
80
#include "itkUnaryFunctorImageFilter.h"
Emmanuel Christophe's avatar
Emmanuel Christophe committed
81 82 83 84
#include "itkRescaleIntensityImageFilter.h"

#include "otbKullbackLeiblerDistanceImageFilter.h"

85
int main(int argc, char* argv[])
Emmanuel Christophe's avatar
Emmanuel Christophe committed
86 87
{
  try
88
  {
Emmanuel Christophe's avatar
Emmanuel Christophe committed
89
    if (argc != 5)
90 91 92
    {
      std::cerr << "Change detection through a Kullback-Leibler measure (which is a distance between local distributions)\n";
      std::cerr << "Kullback-Leibler measure is optimized by a Edgeworth series expansion\n";
Emmanuel Christophe's avatar
Emmanuel Christophe committed
93 94
      std::cerr << argv[0] << " imgAv imgAp imgResu winSize\n";
      return 1;
95
    }
Emmanuel Christophe's avatar
Emmanuel Christophe committed
96

97 98 99 100
    char* fileName1   = argv[1];
    char* fileName2   = argv[2];
    char* fileNameOut = argv[3];
    int   winSize     = atoi(argv[4]);
Emmanuel Christophe's avatar
Emmanuel Christophe committed
101

102
    const unsigned int    Dimension = 2;
Emmanuel Christophe's avatar
Emmanuel Christophe committed
103 104 105 106 107 108 109 110 111 112 113 114 115
    typedef double        PixelType;
    typedef unsigned char OutputPixelType;

    typedef otb::Image<PixelType, Dimension>       ImageType;
    typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

    //  The \doxygen{otb}{KullbackLeiblerDistanceImageFilter} is templated over
    //  the types of the two input images and the type of the generated change
    //  image, in a similar way as the \doxygen{otb}{MeanRatioImageFilter}. It is
    //  the only line to be changed from the ratio of means change detection
    //  example to perform a change detection through a distance between
    //  distributions...

116
    typedef otb::KullbackLeiblerDistanceImageFilter<ImageType, ImageType, ImageType> FilterType;
Emmanuel Christophe's avatar
Emmanuel Christophe committed
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

    //  The different elements of the pipeline can now be instantiated. Follow the
    //  ratio of means change detector example.

    typedef otb::ImageFileReader<ImageType>       ReaderType;
    typedef otb::ImageFileWriter<OutputImageType> WriterType;

    ReaderType::Pointer reader1 = ReaderType::New();
    reader1->SetFileName(fileName1);

    ReaderType::Pointer reader2 = ReaderType::New();
    reader2->SetFileName(fileName2);

    //  The only parameter for this change detector is the radius of
    //  the window used for computing the cumulants.

    FilterType::Pointer filter = FilterType::New();
    filter->SetRadius((winSize - 1) / 2);

    //  The pipeline is built by plugging all the elements together.

    filter->SetInput1(reader1->GetOutput());
    filter->SetInput2(reader2->GetOutput());

141 142
    typedef itk::RescaleIntensityImageFilter<ImageType, OutputImageType> RescaleFilterType;
    RescaleFilterType::Pointer                                           rescaler = RescaleFilterType::New();
Emmanuel Christophe's avatar
Emmanuel Christophe committed
143 144 145 146 147 148 149 150 151

    rescaler->SetInput(filter->GetOutput());
    rescaler->SetOutputMinimum(0);
    rescaler->SetOutputMaximum(255);

    WriterType::Pointer writer = WriterType::New();
    writer->SetFileName(fileNameOut);
    writer->SetInput(rescaler->GetOutput());
    writer->Update();
152
  }
Emmanuel Christophe's avatar
Emmanuel Christophe committed
153 154

  catch (itk::ExceptionObject& err)
155
  {
Emmanuel Christophe's avatar
Emmanuel Christophe committed
156 157 158
    std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
    std::cout << err << std::endl;
    return EXIT_FAILURE;
159
  }
Emmanuel Christophe's avatar
Emmanuel Christophe committed
160 161

  catch (...)
162
  {
Emmanuel Christophe's avatar
Emmanuel Christophe committed
163 164
    std::cout << "Unknown exception thrown !" << std::endl;
    return EXIT_FAILURE;
165
  }
Emmanuel Christophe's avatar
Emmanuel Christophe committed
166 167 168 169 170 171 172 173 174 175 176 177 178 179

  // Figure \ref{fig:RESKLDCHDET} shows the result of the change
  // detection by computing the Kullback-Leibler distance between
  // local pdf through an Edgeworth approximation.
  // \begin{figure}
  // \center
  // \includegraphics[width=0.35\textwidth]{KLdistanceChDet.eps}
  // \itkcaption[Kullback-Leibler Change Detection Results]{Result of the
  // Kullback-Leibler change detector}
  // \label{fig:RESKLDCHDET}
  // \end{figure}

  return EXIT_SUCCESS;
}