Commit 0564820e authored by Julien Michel's avatar Julien Michel

DOC: Adding an example for MAD filter (work in progress)

parent 86caf57f
......@@ -38,6 +38,8 @@ TARGET_LINK_LIBRARIES(KullbackLeiblerProfileChDet OTBIO OTBCommon ITKStatistics
ADD_EXECUTABLE(KullbackLeiblerSupervizedDistanceChDet KullbackLeiblerSupervizedDistanceChDet.cxx)
TARGET_LINK_LIBRARIES(KullbackLeiblerSupervizedDistanceChDet OTBIO OTBCommon ITKStatistics)
ADD_EXECUTABLE(MultivariateAlterationDetector otbMultivariateAlterationDetectorExample.cxx)
TARGET_LINK_LIBRARIES(MultivariateAlterationDetector OTBIO OTBCommon ITKStatistics)
IF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
......
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbStreamingImageFileWriter.h"
#include "otbPrintableImageFilter.h"
// Software Guide : BeginCommandLineArgs
// INPUTS: {/Spot5-Gloucester-before.tif}, {Spot5-Gloucester-after.tif}
// OUTPUTS: {MADOutput.tif}, {mad-input1.png}, {mad-input2.png}, {mad-output.png}
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex This example illustrates the class
// \doxygen{otb}{MultivariateAlterationChangeDetectorImageFilter},
// which implements the Multivariate Alteration Change Detector
// algorithm (TODO: Add citation). This algorihtm allows to perform
// change detection from a pair multi-band images, including images
// with different number of bands or modalities. Its output is a a
// multi-band image of change maps, each one being unccorrelated with
// the remaining. The number of bands of the output image is the
// minimum number of bands between the two input images.
//
// The algorithm works as follows. It tries to find two linear
// combinations of bands (one for each input images) which maximize
// correlation, and subtract these two linear combinitation, leading
// to the first change map. Then, it looks for a second set of linear
// combinations which are orthogonal to the first ones, a which
// maximize correlation, and use it as the second change map. This
// process is iterated until no more orthogonal linear combinations
// can be found.
//
// This algorithms has numerous advantages, such as radiometry scaling
// and shifting invariance and absence of parameters, but it can not
// be used on a pair of single band images (in this case the output is
// simply the difference between the two images).
//
// We start by including the corresponding header file.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "otbMultivariateAlterationDetectorImageFilter.h"
// Software Guide : EndCodeSnippet
int main(int argc, char* argv[])
{
if (argc < 4)
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " inputImageFile1 inputImageFile2 outIn1Pretty outIn2Pretty outPretty"
<< "outputImageFile" << std::endl;
return -1;
}
// Define the dimension of the images
const unsigned int Dimension = 2;
// Software Guide : BeginLatex
// We then define the types for the input images, of the
// change image and of the image to be stored in a file for visualization.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef unsigned short InputPixelType;
typedef float OutputPixelType;
typedef otb::VectorImage<InputPixelType, Dimension> InputImageType;
typedef otb::VectorImage<OutputPixelType, Dimension> OutputImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We can now declare the types for the reader. Since the images
// can be vey large, we will force the pipeline to use
// streaming. For this purpose, the file writer will be
// streamed. This is achieved by using the
// \doxygen{otb}{StreamingImageFileWriter} class.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef otb::ImageFileReader<InputImageType> ReaderType;
typedef otb::StreamingImageFileWriter<OutputImageType> WriterType;
// Software Guide : EndCodeSnippet
// This is for rendering in software guide
typedef otb::PrintableImageFilter<InputImageType,InputImageType> InputPrintFilterType;
typedef otb::PrintableImageFilter<OutputImageType,OutputImageType> OutputPrintFilterType;
typedef InputPrintFilterType::OutputImageType VisuImageType;
typedef otb::StreamingImageFileWriter<VisuImageType> VisuWriterType;
// The \doxygen{otb}{MultivariateAlterationDetectorImageFilter} is templated over
// the type of the input images and the type of the generated change
// image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef otb::MultivariateAlterationDetectorImageFilter<
InputImageType,OutputImageType> MADFilterType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The different elements of the pipeline can now be instantiated.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ReaderType::Pointer reader1 = ReaderType::New();
ReaderType::Pointer reader2 = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
MADFilterType::Pointer madFilter = MADFilterType::New();
// Software Guide : EndCodeSnippet
const char * inputFilename1 = argv[1];
const char * inputFilename2 = argv[2];
const char * outputFilename = argv[3];
const char * in1pretty = argv[4];
const char * in2pretty = argv[5];
const char * outpretty = argv[6];
// Software Guide : BeginLatex
//
// We set the parameters of the different elements of the pipeline.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
reader1->SetFileName(inputFilename1);
reader2->SetFileName(inputFilename2);
writer->SetFileName(outputFilename);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We build the pipeline by plugging all the elements together.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
madFilter->SetInput1(reader1->GetOutput());
madFilter->SetInput2(reader2->GetOutput());
writer->SetInput(madFilter->GetOutput());
// Software Guide : EndCodeSnippet
try
{
// Software Guide : BeginLatex
//
// And then we can trigger the pipeline update, as usual.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
writer->Update();
// Software Guide : EndCodeSnippet
}
catch (itk::ExceptionObject& err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
// Here we generate the figures
InputPrintFilterType::Pointer input1PrintFilter = InputPrintFilterType::New();
InputPrintFilterType::Pointer input2PrintFilter = InputPrintFilterType::New();
OutputPrintFilterType::Pointer outputPrintFilter = OutputPrintFilterType::New();
VisuWriterType::Pointer input1VisuWriter = VisuWriterType::New();
VisuWriterType::Pointer input2VisuWriter = VisuWriterType::New();
VisuWriterType::Pointer outputVisuWriter = VisuWriterType::New();
input1PrintFilter->SetInput(reader1->GetOutput());
input1PrintFilter->SetChannel(3);
input1PrintFilter->SetChannel(2);
input1PrintFilter->SetChannel(1);
input2PrintFilter->SetInput(reader2->GetOutput());
input2PrintFilter->SetChannel(3);
input2PrintFilter->SetChannel(2);
input2PrintFilter->SetChannel(1);
outputPrintFilter->SetInput(madFilter->GetOutput());
outputPrintFilter->SetChannel(3);
outputPrintFilter->SetChannel(2);
outputPrintFilter->SetChannel(1);
input1VisuWriter->SetInput(input1PrintFilter->GetOutput());
input2VisuWriter->SetInput(input2PrintFilter->GetOutput());
outputVisuWriter->SetInput(outputPrintFilter->GetOutput());
input1VisuWriter->SetFileName(in1pretty);
input2VisuWriter->SetFileName(in2pretty);
outputVisuWriter->SetFileName(outpretty);
input1VisuWriter->Update();
input2VisuWriter->Update();
outputVisuWriter->Update();
// Software Guide : BeginLatex
// Figure \ref{fig:MADCHDET} shows the result of the change
// detection by Multivariate Alteration Detector.
// \begin{figure}
// \center
// \includegraphics[width=0.32\textwidth]{outIn1Pretty.eps}
// \includegraphics[width=0.32\textwidth]{outIn2Pretty.eps}
// \includegraphics[width=0.32\textwidth]{outPretty.eps}
// \itkcaption[CorrelationMultivariate Alteration Detection Results]{Result of the
// Multivariate Alteration Detector results on SPOT5 data before and
// after flooding.}
// \label{fig:RESCORRCHDET}
// \end{figure}
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
#include "otbImage.h"
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbStreamingImageFileWriter.h"
#include "otbMeanShiftVectorImageFilter.h"
//Code adapted from submission from Christophe Simler
// http://bugs.orfeo-toolbox.org/view.php?id=41
int main(int argc, char *argv[])
{
if (argc < 9)
{
std::cout <<
"Usage : inputImage rangeRadius spatialRadius minRegionSize outfilenamefiltered outfilenamesegmented outfilenamelabeled outfilenameboundary"
<< std::endl;
return EXIT_FAILURE;
}
char * filename = argv[1];
int rangeRadius = atoi(argv[2]);
int spatialRadius = atoi(argv[3]);
int minRegionSize = atoi(argv[4]);
char * outfilenamefiltered = argv[5];
char * outfilenamesegmented = argv[6];
char * outfilenamelabeled = argv[7];
char * outfilenameboundary = argv[8];
typedef otb::VectorImage<unsigned char, 2> ImageType; // image d'entree, image filtree et image segmente
typedef otb::Image<int, 2> TLabeledOutput; // image labelisee, image des contours (de l'image labellisee)
// lecture de l'image d'entree a partir d'un fichier
typedef otb::ImageFileReader<ImageType> ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(filename);
// traitement avec le filtre
typedef otb::MeanShiftVectorImageFilter<ImageType, ImageType, TLabeledOutput> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetRangeRadius(rangeRadius);
filter->SetSpatialRadius(spatialRadius);
filter->SetMinimumRegionSize(minRegionSize);
// sauvegarde de l'image filtree,
typedef otb::StreamingImageFileWriter<ImageType> WriterType1;
WriterType1::Pointer writer1 = WriterType1::New();
writer1->SetFileName(outfilenamefiltered);
// sauvegarde de l'image segmente,
typedef otb::StreamingImageFileWriter<ImageType> WriterType2;
WriterType2::Pointer writer2 = WriterType2::New();
writer2->SetFileName(outfilenamesegmented);
// sauvegarde de l'image labelisee
typedef otb::StreamingImageFileWriter<TLabeledOutput> WriterType3;
WriterType3::Pointer writer3 = WriterType3::New();
writer3->SetFileName(outfilenamelabeled);
// sauvegarde de l'image de contours
typedef otb::StreamingImageFileWriter<TLabeledOutput> WriterType4;
WriterType4::Pointer writer4 = WriterType4::New();
writer4->SetFileName(outfilenameboundary);
// construction du pipeline
filter->SetInput(reader->GetOutput());
writer1->SetInput(filter->GetOutput()); // image filtree (*)
writer2->SetInput(filter->GetClusteredOutput()); // image segmente (clusterisee avec coherence spaciale ?) (*)
writer3->SetInput(filter->GetLabeledClusteredOutput()); // image des labels des clusters
writer4->SetInput(filter->GetClusterBoundariesOutput()); // image des contours des clusters (contours de l'image labelisee)
writer1->Update();
writer2->Update();
writer3->Update();
writer4->Update();
return 0;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
//Code adapted from submission from Jordi INGLADA
// http://bugs.orfeo-toolbox.org/view.php?id=132
int main(int argc, char *argv[])
{
if (argc < 1)
{
std::cout << "Usage : inputImage" << std::endl;
return EXIT_FAILURE;
}
char * filename = argv[1];
typedef double PixelType;
typedef otb::VectorImage<PixelType> ImageType;
typedef otb::ImageFileReader<ImageType> ReaderType;
// check for input images
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(filename);
reader->UpdateOutputInformation();
std::cout << reader << std::endl;
return EXIT_SUCCESS;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "itkMacro.h"
#include "itkFFTComplexToComplexImageFilter.h"
int main(int argc, char * argv[])
{
typedef itk::FFTComplexToComplexImageFilter <float, 2> FFTFilterType;
FFTFilterType::Pointer fftFilter = FFTFilterType::New();
fftFilter->DebugOn();
return EXIT_SUCCESS;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "itkMacro.h"
#include "itkListSample.h"
#include <iostream>
#include "otbSVMSampleListModelEstimator.h"
#include "otbSVMClassifier.h"
#include "otbSVMKernels.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "otbSVMClassifier.h"
#include "otbConfusionMatrixCalculator.h"
#include "otbMath.h"
#include <fstream>
/*
This test show a problem with the SVM library using the probability
estimation : bug 209.
If the probability estimation is activated, the classifier isn't
abble to find the hyperplan even if the sample are linearly
seperable.
cf. test leTvBug209_SVMValidationLinearlySeparableWithoutProbEstimate
=> OK
and leTvBug209_SVMValidationLinearlySeparableWithProbEstimate => KO
http://bugs.orfeo-toolbox.org/view.php?id=209
*/
int main(int argc, char* argv[])
{
if(argc != 14)
{
std::cerr<<"Usage: "<<argv[0]<<" nbTrainingSamples nbValidationSamples positiveCenterX positiveCenterY negativeCenterX negativeCenterY positiveRadiusMin positiveRadiusMax negativeRadiusMin negativeRadiusMax kernel probEstimate"<<std::endl;
return EXIT_FAILURE;
}
unsigned int nbTrainingSamples = atoi(argv[2]);
unsigned int nbValidationSamples = atoi(argv[3]);
double cpx = atof(argv[4]);
double cpy = atof(argv[5]);
double cnx = atof(argv[6]);
double cny = atof(argv[7]);
double prmin = atof(argv[8]);
double prmax = atof(argv[9]);
double nrmin = atof(argv[10]);
double nrmax = atof(argv[11]);
unsigned int kernel = atoi(argv[12]);
bool probEstimate = atoi(argv[13]);
typedef double InputPixelType;
typedef unsigned short LabelType;
typedef itk::VariableLengthVector<InputPixelType> SampleType;
typedef itk::Statistics::ListSample<SampleType> ListSampleType;
typedef itk::FixedArray<LabelType, 1> TrainingSampleType;
typedef itk::Statistics::ListSample<TrainingSampleType> TrainingListSampleType;
typedef itk::Statistics::MersenneTwisterRandomVariateGenerator RandomGeneratorType;
typedef otb::SVMSampleListModelEstimator<ListSampleType, TrainingListSampleType> EstimatorType;
typedef otb::SVMClassifier<ListSampleType, LabelType> ClassifierType;
typedef ClassifierType::OutputType ClassifierOutputType;
typedef otb::ConfusionMatrixCalculator
<TrainingListSampleType, TrainingListSampleType> ConfusionMatrixCalculatorType;
RandomGeneratorType::Pointer random = RandomGeneratorType::New();
random->SetSeed((unsigned int)0);
// First, generate training and validation sets
ListSampleType::Pointer trainingSamples = ListSampleType::New();
TrainingListSampleType::Pointer trainingLabels = TrainingListSampleType::New();
ListSampleType::Pointer validationSamples = ListSampleType::New();
TrainingListSampleType::Pointer validationLabels = TrainingListSampleType::New();
// Generate training set
//std::ofstream training("training.csv");
for(unsigned int i =0; i < nbTrainingSamples; ++i)
{
// Generate a positive sample
double angle = random->GetVariateWithOpenUpperRange( otb::CONST_2PI );
double radius = random->GetUniformVariate(prmin, prmax);
SampleType pSample(2);
pSample[0] = cpx+radius*vcl_sin(angle);
pSample[1] = cpy+radius*vcl_cos(angle);
TrainingSampleType label;
label[0]=1;
trainingSamples->PushBack(pSample);
trainingLabels->PushBack(label);
//training<<"1 1:"<<pSample[0]<<" 2:"<<pSample[1]<<std::endl;
// Generate a negative sample
angle = random->GetVariateWithOpenUpperRange( otb::CONST_2PI );
radius = random->GetUniformVariate(nrmin, nrmax);
SampleType nSample(2);
nSample[0] = cnx+radius*vcl_sin(angle);
nSample[1] = cny+radius*vcl_cos(angle);
label[0]=2;
trainingSamples->PushBack(nSample);
trainingLabels->PushBack(label);
//training<<"2 1:"<<nSample[0]<<" 2:"<<nSample[1]<<std::endl;
}
//training.close();
// Generate validation set
std::ofstream validation("validation.csv");
for(unsigned int i =0; i < nbValidationSamples; ++i)
{
// Generate a positive sample
double angle = random->GetVariateWithOpenUpperRange( otb::CONST_2PI );
double radius = random->GetUniformVariate(prmin, prmax);
SampleType pSample(2);
pSample[0] = cpx+radius*vcl_sin(angle);
pSample[1] = cpy+radius*vcl_cos(angle);
TrainingSampleType label;
label[0]=1;
validationSamples->PushBack(pSample);
validationLabels->PushBack(label);
//validation<<"1 1:"<<pSample[0]<<" 2:"<<pSample[1]<<std::endl;
// Generate a negative sample
angle = random->GetVariateWithOpenUpperRange( otb::CONST_2PI );
radius = random->GetUniformVariate(nrmin, nrmax);
SampleType nSample(2);
nSample[0] = cnx+radius*vcl_sin(angle);
nSample[1] = cny+radius*vcl_cos(angle);
label[0]=2;
validationSamples->PushBack(nSample);
validationLabels->PushBack(label);
//validation<<"2 1:"<<nSample[0]<<" 2:"<<nSample[1]<<std::endl;
}
//validation.close();
// Learn
EstimatorType::Pointer estimator = EstimatorType::New();
estimator->SetInputSampleList(trainingSamples);
estimator->SetTrainingSampleList(trainingLabels);
estimator->SetKernelType(kernel);
estimator->DoProbabilityEstimates(probEstimate);
// estimator->SetParametersOptimization(true);
estimator->Update();
//estimator->SaveModel("model.svm");
// Classify
ClassifierType::Pointer validationClassifier = ClassifierType::New();
validationClassifier->SetSample(validationSamples);
validationClassifier->SetNumberOfClasses(2);
validationClassifier->SetModel(estimator->GetModel());
validationClassifier->Update();
// Confusion
ClassifierOutputType::ConstIterator it = validationClassifier->GetOutput()->Begin();
ClassifierOutputType::ConstIterator itEnd = validationClassifier->GetOutput()->End();
TrainingListSampleType::Pointer classifierListLabel = TrainingListSampleType::New();
while (it != itEnd)
{
classifierListLabel->PushBack(it.GetClassLabel());
++it;
}
ConfusionMatrixCalculatorType::Pointer confMatCalc = ConfusionMatrixCalculatorType::New();
confMatCalc->SetReferenceLabels(validationLabels);
confMatCalc->SetProducedLabels(classifierListLabel);
confMatCalc->Update();
std::cout<<std::endl;
std::cout<<"Confusion matrix: "<<std::endl<< confMatCalc->GetConfusionMatrix()<<std::endl<<std::endl;
std::cout<<"Kappa Index: "<<std::endl<< confMatCalc->GetKappaIndex()<<std::endl<<std::endl;
if(confMatCalc->GetKappaIndex()!=1)
{
std::cerr<<"Kappa index should be 1."<<std::endl;
return EXIT_FAILURE;
}
else
{
return EXIT_SUCCESS;
}
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbStreamingImageFileWriter.h"
#include "itkCastImageFilter.h"
/**
* This test reproduces a problem encountered with CastImageFilter and a streaming
* writer. When the image to write is already loaded in memory and InPlaceOn() is
* called, only the first stripe is repeated across the output image.
* Issue number 0000428
*/
int main(int argc, char* argv[])
{
if (argc != 3)
{
std::cout << argv[0] << " <input image> <output image>" << std::endl;
return EXIT_FAILURE;
}
typedef otb::Image<float, 2> ImageType;
typedef otb::ImageFileReader<ImageType> ReaderType;
typedef otb::StreamingImageFileWriter<ImageType> WriterType;
typedef itk::CastImageFilter<ImageType, ImageType> CastFilterType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(argv[1]);
reader->Update();
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput( reader->GetOutput() );
caster->InPlaceOn();
WriterType::Pointer writer = WriterType::New();
writer->SetInput(caster->GetOutput());
writer->SetFileName(argv[2]);
//writer->SetAutomaticTiledStreaming(1024);
writer->SetNumberOfDivisionsStrippedStreaming(10);
writer->Update();
return EXIT_SUCCESS;
}
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "otbImage.h"
#include "otbLineSegmentDetector.h"
#include "otbImageFileReader.h"
#include "otbVectorDataFileWriter.h"
#include "otbImageFileWriter.h"