KullbackLeiblerDistanceChDet.cxx 6.95 KB
 Sébastien Dinot committed Mar 08, 2017 1 /*  Julien Michel committed Jan 14, 2019 2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)  Sébastien Dinot committed Mar 08, 2017 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 committed Sep 17, 2010 20 21   Victor Poughon committed Feb 20, 2019 22 23 24 25 /* Example usage: ./KullbackLeiblerDistanceChDet Input/GomaAvant.png Input/GomaApres.png Output/KLdistanceChDet.png 35 */  Emmanuel Christophe committed Sep 17, 2010 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 committed Jan 31, 2011 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 committed Nov 22, 2015 40 // \right)  OTB Bot committed Jan 31, 2011 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 committed Sep 17, 2010 44 // \left(  OTB Bot committed Jan 31, 2011 45 // c_6 - 6 \frac{c_4}{\kappa_{X_1; 2}} + 9 \frac{c_2}{\kappa_{X_2; 2}^2}  Rashad Kanavath committed Nov 22, 2015 46 // \right)  OTB Bot committed Jan 31, 2011 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 committed Sep 17, 2010 50 51 52 // \end{multline} // where // \begin{align*}  Rashad Kanavath committed Nov 22, 2015 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 committed Jan 31, 2011 55 // a_3 &= c_6 - 15\frac{c_4}{\kappa_{X_2; 2}} + 45\frac{c_2}{\kappa_{X_2; 2}^2} -  Rashad Kanavath committed Nov 22, 2015 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 committed Jan 31, 2011 62 // \beta &= \frac{ \kappa_{X_1; 2}^{1/2} }{\kappa_{X_2; 2}}.  Emmanuel Christophe committed Sep 17, 2010 63 // \end{align*}  OTB Bot committed Jan 31, 2011 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 committed Sep 17, 2010 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.  Emmanuel Christophe committed Jul 22, 2011 76 #include "itkMacro.h"  Emmanuel Christophe committed Sep 17, 2010 77 78 79 #include "otbImage.h" #include "otbImageFileReader.h" #include "otbImageFileWriter.h"  Guillaume Pasero committed May 04, 2015 80 #include "itkUnaryFunctorImageFilter.h"  Emmanuel Christophe committed Sep 17, 2010 81 82 83 84 #include "itkRescaleIntensityImageFilter.h" #include "otbKullbackLeiblerDistanceImageFilter.h"  Victor Poughon committed Feb 20, 2019 85 int main(int argc, char* argv[])  Emmanuel Christophe committed Sep 17, 2010 86 87 { try  Victor Poughon committed Feb 20, 2019 88  {  Emmanuel Christophe committed Sep 17, 2010 89  if (argc != 5)  Victor Poughon committed Feb 20, 2019 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 committed Sep 17, 2010 93 94  std::cerr << argv[0] << " imgAv imgAp imgResu winSize\n"; return 1;  Victor Poughon committed Feb 20, 2019 95  }  Emmanuel Christophe committed Sep 17, 2010 96   Victor Poughon committed Feb 20, 2019 97 98 99 100  char* fileName1 = argv[1]; char* fileName2 = argv[2]; char* fileNameOut = argv[3]; int winSize = atoi(argv[4]);  Emmanuel Christophe committed Sep 17, 2010 101   Victor Poughon committed Feb 20, 2019 102  const unsigned int Dimension = 2;  Emmanuel Christophe committed Sep 17, 2010 103 104 105 106 107 108 109 110 111 112 113 114 115  typedef double PixelType; typedef unsigned char OutputPixelType; typedef otb::Image ImageType; typedef otb::Image 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...  Victor Poughon committed Feb 20, 2019 116  typedef otb::KullbackLeiblerDistanceImageFilter FilterType;  Emmanuel Christophe committed Sep 17, 2010 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 ReaderType; typedef otb::ImageFileWriter 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());  Victor Poughon committed Feb 20, 2019 141 142  typedef itk::RescaleIntensityImageFilter RescaleFilterType; RescaleFilterType::Pointer rescaler = RescaleFilterType::New();  Emmanuel Christophe committed Sep 17, 2010 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();  Victor Poughon committed Feb 20, 2019 152  }  Emmanuel Christophe committed Sep 17, 2010 153 154  catch (itk::ExceptionObject& err)  Victor Poughon committed Feb 20, 2019 155  {  Emmanuel Christophe committed Sep 17, 2010 156 157 158  std::cout << "Exception itk::ExceptionObject thrown !" << std::endl; std::cout << err << std::endl; return EXIT_FAILURE;  Victor Poughon committed Feb 20, 2019 159  }  Emmanuel Christophe committed Sep 17, 2010 160 161  catch (...)  Victor Poughon committed Feb 20, 2019 162  {  Emmanuel Christophe committed Sep 17, 2010 163 164  std::cout << "Unknown exception thrown !" << std::endl; return EXIT_FAILURE;  Victor Poughon committed Feb 20, 2019 165  }  Emmanuel Christophe committed Sep 17, 2010 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; }