Skip to content
GitLab
Select Git revision
  • master default
  • geosud
  • no_toolsmodule
  • v3.0
  • V2
Find file BlameHistoryPermalink
otbClearCutsDetection.cxx 14.43 KiB
1 
/*=========================================================================
2
3 
  Copyright (c) IRSTEA. All rights reserved.
4
5
6 
     This software is distributed WITHOUT ANY WARRANTY; without even
7 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 
     PURPOSE.  See the above copyright notices for more information.
9
10 
=========================================================================*/
11 
#include "otbDeltaNDVIFunctor.h"
12 
#include "itkFixedArray.h"
13 
#include "itkObjectFactory.h"
14
15 
// Elevation handler
16 
#include "otbWrapperElevationParametersHandler.h"
17 
#include "otbWrapperApplicationFactory.h"
18
19 
// Application engine
20 
#include "otbStandardFilterWatcher.h"
21 
#include "itkFixedArray.h"
22
23 
// Filters
24 
#include "itkMaskImageFilter.h"
25 
#include "otbStreamingStatisticsImageFilter.h"
26 
#include "otbStreamingResampleImageFilter.h"
27 
#include "otbMultiChannelExtractROI.h"
28 
#include "otbDeltaNDVILabelerFilter.h"
29 
#include "itkAndImageFilter.h"
30
31 
// Helper
32 
#include "otbRegionComparator.h"
33
34 
// Mask handler
35 
#include "otbMosaicFromDirectoryHandler.h"
36
37 
namespace otb
38 
{
39
40 
namespace Wrapper
41 
{
42
43 
class ClearCutsDetection : public Application
44 
{
45 
public:
46 
  /** Standard class typedefs. */
47 
  typedef ClearCutsDetection            Self;
48 
  typedef Application                   Superclass;
49 
  typedef itk::SmartPointer<Self>       Pointer;
50 
  typedef itk::SmartPointer<const Self> ConstPointer;
51
52 
  /** Standard macro */
53 
  itkNewMacro(Self);
54 
  itkTypeMacro(ClearCutsDetection, Application);
55
56 
  /** Filters */
57 
  typedef otb::StreamingResampleImageFilter<FloatVectorImageType,
58 
      FloatVectorImageType> ResampleImageFilterType;
59 
  typedef otb::StreamingResampleImageFilter<UInt8ImageType,
60 
      UInt8ImageType> ResampleMaskImageFilterType;
61 
  typedef itk::NearestNeighborInterpolateImageFunction<
62 
      FloatVectorImageType> NNInterpolatorType;
63 
  typedef itk::NearestNeighborInterpolateImageFunction<
64 
      UInt8ImageType> NNMaskImageInterpolatorType;
65 
  typedef otb::Functor::DeltaNDVIFromChannels<FloatVectorImageType::PixelType,
66 
      FloatImageType::PixelType> DeltaNDVIFunctorType;
67 
  typedef itk::BinaryFunctorImageFilter<FloatVectorImageType,
68 
      FloatVectorImageType, FloatImageType, DeltaNDVIFunctorType> DeltaNDVIFilterType;
69 
  typedef itk::MaskImageFilter<FloatImageType, UInt8ImageType,
70 
      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 )