File: otbBandMathX.cxx

package info (click to toggle)
otb 5.8.0%2Bdfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 38,496 kB
  • ctags: 40,282
  • sloc: cpp: 306,573; ansic: 3,575; python: 450; sh: 214; perl: 74; java: 72; makefile: 70
file content (328 lines) | stat: -rw-r--r-- 14,424 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
/*=========================================================================

 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.

 =========================================================================*/
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "otbObjectList.h"

#include "otbBandMathXImageFilter.h"
#include "otbMultiChannelExtractROI.h"
#include "itksys/SystemTools.hxx"

namespace otb
{

namespace Wrapper
{

class BandMathX : public Application
{

public:
  /** Standard class typedefs. */
  typedef BandMathX                      Self;
  typedef Application                    Superclass;
  typedef itk::SmartPointer<Self>        Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

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

  itkTypeMacro(BandMathX, otb::Application);

  typedef otb::MultiChannelExtractROI
    <FloatVectorImageType::InternalPixelType,
     FloatVectorImageType::InternalPixelType> MultiChannelExtractorType;
  typedef otb::BandMathXImageFilter<FloatVectorImageType>   BandMathImageFilterType;


private:
  void DoInit() ITK_OVERRIDE
  {
    SetName("BandMathX");
    SetDescription("This application performs mathematical operations on multiband images.\n"
      "Mathematical formula interpretation is done via muParserX library : http://articles.beltoforion.de/article.php?a=muparserx");

    SetDocName("Band Math X");
    SetDocLongDescription("The goal of this documentation is to give the user some hints about the syntax used in this application.\n"
      "The syntax is mainly constrained by the muparserx library, which can be considered as the core of the application.\n"
      "\n\n"
      "- Fundamentals:\n\n"
      "The default prefix name for variables related to the ith input is 'im(i+1)'(note the indexing from 1 to N, for N inputs). \n"
      "The following list summaries the available variables for input #0 (and so on for every input): \n"
      "\n"
      "im1                                -->   a pixel from first input, made of n components (n bands)\n"
      "im1bj                             -->   jth component of a pixel from first input (first band is indexed by 1)\n"
      "im1bjNkxp                     -->   a neighbourhood ('N') of pixels of the jth component from first input, of size kxp\n"
      "im1PhyX and im1PhyY  -->   spacing of first input in X and Y directions (horizontal and vertical)\n"
      "im1bjMean im1bjMin im1bjMax im1bjSum im1bjVar  -->   mean, min, max, sum, variance of jth band from first input (global statistics)\n"
      "\nMoreover, we also have the following generic variables:\n"
      "idxX and idxY        -->   indices of the current pixel\n"
      "\n"
      "Always keep in mind that this application only addresses mathematically well-defined formulas.\n"
      "For instance, it is not possible to add vectors of different dimensions (this implies the addition of a row vector with a column vector),\n"
      "or add a scalar to a vector or a matrix, or divide two vectors, and so on...\n"
      "Thus, it is important to remember that a pixel of n components is always represented as a row vector.\n"
      "\n"
      "Example :\n\n"
      "                   im1 + im2       (1)\n"
      "\nrepresents the addition of pixels from first and second inputs. This expression is consistent only if\n"
      "both inputs have the same number of bands.\n"
      "Note that it is also possible to use the following expressions to obtain the same result:\n"
      "\n"
      "                   im1b1 + im2b1 \n"
      "                   im1b2 + im2b2       (2)\n"
      "                   ..."
      "\n\nNevertheless, the first expression is by far much pleaseant. We call this new functionality the 'batch mode'\n"
      "(performing the same operation in a band-to-band fashion).\n"

      "\n\n"
      "- Operations involving neighborhoods of pixels:\n\n"
      "Another new fonctionnality is the possibility to perform operations that involve neighborhoods of pixels.\n"
      "Variable related to such neighborhoods are always defined following the pattern imIbJNKxP, where: \n"
      "- I is an number identifying the image input (remember, input #0 = im1, and so on)\n"
      "- J is an number identifying the band (remember, first band is indexed by 1)\n"
      "- KxP are two numbers that represent the size of the neighborhood (first one is related to the horizontal direction)\n"
      "All neighborhood are centered, thus K and P must be odd numbers.\n"
      "Many operators come with this new functionality: dotpr, mean var median min max...\n"
      "For instance, if im1 represents the pixel of 3 bands image:\n\n"
      "               im1 - mean(im1b1N5x5,im1b2N5x5,im1b3N5x5)       (3)\n"
      "\ncould represent a high pass filter (Note that by implying three neighborhoods, the operator mean returns a row vector of three components.\n"
      "It is a typical behaviour for many operators of this application).\n"

      "\n\n"
      "- Operators:\n\n"
      "In addition to the previous operators, other operators are available:\n"
      "- existing operators/functions from muParserX, that were not originally defined for vectors and\n"
      "matrices (for instance cos, sin, ...). These new operators/ functions keep the original names to which we added the prefix 'v' for vector (vcos, vsin, ...).\n"
      "- mult, div and pow operators, that perform element-wise multiplication, division or exponentiation of vector/matrices (for instance im1 div im2)\n"
      "- mlt, dv and pw operators, that perform multiplication, division or exponentiation of vector/matrices by a scalar (for instance im1 dv 2.0)\n"
      "- bands, which is a very useful operator. It allows selecting specific bands from an image, and/or to rearrange them in a new vector;\n"
      "for instance bands(im1,{1,2,1,1}) produces a vector of 4 components made of band 1, band 2, band 1 and band 1 values from the first input.\n"
      "Note that curly brackets must be used in order to select the desired band indices.\n"
      "... and so on.\n"

      "\n\n"
      "- Application itself:\n\n"
      "The application takes the following parameters :\n"
      "- Setting the list of inputs can be done with the 'il' parameter.\n"
      "- Setting expressions can be done with the 'exp' parameter (see also limitations section below).\n"
      "- Setting constants can be done with the 'incontext' parameter. User must provide a txt file with a specific syntax: #type name value\n"
      "An example of such a file is given below:\n\n"
      "#F expo 1.1\n"
      "#M kernel1 { 0.1 , 0.2 , 0.3; 0.4 , 0.5 , 0.6; 0.7 , 0.8 , 0.9; 1 , 1.1 , 1.2; 1.3 , 1.4 , 1.5 }\n"
      "\nAs we can see,  #I/#F allows the definition of an integer/float constant, whereas #M allows the definition of a vector/matrix.\n"
      "In the latter case, elements of a row must be separated by commas, and rows must be separated by semicolons.\n"
      "It is also possible to define expressions within the same txt file, with the pattern #E expr. For instance (two expressions; see also limitations section below):\n\n"
      "#E $dotpr(kernel1,im1b1N3x5); im2b1^expo$\n"
      "\n- The 'outcontext' parameter allows saving user's constants and expressions (context).\n"
      "- Setting the output image can be done with the 'out' parameter (multi-outputs is not implemented yet).\n"
      "\n\n"
      "Finally, we strongly recommend that the reader takes a look at the cookbook, where additional information can be found (http://www.orfeo-toolbox.org/packages/OTBCookBook.pdf).\n"

);

    SetDocLimitations("The application is currently unable to produce one output image per expression, contrary to otbBandMathXImageFilter.\n"
                      "Separating expressions by semi-colons (; ) will concatenate their results into a unique multiband output image. ");
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");
    AddDocTag("Util");

    AddParameter(ParameterType_InputImageList,  "il",   "Input image list");
    SetParameterDescription("il", "Image list to perform computation on.");

    AddParameter(ParameterType_OutputImage, "out", "Output Image");
    SetParameterDescription("out","Output image.");

    AddRAMParameter();

    AddParameter(ParameterType_String, "exp", "Expressions");
    SetParameterDescription("exp",
                            "Mathematical expression to apply.");

    AddParameter(ParameterType_InputFilename, "incontext", "Import context");
    SetParameterDescription("incontext",
                            "A txt file containing user's constants and expressions.");
    MandatoryOff("incontext");

    AddParameter(ParameterType_OutputFilename, "outcontext", "Export context");
    SetParameterDescription("outcontext",
                            "A txt file where to save user's constants and expressions.");
    MandatoryOff("outcontext");

    // Doc example parameter settings
    SetDocExampleParameterValue("il", "verySmallFSATSW_r.tif verySmallFSATSW_nir.tif verySmallFSATSW.tif");
    SetDocExampleParameterValue("out", "apTvUtBandMathOutput.tif");
    SetDocExampleParameterValue("exp", "\"cos(im1b1)+im2b1*im3b1-im3b2+ndvi(im3b3, im3b4)\"");
  }

