otbMorphologicalProfilesAnalysis.cxx 15.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42



#include <itkComposeImageFilter.h>
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "itkBinaryBallStructuringElement.h"
#include "itkBinaryCrossStructuringElement.h"

#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryMorphologicalOpeningImageFilter.h"
#include "itkBinaryMorphologicalClosingImageFilter.h"

#include "otbMultiToMonoChannelExtractROI.h"
#include "otbImageList.h"
#include "otbImageListToVectorImageFilter.h"

#include "otbConvexOrConcaveClassificationFilter.h"
#include "otbMorphologicalProfilesSegmentationFilter.h"
#include "otbGeodesicMorphologyIterativeDecompositionImageFilter.h"

43 44 45 46
namespace otb
{
namespace Wrapper
{
47

48 49
class MorphologicalProfilesAnalysis : public Application
{
50 51 52 53 54 55 56 57 58
public:
/** Standard class typedefs. */
  typedef MorphologicalProfilesAnalysis Self;
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  typedef FloatVectorImageType::InternalPixelType InputPixelType;

59
  typedef unsigned short LabeledPixelType;
60 61
  typedef otb::Image<LabeledPixelType, 2> LabeledImageType;

62
  typedef otb::MultiToMonoChannelExtractROI<InputPixelType, InputPixelType> ExtractorFilterType;
63 64 65 66 67 68 69 70 71 72 73

  typedef itk::BinaryBallStructuringElement<InputPixelType, 2> BallStructuringElementType;
  typedef itk::BinaryCrossStructuringElement<InputPixelType, 2> CrossStructuringElementType;

/** Standard macro */
  itkNewMacro( Self );

