otbMorphologicalProfilesAnalysis.cxx 15.6 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
43



#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 "itkTimeProbe.h"
#include "otbConvexOrConcaveClassificationFilter.h"
#include "otbMorphologicalProfilesSegmentationFilter.h"
#include "otbGeodesicMorphologyIterativeDecompositionImageFilter.h"

44
45
46
47
namespace otb
{
namespace Wrapper
{
48

49
50
class MorphologicalProfilesAnalysis : public Application
{
51
52
53
54
55
56
57
58
59
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;

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

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

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

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

  itkTypeMacro( MorphologicalProfilesAnalysis, otb::Application );

private:

  void DoInit() ITK_OVERRIDE
  {
    SetName( "MorphologicalProfilesAnalysis" );
78
    SetDescription( "Performs morphological profiles analysis on an input image channel." );
79
80
81

    // Documentation
    SetDocName( "Morphological Profiles Analysis" );
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    SetDocLongDescription( "This algorithm is derived from the following publication:\n"
                                   "\n"
                                   "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"
                                   "February 2001, p. 309-320.\n"
                                   "\n"
                                   "Depending of the profile selection, the application provides::\n\n"
                                   "- 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"
                                   "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"
Guillaume Pasero's avatar
Guillaume Pasero committed
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
103
104
105
106
                                   "\n"
                                   "The output image can be :"
                                   "- 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"
                                   "- A labeled image for the classification\n" );
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
125
126

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

    AddRAMParameter();

127
    // Structuring Element (Ball | Cross)
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
    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 );

    AddParameter( ParameterType_Int, "step", "Radius step." );
    SetParameterDescription( "step", "Radius step along the profile (in pixels)" );
    SetDefaultParameterInt( "step", 1 );
    SetMinimumParameterIntValue( "step", 1 );


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

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

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
175
176
177
178
179
180
181
182
183
  }

  void DoUpdateParameters() ITK_OVERRIDE
  {
    // Nothing to do here : all parameters are independent
  }

  void DoExecute() ITK_OVERRIDE
  {

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

185
186
187
    int nBComp = inImage->GetNumberOfComponentsPerPixel();
    int selectedChannel = GetParameterInt( "channel" );

188
    if( selectedChannel > nBComp )
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
      {
      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 );
211
212
      }
    else // Cross
213
214
215
216
217
218
219
220
221
222
      {
      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) {

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

227
228
    typedef otb::MultiScaleConvexOrConcaveClassificationFilter<FloatImageType, LabeledImageType> MultiScaleClassificationFilterType;
    typedef otb::ProfileDerivativeToMultiScaleCharacteristicsFilter<FloatImageType, FloatImageType, LabeledImageType> MultiScaleCharacteristicsFilterType;
229
230

    // Instantiation
231
232
233
234
235
236
237
    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();
238
239
240
241
242
243

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

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

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

    if ( doClosing )
      {
      performOperations<ClosingProfileFilterType, DerivativeFilterType, MultiScaleCharacteristicsFilterType>(
              cprofileFilter, cderivativeFilter, cmsCharFilter,
264
              closing, derivativeClosing, characClosing,
265
266
267
268
269
270
271
272
273
274
275
276
              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) );
277
    AddProcess(classificationFilter, "Classification");
278
279
280
281
282
283
284
285
286
287
    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,
288
                    bool profile, bool derivative, bool characteristics,
289
290
                    unsigned int profileSize, unsigned short initValue, unsigned short step) {

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

    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() );
303
      AddProcess(listToVectorImageFilter, "Profile");
304
305
306
307
308
309
310
      listToVectorImageFilter->Update();
      SetParameterOutputImage( "out", listToVectorImageFilter->GetOutput() );
      return;
      }

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

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

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

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

  ExtractorFilterType::Pointer m_ExtractorFilter;

};
}
}

OTB_APPLICATION_EXPORT( otb::Wrapper::MorphologicalProfilesAnalysis )