  void DoUpdateParameters() ITK_OVERRIDE
  {
    // check if input context should be used
    bool useContext = this->ContextCheck();
    // Check if the expression is correctly set
    if (HasValue("il") && HasValue("exp"))
      {
      this->LiveCheck(useContext);
      }
  }

  bool ContextCheck(void)
    {
    bool useContext = false;
    if (IsParameterEnabled("incontext") && HasValue("incontext"))
      {
      std::string contextPath = GetParameterString("incontext");
      // check that file exists
      if (itksys::SystemTools::FileExists(contextPath.c_str(),true))
        {
        BandMathImageFilterType::Pointer dummyFilter =
          BandMathImageFilterType::New();
        dummyFilter->SetManyExpressions(false);
        try
          {
          dummyFilter->ImportContext(contextPath);
          useContext = true;
          }
        catch(itk::ExceptionObject& err)
          {
          //trick to prevent unreferenced local variable warning on MSVC
          (void)err;
          // silent catch
          useContext = false;
          }
        if (useContext)
          {
          // only set the first expression, 'ManyExpression' is disabled.
          this->SetParameterString("exp",dummyFilter->GetExpression(0));
          }
        }
      }
    return useContext;
    }

  void LiveCheck(bool useContext=false)
    {
    BandMathImageFilterType::Pointer dummyFilter =
      BandMathImageFilterType::New();
    dummyFilter->SetManyExpressions(false);

    std::vector<MultiChannelExtractorType::Pointer> extractors;
    FloatVectorImageListType::Pointer inList = GetParameterImageList("il");
    for (unsigned int i = 0; i < inList->Size(); i++)
      {
      FloatVectorImageType::Pointer inImg = inList->GetNthElement(i);
      FloatVectorImageType::RegionType  largestRegion = inImg->GetLargestPossibleRegion();
      unsigned int nbChannels = inImg->GetNumberOfComponentsPerPixel();

      MultiChannelExtractorType::Pointer extract = MultiChannelExtractorType::New();
      extractors.push_back(extract);
      extract->SetInput(inImg);
      extract->SetStartX(largestRegion.GetIndex(0));
      extract->SetStartY(largestRegion.GetIndex(1));
      // Set extract size to 1 in case of global stats computation
      extract->SetSizeX(1);
      extract->SetSizeY(1);
      for (unsigned int j=0 ; j<nbChannels ; ++j)
        {
        extract->SetChannel(j+1);
        }
      dummyFilter->SetNthInput(i,extract->GetOutput());
      }
    if (useContext)
      {
      dummyFilter->ImportContext(GetParameterString("incontext"));
      }
    else
      {
      dummyFilter->SetExpression(GetParameterString("exp"));
      }
    try
      {
      dummyFilter->UpdateOutputInformation();
      SetParameterDescription("exp", "Valid expression");
      }
    catch(itk::ExceptionObject& err)
      {
      // Change the parameter description to be able to have the
      // parser errors in the tooltip
      SetParameterDescription("exp", err.GetDescription());
      }
    }