  itkTypeMacro( MorphologicalProfilesAnalysis, otb::Application );

private:

74
  void DoInit() override
75 76
  {
    SetName( "MorphologicalProfilesAnalysis" );
77
    SetDescription( "Performs morphological profiles analysis on an input image channel." );
78 79 80

    // Documentation
    SetDocName( "Morphological Profiles Analysis" );
Victor Poughon's avatar
Victor Poughon committed
81 82
    SetDocLongDescription( "This algorithm is derived from the following publication:\n\n"

83 84 85
                                   "Martino Pesaresi and Jon Alti Benediktsson, Member, IEEE: A new approach\n"
                                   "for the morphological segmentation of high resolution satellite imagery.\n"
                                   "IEEE Transactions on geoscience and remote sensing, vol. 39, NO. 2,\n"
Victor Poughon's avatar
Victor Poughon committed
86 87
                                   "February 2001, p. 309-320.\n\n"

88
                                   "Depending of the profile selection, the application provides::\n\n"
Victor Poughon's avatar
Victor Poughon committed
89 90 91 92 93 94 95 96

                                   "* The multi scale geodesic morphological opening or closing profile of the input image.\n"
                                   "* The multi scale derivative of the opening or closing profile.\n"
                                   "* The parameter (called characteristic) of the maximum derivative value of the multi scale closing or opening profile for which this maxima occurs.\n"
                                   "* The labeled classification of the input image.\n\n"

                                   "The behavior of the classification is:\n\n"

97 98 99 100
                                   "Given :math:`x_1` and :math:`x_2` two membership values,\n"
                                   ":math:`L_1, L_2` two labels associated, and :math:`\\sigma` a tolerance\n"
                                   "value, the following decision rule is applied:\n"
                                   "\n"
101
                                   ":math:`L = \\begin{cases}  L_{1} & : x_{1}>x_{2} \\quad and \\quad x_{1}>\\sigma \\\\ L_{2} & : x_{2}>x_{1} \\quad and \\quad x_{2}>\\sigma \\\\ 0   & : otherwise. \\end{cases}` \n"
102
                                   "\n"
Victor Poughon's avatar
Victor Poughon committed
103
                                   "The output image can be:"
104 105
                                   "- A :math:`N` multi band image for the opening/closing normal or derivative profiles.\n"
                                   "- A mono band image for the opening/closing characteristics.\n"
106
                                   "- A labeled image for the classification." );
107
    SetDocLimitations( "Generation of the morphological profile is not streamable, pay attention to this fact when setting the radius initial size and step of the structuring element." );
108
    SetDocAuthors( "OTB-Team" );
109
    SetDocSeeAlso( "otbMorphologicalOpeningProfileFilter, otbMorphologicalClosingProfileFilter, otbProfileToProfileDerivativeFilter, otbProfileDerivativeToMultiScaleCharacteristicsFilter, otbMultiScaleConvexOrConcaveClassificationFilter, classes" );
110

111 112
    AddDocTag(Tags::FeatureExtraction);
    AddDocTag("Morphology");
113 114

    AddParameter( ParameterType_InputImage, "in", "Input Image" );
115
    SetParameterDescription( "in", "The input image." );
116 117

    AddParameter( ParameterType_OutputImage, "out", "Output Image" );
118
    SetParameterDescription( "out", "The output image." );
119 120 121 122 123 124

    AddParameter( ParameterType_Int, "channel", "Selected Channel" );
    SetParameterDescription( "channel", "The selected channel index for input image" );
    SetDefaultParameterInt( "channel", 1 );
    SetMinimumParameterIntValue( "channel", 1 );

125
    // Structuring Element (Ball | Cross)
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
    AddParameter( ParameterType_Choice, "structype", "Structuring Element Type" );
    SetParameterDescription( "structype", "Choice of the structuring element type" );
    AddChoice( "structype.ball", "Ball" );
    AddChoice( "structype.cross", "Cross" );

    AddParameter( ParameterType_Int, "size", "Profile Size" );
    SetParameterDescription( "size", "Size of the profiles" );
    SetDefaultParameterInt( "size", 5 );
    SetMinimumParameterIntValue( "size", 2 );

    AddParameter( ParameterType_Int, "radius", "Initial radius" );
    SetParameterDescription( "radius", "Initial radius of the structuring element (in pixels)" );
    SetDefaultParameterInt( "radius", 5 );
    SetMinimumParameterIntValue( "radius", 1 );

Victor Poughon's avatar
Victor Poughon committed
141
    AddParameter( ParameterType_Int, "step", "Radius step" );
142 143 144 145 146
    SetParameterDescription( "step", "Radius step along the profile (in pixels)" );
    SetDefaultParameterInt( "step", 1 );
    SetMinimumParameterIntValue( "step", 1 );


147
    AddParameter( ParameterType_Choice, "profile", "Profile" );
148 149 150
    SetParameterDescription( "profile", "" );
    AddChoice( "profile.opening", "opening" );
    AddChoice( "profile.closing", "closing" );
151 152
    AddChoice( "profile.derivativeopening", "derivativeopening" );
    AddChoice( "profile.derivativeclosing", "derivativeclosing" );
153 154 155 156
    AddChoice( "profile.openingcharacteristics", "openingcharacteristics" );
    AddChoice( "profile.closingcharacteristics", "closingcharacteristics" );
    AddChoice( "profile.classification", "classification" );

157 158
    AddParameter( ParameterType_Float, "profile.classification.sigma", "Sigma value for leveling tolerance" );
    SetParameterDescription( "profile.classification.sigma", "Sigma value for leveling tolerance" );
159 160 161
    SetDefaultParameterFloat( "profile.classification.sigma", 1 );
    SetMinimumParameterFloatValue( "profile.classification.sigma", 0 );

Victor Poughon's avatar
Victor Poughon committed
162 163
    AddRAMParameter();

164 165 166 167 168 169 170 171 172
    SetDocExampleParameterValue( "in", "ROI_IKO_PAN_LesHalles.tif" );
    SetDocExampleParameterValue( "channel", "1" );
    SetDocExampleParameterValue( "structype", "ball" );
    SetDocExampleParameterValue( "profile", "classification" );
    SetDocExampleParameterValue( "size", "5" );
    SetDocExampleParameterValue( "radius", "1" );
    SetDocExampleParameterValue( "step", "1" );
    SetDocExampleParameterValue( "profile.classification.sigma", "1" );
    SetDocExampleParameterValue( "out", "output.tif" );
173

174
    SetOfficialDocLink();
175 176
  }

177
  void DoUpdateParameters() override
178
  {
Victor Poughon's avatar
Victor Poughon committed
179
    // Nothing to do here: all parameters are independent
180 181
  }

182
  void DoExecute() override
183 184 185
  {

    FloatVectorImageType::Pointer inImage = GetParameterImage( "in" );
186

187 188 189
    int nBComp = inImage->GetNumberOfComponentsPerPixel();
    int selectedChannel = GetParameterInt( "channel" );

190
    if( selectedChannel > nBComp )
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
      {
      itkExceptionMacro( << "The specified channel index for input image is invalid." );
      }

    m_ExtractorFilter = ExtractorFilterType::New();
    m_ExtractorFilter->SetInput( inImage );
    m_ExtractorFilter->SetStartX( static_cast<unsigned int>(inImage->GetLargestPossibleRegion().GetIndex( 0 )) );
    m_ExtractorFilter->SetStartY( static_cast<unsigned int>(inImage->GetLargestPossibleRegion().GetIndex( 1 )) );
    m_ExtractorFilter->SetSizeX( inImage->GetLargestPossibleRegion().GetSize( 0 ) );
    m_ExtractorFilter->SetSizeY( inImage->GetLargestPossibleRegion().GetSize( 1 ) );
    m_ExtractorFilter->SetChannel( static_cast<unsigned int>(GetParameterInt( "channel" )) );

    unsigned int profileSize = static_cast<unsigned int>(GetParameterInt( "size" ));
    unsigned short initValue = static_cast<unsigned short>(GetParameterInt( "radius" ));
    unsigned short step = static_cast<unsigned short>(GetParameterInt( "step" ));
    float sigma = GetParameterFloat( "profile.classification.sigma" );
    std::string profile = GetParameterString( "profile" );


    if ( GetParameterString( "structype" ) == "ball" )
      {
      performProfileAnalysis<BallStructuringElementType>( profile, profileSize, initValue, step, sigma );
213 214
      }
    else // Cross
215 216 217 218 219 220 221 222 223 224
      {
      performProfileAnalysis<CrossStructuringElementType>( profile, profileSize, initValue, step, sigma );
      }
  }

