An error occurred while loading the file. Please try again.
-
Guillaume Perréal authored8728b8a5
/*========================================================================= Copyright (c) IRSTEA. All rights reserved. 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 "otbDeltaNDVIFunctor.h"#include "itkFixedArray.h"#include "itkObjectFactory.h"// Elevation handler#include "otbWrapperElevationParametersHandler.h"#include "otbWrapperApplicationFactory.h"// Application engine#include "otbStandardFilterWatcher.h"#include "itkFixedArray.h"// Filters#include "itkMaskImageFilter.h"#include "otbStreamingStatisticsImageFilter.h"#include "otbStreamingResampleImageFilter.h"#include "otbMultiChannelExtractROI.h"#include "otbDeltaNDVILabelerFilter.h"#include "itkAndImageFilter.h"// Helper#include "otbRegionComparator.h"// Mask handler#include "otbMosaicFromDirectoryHandler.h"namespace otb{namespace Wrapper{class ClearCutsDetection : public Application{public: /** Standard class typedefs. */ typedef ClearCutsDetection Self; typedef Application Superclass; typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<const Self> ConstPointer; /** Standard macro */ itkNewMacro(Self); itkTypeMacro(ClearCutsDetection, Application); /** Filters */ typedef otb::StreamingResampleImageFilter<FloatVectorImageType, FloatVectorImageType> ResampleImageFilterType; typedef otb::StreamingResampleImageFilter<UInt8ImageType, UInt8ImageType> ResampleMaskImageFilterType; typedef itk::NearestNeighborInterpolateImageFunction< FloatVectorImageType> NNInterpolatorType; typedef itk::NearestNeighborInterpolateImageFunction< UInt8ImageType> NNMaskImageInterpolatorType; typedef otb::Functor::DeltaNDVIFromChannels<FloatVectorImageType::PixelType, FloatImageType::PixelType> DeltaNDVIFunctorType; typedef itk::BinaryFunctorImageFilter<FloatVectorImageType, FloatVectorImageType, FloatImageType, DeltaNDVIFunctorType> DeltaNDVIFilterType; typedef itk::MaskImageFilter<FloatImageType, UInt8ImageType, FloatImageType> MaskImageFilterType;7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
typedef itk::MaskImageFilter<FloatVectorImageType, UInt8ImageType,
FloatVectorImageType> MaskVectorImageFilterType;
typedef otb::DeltaNDVILabelerFilter<FloatImageType,
FloatImageType> NDVILabelImageFilterType;
typedef otb::StreamingStatisticsImageFilter<FloatImageType> StatsFilterType;
typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
FloatVectorImageType::InternalPixelType> ExtractROIFilterType;
typedef otb::MosaicFromDirectoryHandler<UInt8ImageType> MaskHandlerType;
typedef itk::AndImageFilter<UInt8ImageType, UInt8ImageType> AndFilterType;
void DoUpdateParameters()
{
// Nothing to do here : all parameters are independent
}
void DoInit()
{
SetName("ClearCutsDetection");
SetDescription("This application performs cuts detection, "
"from two input images and an optional forest mask");
// Documentation
SetDocName("ClearCutsDetection");
SetDocLimitations("None");
SetDocAuthors("RemiCresson");
SetDocLongDescription(" This filter implements the clear cut detection method, based "
"on the work of Kenji Ose and Michel Deshayes at IRSTEA. "
"Steps of the process are the following: 1) Compute the difference "
"between NDVI of dates t1 and t0. 2) compute mean and std of the "
"dNDVI on masked pixels. 3) Apply multiple thresholds based on "
"computed stats on the dNDVI image to obtain a label image which "
"quantize NDVI decrease. Input images must be in the same projection "
"but might have different origin/size (A NN-interpolation is done). "
"Outputs are dNDVI (out) and classes labels (outlabel). Classes "
"are the following: 0 is no-data, 1 is no vegetation decrease, "
"2 is little vegetation decrease, 3 is medium vegetation decrease, "
"and 4 is high vegetation decrease, that is, clear cuts.");
AddDocTag(Tags::FeatureExtraction);
// Input images
AddParameter(ParameterType_InputImage, "inb", "Input T0 Image (Before)");
AddParameter(ParameterType_InputImage, "ina", "Input T1 Image (After)");
// Input masks
AddParameter(ParameterType_InputImage, "inbmask", "Input T0 Image mask (Before)");
MandatoryOff("inbmask");
AddParameter(ParameterType_InputImage, "inamask", "Input T1 Image mask (After)");
MandatoryOff("inamask");
AddParameter(ParameterType_InputImage, "mask", "Input mask (used for both input images)");
MandatoryOff("mask");
AddParameter(ParameterType_Directory, "masksdir", "Vegetation masks directory");
MandatoryOff("masksdir");
// Input images band indices
AddParameter(ParameterType_Int, "nirb", "near infrared band index for input T0 image" );
SetParameterDescription("nirb","index of near infrared band of image b");
SetMinimumParameterIntValue("nirb", 1);
SetDefaultParameterInt("nirb",4);
AddParameter(ParameterType_Int, "redb", "red band index for input T0 image" );
SetParameterDescription("redb","index of red band of image b");
SetMinimumParameterIntValue("redb", 1);
SetDefaultParameterInt("redb",1);
AddParameter(ParameterType_Int, "nira", "near infrared band index for input T1 image" );
SetParameterDescription("nira","index of near infrared band of image a");
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
SetMinimumParameterIntValue("nira", 1);
SetDefaultParameterInt("nira",4);
AddParameter(ParameterType_Int, "reda", "red band index for input T1 image" );
SetParameterDescription("reda","index of red band of image a");
SetMinimumParameterIntValue("reda", 1);
SetDefaultParameterInt("reda",1);
// Output image
AddParameter(ParameterType_OutputImage, "outdndvi", "Output d(NDVI) image");
MandatoryOff("outdndvi");
SetDefaultOutputPixelType("outdndvi",ImagePixelType_float);
AddParameter(ParameterType_OutputImage, "outlabel", "Output labeled image");
SetDefaultOutputPixelType("outlabel",ImagePixelType_uint8);
AddRAMParameter();
}
void PrepareFilters(FloatVectorImageType * &imageToResample,
FloatVectorImageType * &imageToExtract,
FloatVectorImageType::RegionType imageToExtractRegion)
{
// Initialize resample filter
NNInterpolatorType::Pointer interpolator = NNInterpolatorType::New();
m_ResampleFilter = ResampleImageFilterType::New();
m_ResampleFilter->SetInput(imageToResample);
m_ResampleFilter->SetInterpolator(interpolator);
// Initialize roi extract filter
m_ExtractROIFilter = ExtractROIFilterType::New();
m_ExtractROIFilter->SetInput(imageToExtract);
m_ExtractROIFilter->SetExtractionRegion(imageToExtractRegion);
m_ExtractROIFilter->UpdateOutputInformation();
// Set the resample filter with extracted image origin, spacing, and size
m_ResampleFilter->SetOutputOrigin(m_ExtractROIFilter->GetOutput()->GetOrigin());
m_ResampleFilter->SetOutputSpacing(m_ExtractROIFilter->GetOutput()->GetSpacing());
m_ResampleFilter->SetOutputSize(
m_ExtractROIFilter->GetOutput()->GetLargestPossibleRegion().GetSize());
m_ResampleFilter->UpdateOutputInformation();
}
void DoExecute()
{
// Get input images pointers
FloatVectorImageType* t0 = GetParameterImage("inb");
FloatVectorImageType* t1 = GetParameterImage("ina");
// Use input image mask for image b (t0)
if (HasValue("inbmask"))
{
m_MaskImageFilter0 = MaskVectorImageFilterType::New();
m_MaskImageFilter0->SetMaskImage(GetParameterUInt8Image("inbmask"));
m_MaskImageFilter0->SetInput(GetParameterFloatVectorImage("inb"));
m_MaskImageFilter0->UpdateOutputInformation();
t0 = m_MaskImageFilter0->GetOutput();
}
// Use input image mask for image a (t1)
if (HasValue("inamask"))
{
m_MaskImageFilter1 = MaskVectorImageFilterType::New();
m_MaskImageFilter1->SetMaskImage(GetParameterUInt8Image("inamask"));
m_MaskImageFilter1->SetInput(GetParameterFloatVectorImage("ina"));
m_MaskImageFilter1->UpdateOutputInformation();
t1 = m_MaskImageFilter1->GetOutput();
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
}
// Compute rasters intersection region, check overlap
otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator;
comparator.SetImage1(t0);
comparator.SetImage2(t1);
if (!comparator.DoesOverlap())
{
otbAppLogFATAL("Input rasters do not overlap");
}
// Detect which input image (t0 or t1) have the smallest pixel
FloatVectorImageType::Pointer inputExtractedT0;
FloatVectorImageType::Pointer inputExtractedT1;
FloatVectorImageType::SpacingType spacingT0 = t0->GetSpacing();
FloatVectorImageType::SpacingType spacingT1 = t1->GetSpacing();
double pixelAreaTO = vnl_math_abs(spacingT0[0]*spacingT0[1]);
double pixelAreaT1 = vnl_math_abs(spacingT1[0]*spacingT1[1]);
if (pixelAreaTO > pixelAreaT1)
{
// Resample t0 over t1 and extract ROI (overlap) of t1
otbAppLogINFO("inb-->resampled");
otbAppLogINFO("ina-->extracted");
PrepareFilters(t0, t1, comparator.GetOverlapInImage2Indices());
inputExtractedT0 = m_ResampleFilter->GetOutput();
inputExtractedT1 = m_ExtractROIFilter->GetOutput();
}
else
{
// Resample t1 over t0 and extract ROI (overlap) of t0
otbAppLogINFO("inb-->extracted");
otbAppLogINFO("ina-->resampled");
PrepareFilters(t1, t0, comparator.GetOverlapInImage1Indices());
inputExtractedT0 = m_ExtractROIFilter->GetOutput();
inputExtractedT1 = m_ResampleFilter->GetOutput();
}
// Compute Delta NDVI
m_DeltaNDVIFilter = DeltaNDVIFilterType::New();
m_DeltaNDVIFilter->SetInput1(inputExtractedT0);
m_DeltaNDVIFilter->SetInput2(inputExtractedT1);
m_DeltaNDVIFilter->GetFunctor().SetNIRChannelT0(GetParameterInt("nirb"));
m_DeltaNDVIFilter->GetFunctor().SetRedChannelT0(GetParameterInt("redb"));
m_DeltaNDVIFilter->GetFunctor().SetNIRChannelT1(GetParameterInt("nira"));
m_DeltaNDVIFilter->GetFunctor().SetRedChannelT1(GetParameterInt("reda"));
// Stats filter
m_StatsFilter = StatsFilterType::New();
m_StatsFilter->SetIgnoreUserDefinedValue(true);
m_StatsFilter->SetUserIgnoredValue(m_DeltaNDVIFilter->GetFunctor().GetNoDataValue());
FloatImageType * deltaNDVIImage = m_DeltaNDVIFilter->GetOutput();
// Mask directory
if (HasValue("masksdir") || HasValue("mask"))
{
// Mask Delta NDVI image
m_MaskImageFilter = MaskImageFilterType::New();
m_MaskImageFilter->SetInput(m_DeltaNDVIFilter->GetOutput());
m_MaskImageFilter->SetOutsideValue(3.0);
deltaNDVIImage = m_MaskImageFilter->GetOutput();
// If there is a directory specified for masks, we instanciate the handler
if (HasValue("masksdir"))
{
otbAppLogINFO("Using specified directory for vegetation masks");
// Instanciate a mosaic from directory handler
m_MaskHandler = MaskHandlerType::New();
m_MaskHandler->SetDirectory(GetParameterAsString("masksdir"));
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
m_MaskHandler->SetOutputOrigin(m_ExtractROIFilter->GetOutput()->GetOrigin());
m_MaskHandler->SetOutputSpacing(m_ExtractROIFilter->GetOutput()->GetSpacing());
m_MaskHandler->SetOutputSize(m_ExtractROIFilter->GetOutput()->GetLargestPossibleRegion().GetSize());
m_MaskHandler->UpdateOutputInformation();
}
// If there is an input mask image, we need to resample it
if (HasValue("mask"))
{
otbAppLogINFO("Using additional input mask");
// Resample the mask
NNMaskImageInterpolatorType::Pointer maskInterpolator = NNMaskImageInterpolatorType::New();
m_MaskResampleFilter = ResampleMaskImageFilterType::New();
m_MaskResampleFilter->SetInput(GetParameterUInt8Image("mask"));
m_MaskResampleFilter->SetOutputOrigin(m_ExtractROIFilter->GetOutput()->GetOrigin());
m_MaskResampleFilter->SetOutputSpacing(m_ExtractROIFilter->GetOutput()->GetSpacing());
m_MaskResampleFilter->SetOutputSize(m_ExtractROIFilter->GetOutput()->GetLargestPossibleRegion().GetSize());
m_MaskResampleFilter->SetInterpolator(maskInterpolator);
m_MaskResampleFilter->UpdateOutputInformation();
}
// Wiring final mask to mask image filter
if (HasValue("masksdir") && HasValue("mask"))
{
// If both masks types are specified, we need to instantiate a AndImageFilter
m_AndFilter = AndFilterType::New();
m_AndFilter->SetInput1(m_MaskHandler->GetOutput());
m_AndFilter->SetInput2(m_MaskResampleFilter->GetOutput());
m_AndFilter->UpdateOutputInformation();
m_MaskImageFilter->SetMaskImage(m_AndFilter->GetOutput());
}
else if (HasValue("masksdir"))
{
m_MaskImageFilter->SetMaskImage(m_MaskHandler->GetOutput());
}
else // HasValue("mask")
{
m_MaskImageFilter->SetMaskImage(m_MaskResampleFilter->GetOutput());
}
}
// Wire stats filter
m_StatsFilter->SetInput(deltaNDVIImage);
// Compute stats
AddProcess(m_StatsFilter->GetStreamer(),"Computing dNDVI statistics");
m_StatsFilter->Update();
if (HasValue("outdndvi"))
{
// Connect delta NDVI to output
SetParameterOutputImage("outdndvi", deltaNDVIImage);
}
// Label image output
m_NDVILabelFilter = NDVILabelImageFilterType::New();
m_NDVILabelFilter->SetInput(deltaNDVIImage);
m_NDVILabelFilter->SetInputMeanObject(m_StatsFilter->GetMeanOutput());
m_NDVILabelFilter->SetInputSigmaObject(m_StatsFilter->GetSigmaOutput());
m_NDVILabelFilter->SetInputNoDataValue(m_DeltaNDVIFilter->GetFunctor().GetNoDataValue());
SetParameterOutputImage("outlabel", m_NDVILabelFilter->GetOutput());
}
ResampleImageFilterType::Pointer m_ResampleFilter;
ExtractROIFilterType::Pointer m_ExtractROIFilter;
ResampleMaskImageFilterType::Pointer m_MaskResampleFilter;
351352353354355356357358359360361362363364365
DeltaNDVIFilterType::Pointer m_DeltaNDVIFilter;
MaskImageFilterType::Pointer m_MaskImageFilter;
MaskVectorImageFilterType::Pointer m_MaskImageFilter0;
MaskVectorImageFilterType::Pointer m_MaskImageFilter1;
NDVILabelImageFilterType::Pointer m_NDVILabelFilter;
StatsFilterType::Pointer m_StatsFilter;
MaskHandlerType::Pointer m_MaskHandler;
AndFilterType::Pointer m_AndFilter;
};
}
}
OTB_APPLICATION_EXPORT( otb::Wrapper::ClearCutsDetection )