  void DoExecute() ITK_OVERRIDE
  {
    // Get the input image list
    FloatVectorImageListType::Pointer inList = GetParameterImageList("il");

    // checking the input images list validity
    const unsigned int nbImages = inList->Size();

    if (nbImages == 0)
    {
     itkExceptionMacro("No input Image set...; please set at least one input image");
    }

    if ( (!IsParameterEnabled("exp")) && (!IsParameterEnabled("incontext")) )
    {
     itkExceptionMacro("No expression set...; please set and enable at least one one expression");
    }

    m_Filter               = BandMathImageFilterType::New();
    m_Filter->SetManyExpressions(false);

    for (unsigned int i = 0; i < nbImages; i++)
      {
        FloatVectorImageType::Pointer currentImage = inList->GetNthElement(i);
        currentImage->UpdateOutputInformation();

        otbAppLogINFO( << "Image #" << i + 1 << " has "
                       << currentImage->GetNumberOfComponentsPerPixel()
                       << " components");

        m_Filter->SetNthInput(i,currentImage);

      }

    bool useContext = this->ContextCheck();
    std::string expStr = GetParameterString("exp");
    if (useContext)
      {
      otbAppLogINFO("Using input context : " << expStr );
      m_Filter->ImportContext(GetParameterString("incontext"));
      }
    else
      {
      otbAppLogINFO("Using expression : " << expStr );
      m_Filter->SetExpression(expStr);
      }

    if ( IsParameterEnabled("outcontext") && HasValue("outcontext") )
      m_Filter->ExportContext(GetParameterString("outcontext"));

    // Set the output image
    SetParameterOutputImage("out", m_Filter->GetOutput());
  }

  BandMathImageFilterType::Pointer  m_Filter;
};

} // namespace Wrapper
} // namespace otb

OTB_APPLICATION_EXPORT(otb::Wrapper::BandMathX)