ExtractRoadExample.cxxΒΆ
Example usage:
./ExtractRoadExample Input/qb_RoadExtract.tif Output/ExtractRoadOutput.png 337 557 432 859 1.0 0.00005 1.0 0.39269 1.0 10.0 25.
Example source code (ExtractRoadExample.cxx):
// The easiest way to use the road extraction filter provided by OTB is to use the composite
// filter. If a modification in the pipeline is required to adapt to a particular situation,
// the step by step example, described in the next section can be adapted.
//
// This example demonstrates the use of the \doxygen{otb}{RoadExtractionFilter}.
// This filter is a composite filter achieving road extraction according to the algorithm
// adapted by E. Christophe and J. Inglada \cite{Christophe2007} from an original method
// proposed in \cite{Lacroix1998}.
//
// The first step toward the use of this filter is the inclusion of the proper header files.
#include "otbPolyLineParametricPathWithValue.h"
#include "otbRoadExtractionFilter.h"
#include "otbDrawPathListFilter.h"
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "otbMath.h"
#include "itkInvertIntensityImageFilter.h"
#include "itkGrayscaleDilateImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
int main(int argc, char* argv[])
{
if (argc != 14)
{
std::cerr << "Usage: " << argv[0];
std::cerr << " inputFileName outputFileName firstPixelComponent secondPixelComponent ";
std::cerr << "thirdPixelComponent fourthPixelComponent alpha amplitudeThrehsold tolerance ";
std::cerr << "angularThreshold-maxAngle firstMeanDistanceThreshold secondMeanDistanceThreshold ";
std::cerr << "distanceThreshold" << std::endl;
return EXIT_FAILURE;
}
const unsigned int Dimension = 2;
// Then we must decide what pixel type to use for the image. We choose to do
// all the computation in floating point precision and rescale the results
// between 0 and 255 in order to export PNG images.
typedef double InputPixelType;
typedef unsigned char OutputPixelType;
// The images are defined using the pixel type and the dimension. Please note that
// the \doxygen{otb}{RoadExtractionFilter} needs an \doxygen{otb}{VectorImage} as input
// to handle multispectral images.
typedef otb::VectorImage<InputPixelType, Dimension> InputVectorImageType;
typedef otb::Image<InputPixelType, Dimension> InputImageType;
typedef otb::Image<OutputPixelType, Dimension> OutputImageType;
// We define the type of the polyline that the filter produces. We use the
// \doxygen{otb}{PolyLineParametricPathWithValue}, which allows the filter to produce
// a likehood value along with each polyline. The filter is able to produce
// \doxygen{itk}{PolyLineParametricPath} as well.
typedef otb::PolyLineParametricPathWithValue<InputPixelType, Dimension> PathType;
// Now we can define the \doxygen{otb}{RoadExtractionFilter} that takes a multi-spectral
// image as input and produces a list of polylines.
typedef otb::RoadExtractionFilter<InputVectorImageType, PathType> RoadExtractionFilterType;
// We also define an \doxygen{otb}{DrawPathListFilter} to draw the output
// polylines on an image, taking their likehood values into account.
typedef otb::DrawPathListFilter<InputImageType, PathType, InputImageType> DrawPathFilterType;
// The intensity rescaling of the results will be carried out by the
// \doxygen{itk}{RescaleIntensityImageFilter} which is templated by the
// input and output image types.
typedef itk::RescaleIntensityImageFilter<InputImageType, OutputImageType> RescalerType;
// An \doxygen{otb}{ImageFileReader} class is also instantiated in order to read
// image data from a file. Then, an \doxygen{otb}{ImageFileWriter}
// is instantiated in order to write the output image to a file.
typedef otb::ImageFileReader<InputVectorImageType> ReaderType;
typedef otb::ImageFileWriter<OutputImageType> WriterType;
// The different filters composing our pipeline are created by invoking their
// \code{New()} methods, assigning the results to smart pointers.
ReaderType::Pointer reader = ReaderType::New();
RoadExtractionFilterType::Pointer roadExtractionFilter = RoadExtractionFilterType::New();
DrawPathFilterType::Pointer drawingFilter = DrawPathFilterType::New();
RescalerType::Pointer rescaleFilter = RescalerType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName(argv[1]);
// The \doxygen{otb}{RoadExtractionFilter} needs to have a reference pixel
// corresponding to the spectral content likely to represent a road. This is done
// by passing a pixel to the filter. Here we suppose that the input image
// has four spectral bands.
InputVectorImageType::PixelType ReferencePixel;
ReferencePixel.SetSize(4);
ReferencePixel.SetElement(0, ::atof(argv[3]));
ReferencePixel.SetElement(1, ::atof(argv[4]));
ReferencePixel.SetElement(2, ::atof(argv[5]));
ReferencePixel.SetElement(3, ::atof(argv[6]));
roadExtractionFilter->SetReferencePixel(ReferencePixel);
// We must also set the alpha parameter of the filter which allows us to tune the width of the roads
// we want to extract. Typical value is $1.0$ and should be working in most situations.
roadExtractionFilter->SetAlpha(atof(argv[7]));
// All other parameter should not influence the results too much in most situation and can
// be kept at the default value.
//
// The amplitude threshold parameter tunes the sensitivity of the vectorization step. A typical
// value is $5 \cdot 10^{-5}$.
roadExtractionFilter->SetAmplitudeThreshold(atof(argv[8]));
// The tolerance threshold tunes the sensitivity of the path simplification step.
// Typical value is $1.0$.
roadExtractionFilter->SetTolerance(atof(argv[9]));
// Roads are not likely to have sharp turns. Therefore we set the max angle parameter,
// as well as the link angular threshold. The value is typically $\frac{\pi}{8}$.
roadExtractionFilter->SetMaxAngle(atof(argv[10]));
roadExtractionFilter->SetAngularThreshold(atof(argv[10]));
// The \doxygen{otb}{RoadExtractionFilter} performs two odd path removing operations at different stage of
// its execution. The first mean distance threshold and the second mean distance threshold set their criterion
// for removal. Path are removed if their mean distance between nodes is to small, since such path coming
// from previous filters are likely to be tortuous. The first removal operation as a typical mean distance
// threshold parameter of $1.0$, and the second of $10.0$.
roadExtractionFilter->SetFirstMeanDistanceThreshold(atof(argv[11]));
roadExtractionFilter->SetSecondMeanDistanceThreshold(atof(argv[12]));
// The \doxygen{otb}{RoadExtractionFilter} is able to link path whose ends are near
// according to an euclidean distance criterion. The threshold for this distance
// to link a path is the distance threshold parameter. A typical value is $25$.
roadExtractionFilter->SetDistanceThreshold(atof(argv[13]));
// We will now create a black background image to draw the resulting polyline on.
// To achieve this we need to know the size of our input image. Therefore we trigger the
// \code{GenerateOutputInformation()} of the reader.
reader->GenerateOutputInformation();
InputImageType::Pointer blackBackground = InputImageType::New();
blackBackground->CopyInformation(reader->GetOutput());
blackBackground->SetRegions(blackBackground->GetLargestPossibleRegion());
blackBackground->Allocate();
blackBackground->FillBuffer(0);
// We tell the \doxygen{otb}{DrawPathListFilter} to try to use the likehood value
// embedded within the polyline as a value for drawing this polyline if possible.
drawingFilter->UseInternalPathValueOn();
// The \code{itk::RescaleIntensityImageFilter} needs to know which
// is the minimum and maximum values of the output generated
// image. Those can be chosen in a generic way by using the
// \code{NumericTraits} functions, since they are templated over
// the pixel type.
rescaleFilter->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
rescaleFilter->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
// Now it is time for some pipeline wiring.
roadExtractionFilter->SetInput(reader->GetOutput());
drawingFilter->SetInput(blackBackground);
drawingFilter->SetInputPath(roadExtractionFilter->GetOutput());
rescaleFilter->SetInput(drawingFilter->GetOutput());
// The update of the pipeline is triggered by the \code{Update()} method
// of the rescale intensity filter.
rescaleFilter->Update();
// output image enhancement
typedef itk::BinaryBallStructuringElement<OutputPixelType, Dimension> StructuringElementType;
typedef itk::GrayscaleDilateImageFilter<OutputImageType, OutputImageType, StructuringElementType> DilateFilterType;
typedef itk::InvertIntensityImageFilter<OutputImageType, OutputImageType> InvertFilterType;
StructuringElementType se;
se.SetRadius(1);
se.CreateStructuringElement();
DilateFilterType::Pointer dilater = DilateFilterType::New();
dilater->SetInput(rescaleFilter->GetOutput());
dilater->SetKernel(se);
InvertFilterType::Pointer invertFilter = InvertFilterType::New();
invertFilter->SetInput(dilater->GetOutput());
writer->SetFileName(argv[2]);
writer->SetInput(invertFilter->GetOutput());
writer->Update();
// Figure~\ref{fig:ROADEXTRACTION_FILTER} shows the result of applying
// the road extraction filter to a fusionned Quickbird image.
// \begin{figure}
// \center
// \includegraphics[width=0.44\textwidth]{qb_ExtractRoad_pretty.eps}
// \includegraphics[width=0.44\textwidth]{ExtractRoadOutput.eps}
// \itkcaption[Road extraction filter application]{Result of applying
// the \doxygen{otb}{RoadExtractionFilter} to a fusionned Quickbird
// image. From left to right : original image, extracted road with their
// likehood values (color are inverted for display).}
// \label{fig:ROADEXTRACTION_FILTER}
// \end{figure}
return EXIT_SUCCESS;
}