  template<typename StructuringElementType>
  void
  performProfileAnalysis(std::string profile, unsigned int profileSize, unsigned short initValue,
                         unsigned short step, float sigma) {

225 226 227
    typedef otb::MorphologicalOpeningProfileFilter<FloatImageType, FloatImageType, StructuringElementType> OpeningProfileFilterType;
    typedef otb::MorphologicalClosingProfileFilter<FloatImageType, FloatImageType, StructuringElementType> ClosingProfileFilterType;
    typedef otb::ProfileToProfileDerivativeFilter<FloatImageType, FloatImageType> DerivativeFilterType;
228

229 230
    typedef otb::MultiScaleConvexOrConcaveClassificationFilter<FloatImageType, LabeledImageType> MultiScaleClassificationFilterType;
    typedef otb::ProfileDerivativeToMultiScaleCharacteristicsFilter<FloatImageType, FloatImageType, LabeledImageType> MultiScaleCharacteristicsFilterType;
231 232

    // Instantiation
233 234 235 236 237 238 239
    typename OpeningProfileFilterType::Pointer oprofileFilter = OpeningProfileFilterType::New();
    typename ClosingProfileFilterType::Pointer cprofileFilter = ClosingProfileFilterType::New();
    typename DerivativeFilterType::Pointer oderivativeFilter = DerivativeFilterType::New();
    typename DerivativeFilterType::Pointer cderivativeFilter = DerivativeFilterType::New();
    typename MultiScaleCharacteristicsFilterType::Pointer omsCharFilter = MultiScaleCharacteristicsFilterType::New();
    typename MultiScaleCharacteristicsFilterType::Pointer cmsCharFilter = MultiScaleCharacteristicsFilterType::New();
    typename MultiScaleClassificationFilterType::Pointer classificationFilter = MultiScaleClassificationFilterType::New();
240 241 242 243 244 245

    bool classify = profile == "classification";
    bool opening = profile == "opening";
    bool closing = profile == "closing";
    bool characOpening = profile == "openingcharacteristics";
    bool characClosing = profile == "closingcharacteristics";
246 247
    bool derivativeOpening = profile == "derivativeopening";
    bool derivativeClosing = profile == "derivativeclosing";
248

249 250
    bool doOpening = classify || opening || derivativeOpening || characOpening;
    bool doClosing = classify || closing || derivativeClosing || characClosing;
251 252 253 254 255

    if ( doOpening )
      {
      performOperations<OpeningProfileFilterType, DerivativeFilterType, MultiScaleCharacteristicsFilterType>(
              oprofileFilter, oderivativeFilter, omsCharFilter,
256
              opening, derivativeOpening, characOpening,
257 258 259 260 261 262 263 264 265
              profileSize, step, initValue );
      if ( !classify )
        return;
      }

    if ( doClosing )
      {
      performOperations<ClosingProfileFilterType, DerivativeFilterType, MultiScaleCharacteristicsFilterType>(
              cprofileFilter, cderivativeFilter, cmsCharFilter,
266
              closing, derivativeClosing, characClosing,
267 268 269 270 271 272 273 274 275 276 277 278
              profileSize, step, initValue );
      if ( !classify )
        return;
      }

    classificationFilter = MultiScaleClassificationFilterType::New();
    classificationFilter->SetOpeningProfileDerivativeMaxima( omsCharFilter->GetOutput() );
    classificationFilter->SetOpeningProfileCharacteristics( omsCharFilter->GetOutputCharacteristics() );
    classificationFilter->SetClosingProfileDerivativeMaxima( cmsCharFilter->GetOutput() );
    classificationFilter->SetClosingProfileCharacteristics( cmsCharFilter->GetOutputCharacteristics() );
    classificationFilter->SetSigma( sigma );
    classificationFilter->SetLabelSeparator( static_cast<unsigned short>(initValue + profileSize * step) );
279
    AddProcess(classificationFilter, "Classification");
280 281 282 283 284 285 286 287 288 289
    classificationFilter->Update();
    SetParameterOutputImage( "out", classificationFilter->GetOutput() );
  }


