diff --git a/Examples/FeatureExtraction/CMakeLists.txt b/Examples/FeatureExtraction/CMakeLists.txt index 515192f052d6bdbfc02fde009e8aeb7d638a1006..9f9c3ee48bb4ce6983ca5ea7531ffa6e24add36b 100644 --- a/Examples/FeatureExtraction/CMakeLists.txt +++ b/Examples/FeatureExtraction/CMakeLists.txt @@ -47,6 +47,11 @@ ADD_EXECUTABLE(ExtractSegmentsExample ExtractSegmentsExample.cxx) TARGET_LINK_LIBRARIES(ExtractSegmentsExample OTBFeatureExtraction OTBCommon OTBIO ITKCommon ITKIO) +ADD_EXECUTABLE(ExtractRoadByStepsExample ExtractRoadByStepsExample.cxx) +TARGET_LINK_LIBRARIES(ExtractRoadByStepsExample OTBIO OTBCommon OTBFeatureExtraction +ITKCommon ITKBasicFilters) + + IF( NOT OTB_DISABLE_CXX_TESTING AND NOT OTB_DISABLE_CXX_EXAMPLES_TESTING ) SET(BASELINE ${OTB_DATA_ROOT}/Baseline/Examples/FeatureExtraction) @@ -176,6 +181,17 @@ ADD_TEST(AlignmentsExampleTest ${EXE_TESTS} ) +ADD_TEST(ExtractRoadByStepsExampleTest ${EXE_TESTS} + --compare-image ${TOL} + ${BASELINE}/ + ${TEMP}/ + ExtractRoadByStepsExampleTest + ${INPUTDATA}/qb_RoadExtract.img.hdr + ${TEMP}/ExtractRoadByStepsExampleOutput.png + 337 557 432 859 0.00005 1.0 +) + + INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}") ADD_EXECUTABLE(otbFeatureExtractionExamplesTests otbFeatureExtractionExamplesTests.cxx) diff --git a/Examples/FeatureExtraction/ExtractRoadByStepsExample.cxx b/Examples/FeatureExtraction/ExtractRoadByStepsExample.cxx new file mode 100755 index 0000000000000000000000000000000000000000..1cd8c7caa53e696f1805b4d799d79130a42e3e32 --- /dev/null +++ b/Examples/FeatureExtraction/ExtractRoadByStepsExample.cxx @@ -0,0 +1,352 @@ +/*========================================================================= + + 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. + +=========================================================================*/ + +#if defined(_MSC_VER) +#pragma warning ( disable : 4786 ) +#endif + +#ifdef __BORLANDC__ +#define ITK_LEAN_AND_MEAN +#endif + +// Software Guide : BeginCommandLineArgs +// INPUTS: {qb_extract.tif} +// OUTPUTS: {qb_extractSegmentExtractionBySteps.png} +// 320 506 367 790 +// Software Guide : EndCommandLineArgs + +// Software Guide : BeginLatex +// +// This example illustrates the detail of the \doxygen{otb}{RoadExtractionFilter}. +// This filter is a composite filter including all the steps below. Individual +// filters can be replaced to design a road detector targeted at SAR images for +// example. +// +// The first step required to use this filter is to include header files. +// +// Software Guide : EndLatex + +// Software Guide : BeginCodeSnippet + +#include "otbSpectralAngleDistanceImageFilter.h" +#include "itkGradientRecursiveGaussianImageFilter.h" +#include "otbNeighborhoodScalarProductFilter.h" +#include "otbRemoveIsolatedByDirectionFilter.h" +#include "otbRemoveWrongDirectionFilter.h" +#include "otbNonMaxRemovalByDirectionFilter.h" +#include "otbVectorizationPathListFilter.h" +#include "otbSimplifyPathListFilter.h" +#include "otbBreakAngularPathListFilter.h" +#include "otbRemoveTortuousPathListFilter.h" +#include "otbLinkPathListFilter.h" +#include "otbLikehoodPathListFilter.h" +#include "otbDrawPathListFilter.h" +#include "otbLikehoodPathListFilter.h" +#include "otbMultiToMonoChannelExtractROI.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkSqrtImageFilter.h" + +// Software Guide : EndCodeSnippet + +#include "itkCovariantVector.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbPolyLineParametricPathWithValue.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" + + +// Software Guide : BeginCommandLineArgs +// INPUTS: {qb_RoadExtract.img.hdr} +// OUTPUTS: {ExtractRoadByStepsExampleOutput.png} +// 337 557 432 859 0.00005 1.0 +// Software Guide : EndCommandLineArgs + +int main( int argc, char * argv[] ) +{ + + const unsigned int Dimension = 2; + typedef double PixelType; + typedef itk::CovariantVector<PixelType,Dimension> VectorPixelType; + typedef otb::Image<PixelType, Dimension> InternalImageType; + typedef otb::VectorImage<PixelType, Dimension> MultiSpectralImageType; + typedef otb::Image<unsigned char, Dimension> OutputImageType; + + typedef otb::Image<VectorPixelType, Dimension > VectorImageType; + + typedef otb::PolyLineParametricPathWithValue< double, Dimension > PathType; + + + typedef otb::ImageFileReader< MultiSpectralImageType > MultispectralReaderType; + MultispectralReaderType::Pointer multispectralReader = MultispectralReaderType::New(); + multispectralReader->SetFileName(argv[1]); + + typedef otb::ImageFileWriter< OutputImageType > FileWriterType; + FileWriterType::Pointer writer = FileWriterType::New(); + + // NB: There might be a better way to pass this parameter (coordinate of the reference ?) + // plan combination with the viewer + // possibility to give 2 parameters (just in future use) + MultiSpectralImageType::PixelType pixelRef; + pixelRef.SetSize(4); + pixelRef[0]=atoi(argv[3]); + pixelRef[1]=atoi(argv[4]); + pixelRef[2]=atoi(argv[5]); + pixelRef[3]=atoi(argv[6]); + + double resolution = 0.6; //to get directly from metadata ? + double alpha = atof(argv[8]); + // Software Guide : BeginLatex + // + // The spectral angle is used to compute a grayscale images from the + // multispectral original image. Pixels corresponding to roads are in + // darker color. + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + typedef otb::SpectralAngleDistanceImageFilter<MultiSpectralImageType,InternalImageType> SAFilterType; + SAFilterType::Pointer saFilter = SAFilterType::New(); + saFilter->SetReferencePixel(pixelRef); + saFilter->SetInput(multispectralReader->GetOutput()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // A square root is applied to the spectral angle image in order to enhance contrast between + // darker pixels (which are pixels of interest). + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + typedef itk::SqrtImageFilter<InternalImageType,InternalImageType> SqrtFilterType; + SqrtFilterType::Pointer sqrtFilter = SqrtFilterType::New(); + sqrtFilter->SetInput(saFilter->GetOutput()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // Use the Gaussian gradient filter + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + double sigma = alpha*(1.2/resolution+1); + typedef itk::GradientRecursiveGaussianImageFilter< InternalImageType, VectorImageType > GradientFilterType; + GradientFilterType::Pointer gradientFilter = GradientFilterType::New(); + gradientFilter->SetSigma(sigma); + gradientFilter->SetInput(sqrtFilter->GetOutput()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // Compute the scalar product of the neighboring pixels and keep the + // minimum value and the direction. This is the line detector described + // in \cite{Lacroix1998}. + // + // Software Guide : EndLatex + + //**** NB: One or two output ?**** + // -> two output, the direction image type can be different + // make sure OTB convention are respected, especially for the direction notation + + + // input output outputdir + // Software Guide : BeginCodeSnippet + typedef otb::NeighborhoodScalarProductFilter<VectorImageType,InternalImageType,InternalImageType> NeighborhoodScalarProductType; + NeighborhoodScalarProductType::Pointer scalarFilter = NeighborhoodScalarProductType::New(); + scalarFilter->SetInput(gradientFilter->GetOutput()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // The resulting image is passed to a filter to remove the pixel + // with no neighbor having the same direction. + // + // Software Guide : EndLatex + + // input inputdir output + // Software Guide : BeginCodeSnippet + typedef otb::RemoveIsolatedByDirectionFilter<InternalImageType,InternalImageType,InternalImageType> RemoveIsolatedByDirectionType; + RemoveIsolatedByDirectionType::Pointer removeIsolatedByDirectionFilter = RemoveIsolatedByDirectionType::New(); + removeIsolatedByDirectionFilter->SetInput(scalarFilter->GetOutput()); + removeIsolatedByDirectionFilter->SetInputDirection(scalarFilter->GetOutputDirection()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // We now remove pixels having a direction corresponding to bright lines + // as we know that after the spectral angle, roads are in darker color. + // + // Software Guide : EndLatex + + // NB: Name can be change according to the conventional values (RemoveNegative ? Positive ?) + // input inputdir output + // Software Guide : BeginCodeSnippet + typedef otb::RemoveWrongDirectionFilter<InternalImageType,InternalImageType,InternalImageType> RemoveWrongDirectionType; + RemoveWrongDirectionType::Pointer removeWrongDirectionFilter = RemoveWrongDirectionType::New(); + removeWrongDirectionFilter->SetInput(removeIsolatedByDirectionFilter->GetOutput()); + removeWrongDirectionFilter->SetInputDirection(scalarFilter->GetOutputDirection()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // We now remove pixels which are not maximum on the direction + // perpendicular to the road direction. + // + // Software Guide : EndLatex + + // input inputdir output + // Software Guide : BeginCodeSnippet + typedef otb::NonMaxRemovalByDirectionFilter<InternalImageType,InternalImageType,InternalImageType> NonMaxRemovalByDirectionType; + NonMaxRemovalByDirectionType::Pointer nonMaxRemovalByDirectionFilter = NonMaxRemovalByDirectionType::New(); + nonMaxRemovalByDirectionFilter->SetInput(removeWrongDirectionFilter->GetOutput()); + nonMaxRemovalByDirectionFilter->SetInputDirection(scalarFilter->GetOutputDirection()); + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // Extracted road are vectorized into polylines + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + typedef otb::VectorizationPathListFilter<InternalImageType,InternalImageType,PathType> VectorizationFilterType; + VectorizationFilterType::Pointer vectorizationFilter = VectorizationFilterType::New(); + vectorizationFilter->SetInput(nonMaxRemovalByDirectionFilter->GetOutput()); + vectorizationFilter->SetInputDirection(scalarFilter->GetOutputDirection()); + vectorizationFilter->SetAmplitudeThreshold(atof(argv[7])); + // Software Guide : EndCodeSnippet + + + // Software Guide : BeginLatex + // + // However, this vectorization is too simple and need to be refined + // to be usable. First, we remove all aligned points to make one segment. + // Then we break the polylines which have sharp angles as they are probably + // not road. Finally we remove path which are too short. + // + // Software Guide : EndLatex + + //NB: This filter works on one path at a time. Possibly better to work on PathType individually ? + // but we must have a simple way to simplify the full list + // Software Guide : BeginCodeSnippet + typedef otb::SimplifyPathListFilter<PathType> SimplifyPathType; + SimplifyPathType::Pointer simplifyPathListFilter = SimplifyPathType::New(); + simplifyPathListFilter->SetTolerance(1.0); + simplifyPathListFilter->SetInput(vectorizationFilter->GetOutput()); + + //NB: This filter need the full list (adding the new created path at the end). + typedef otb::BreakAngularPathListFilter<PathType> BreakAngularPathType; + BreakAngularPathType::Pointer breakAngularPathListFilter = BreakAngularPathType::New(); + breakAngularPathListFilter->SetMaxAngle(M_PI/8.); + breakAngularPathListFilter->SetInput(simplifyPathListFilter->GetOutput()); + + + typedef otb::RemoveTortuousPathListFilter<PathType> RemoveTortuousPathType; + RemoveTortuousPathType::Pointer removeTortuousPathListFilter = RemoveTortuousPathType::New(); + removeTortuousPathListFilter->SetMeanDistanceThreshold(1.0); + removeTortuousPathListFilter->SetInput(breakAngularPathListFilter->GetOutput()); + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // Polylines within a certain range are linked to try to fill gaps due to + // occultations by vehicules, trees, etc. before simplifying polylines and + // removing the shorteest ones. + // + // Software Guide : EndLatex + + + // Software Guide : BeginCodeSnippet + typedef otb::LinkPathListFilter<PathType> LinkPathType; + LinkPathType::Pointer linkPathListFilter = LinkPathType::New(); + linkPathListFilter->SetDistanceThreshold(25.0/resolution);//research area of 25 m + linkPathListFilter->SetAngularThreshold(M_PI/8); + linkPathListFilter->SetInput(removeTortuousPathListFilter->GetOutput()); + + SimplifyPathType::Pointer simplifyPathListFilter2 = SimplifyPathType::New(); + simplifyPathListFilter2->SetTolerance(1.0); + simplifyPathListFilter2->SetInput(linkPathListFilter->GetOutput()); + + RemoveTortuousPathType::Pointer removeTortuousPathListFilter2 = RemoveTortuousPathType::New(); + removeTortuousPathListFilter2->SetMeanDistanceThreshold(10.0); + removeTortuousPathListFilter2->SetInput(simplifyPathListFilter2->GetOutput()); + // Software Guide : EndCodeSnippet + + + // Software Guide : BeginLatex + // + // A value can be associated with each polyline according to pixel values + // under the polyline. A higher value will mean a higher likelihood to be + // a road. Polylines are drawn on the original image. + // + // Software Guide : EndLatex + + //NB: check if PathListWithValue is correct according to conventions + // Software Guide : BeginCodeSnippet + typedef otb::LikehoodPathListFilter<PathType,InternalImageType> PathListToPathListWithValueType; + PathListToPathListWithValueType::Pointer pathListConverter = PathListToPathListWithValueType::New(); + pathListConverter->SetInput(removeTortuousPathListFilter2->GetOutput()); + pathListConverter->SetInputImage(nonMaxRemovalByDirectionFilter->GetOutput()); + + typedef otb::MultiToMonoChannelExtractROI<PixelType,PixelType> ChannelExtractionFilterType; + ChannelExtractionFilterType::Pointer channelExtractor = ChannelExtractionFilterType::New(); + channelExtractor->SetInput(multispectralReader->GetOutput()); + + typedef itk::RescaleIntensityImageFilter<InternalImageType,InternalImageType> RescalerType; + RescalerType::Pointer rescaler = RescalerType::New(); + rescaler->SetOutputMaximum(255); + rescaler->SetOutputMinimum(0); + rescaler->SetInput(channelExtractor->GetOutput()); + + typedef otb::DrawPathListFilter<InternalImageType, PathType, OutputImageType> DrawPathType; + DrawPathType::Pointer drawPathListFilter = DrawPathType::New(); + drawPathListFilter->SetInput(rescaler->GetOutput()); + drawPathListFilter->SetInputPath(pathListConverter->GetOutput()); + drawPathListFilter->SetPathValue(255); + // Software Guide : EndCodeSnippet + + + // Software Guide : BeginLatex + // + // The filter is executed by invoking the \code{Update()} method. If the + // filter is part of a larger image processing pipeline, calling + // \code{Update()} on a downstream filter will also trigger update of this + // filter. + // + // Extracted roads are drawn on the original image and output to a file. + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + writer->SetInput( drawPathListFilter->GetOutput() ); + writer->SetFileName( argv[2] ); + writer->Update(); + // Software Guide : EndCodeSnippet + + return EXIT_SUCCESS; +} + diff --git a/Examples/FeatureExtraction/otbFeatureExtractionExamplesTests.cxx b/Examples/FeatureExtraction/otbFeatureExtractionExamplesTests.cxx index 5c7f89f677fee36333523a786147948f8505d8e8..6db160ead1c32e6590fb0bd8de08ea8abf3f0bb9 100644 --- a/Examples/FeatureExtraction/otbFeatureExtractionExamplesTests.cxx +++ b/Examples/FeatureExtraction/otbFeatureExtractionExamplesTests.cxx @@ -34,6 +34,7 @@ REGISTER_TEST(AssymmetricFusionOfLineDetectorExampleTest); REGISTER_TEST(ExtractSegmentsExampleTest); REGISTER_TEST(RatioLineDetectorExampleTest); REGISTER_TEST(AlignmentsExampleTest); +REGISTER_TEST(ExtractRoadByStepsExampleTest); } #undef main @@ -72,3 +73,7 @@ REGISTER_TEST(AlignmentsExampleTest); #define main AlignmentsExampleTest #include "AlignmentsExample.cxx" +#undef main +#define main ExtractRoadsByStepsExampleTest +#include "ExtractRoadByStepsExample.cxx" + diff --git a/Testing/Code/IO/CMakeLists.txt b/Testing/Code/IO/CMakeLists.txt index 6d333a4acc7906cbb40e72d307be302606bf20a7..5dac56fbf4c6702caadc3f719ec010ccb0b7650b 100755 --- a/Testing/Code/IO/CMakeLists.txt +++ b/Testing/Code/IO/CMakeLists.txt @@ -833,6 +833,17 @@ ADD_TEST(ioTuOtbVectorImageTestSpot5 ${IO_TESTS} +# --- otb::VectorImageReadWriteTest --- +ADD_TEST(ioTvVectorImageFileReaderWriterTest ${IO_TESTS} +# --compare-image ${EPS} +# ${BASELINE_FILES}/ +# ${TEMP}/ + otbVectorImageFileReaderWriterTest + ${IMAGEDATA}/qb_RoadExtract.img.hdr + ${TEMP}/ioTvVectorImageFileReaderWriterTest.hdr +) + + # --- STREAMING : GDAL --- ADD_TEST(ioTvStreamingImageFilterGDAL_ENVI ${IO_TESTS} # --compare-image ${TOL} ${TEMP}/ioStreamingImageFilterTestCDAL_ENVI.png