  template<typename TProfileFilter, typename TDerivativeFilter, typename TCharacteristicsFilter>
  void
  performOperations(typename TProfileFilter::Pointer &profileFilter,
                    typename TDerivativeFilter::Pointer &derivativeFilter,
                    typename TCharacteristicsFilter::Pointer &msCharFilter,
290
                    bool profile, bool derivative, bool characteristics,
291 292
                    unsigned int profileSize, unsigned short initValue, unsigned short step) {

293 294
    typedef ImageList<FloatImageType> TImageList;
    typedef otb::ImageListToVectorImageFilter<TImageList, FloatVectorImageType> TListToVectorImageFilter;
295 296 297 298 299 300 301 302 303 304

    profileFilter->SetInput( m_ExtractorFilter->GetOutput() );
    profileFilter->SetProfileSize( profileSize );
    profileFilter->SetInitialValue( initValue );
    profileFilter->SetStep( step );

    if ( profile )
      {
      TListToVectorImageFilter::Pointer listToVectorImageFilter = TListToVectorImageFilter::New();
      listToVectorImageFilter->SetInput( profileFilter->GetOutput() );
305
      AddProcess(listToVectorImageFilter, "Profile");
306 307 308 309 310 311 312
      listToVectorImageFilter->Update();
      SetParameterOutputImage( "out", listToVectorImageFilter->GetOutput() );
      return;
      }

    derivativeFilter->SetInput( profileFilter->GetOutput() );

313
    if ( derivative )
314 315 316
      {
      TListToVectorImageFilter::Pointer listToVectorImageFilter = TListToVectorImageFilter::New();
      listToVectorImageFilter->SetInput( derivativeFilter->GetOutput() );
317
      AddProcess(listToVectorImageFilter, "Derivative");
318 319 320 321 322 323 324 325 326 327 328
      listToVectorImageFilter->Update();
      SetParameterOutputImage( "out", listToVectorImageFilter->GetOutput() );
      return;
      }

    msCharFilter->SetInput( derivativeFilter->GetOutput() );
    msCharFilter->SetInitialValue( initValue );
    msCharFilter->SetStep( step );

    if ( characteristics )
      {
329
      AddProcess(msCharFilter, "Characteristics");
330
      msCharFilter->Update();
331
      SetParameterOutputImage( "out", msCharFilter->GetOutputCharacteristics() );
332 333 334 335 336 337 338 339 340 341 342
      }
  }

  ExtractorFilterType::Pointer m_ExtractorFilter;

};
}
}

OTB_APPLICATION_EXPORT( otb::Wrapper::MorphologicalProfilesAnalysis )