otbDynamicConvert.cxx 19.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
/*
 * 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.
 */

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "otbVectorRescaleIntensityImageFilter.h"
25
#include "otbFunctorImageFilter.h"
26 27 28 29 30 31 32 33 34
#include "otbStreamingShrinkImageFilter.h"
#include "itkListSample.h"
#include "otbListSampleToHistogramListGenerator.h"
#include "itkImageRegionConstIterator.h"

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

35 36
#include <numeric>

37 38 39 40 41 42
namespace otb
{
namespace Wrapper
{


43 44 45
class DynamicConvert : public Application
{
public:
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  /** Standard class typedefs. */
  typedef DynamicConvert                Self;
  typedef Application                   Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

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

  itkTypeMacro(DynamicConvert, otb::Application);

  /** Filters typedef */
  typedef itk::Statistics::ListSample<FloatVectorImageType::PixelType> ListSampleType;
  typedef itk::Statistics::DenseFrequencyContainer2 DFContainerType;
  typedef ListSampleToHistogramListGenerator<ListSampleType,
61 62 63
    FloatVectorImageType::InternalPixelType,
    DFContainerType> HistogramsGeneratorType;
    
64
  typedef StreamingShrinkImageFilter<FloatVectorImageType,
65 66
    FloatVectorImageType> ShrinkFilterType;

67 68
  typedef StreamingShrinkImageFilter<UInt8ImageType, UInt8ImageType> UInt8ShrinkFilterType;

69 70
private:

71
  void DoInit() override
72 73
  {
    SetName("DynamicConvert");
74
    SetDescription("Change the pixel type and rescale the image's dynamic");
75 76

    SetDocName("Dynamic Conversion");
77 78
    SetDocLongDescription(
      "This application performs an image pixel type "
79 80
      "conversion (short, ushort, uchar, int, uint, float and double types are "
      "handled). The output image is written in the specified format (ie. "
81 82
      "that corresponds to the given extension).\n"
      "The conversion can include a rescale of the data range, by default it's set between the 2nd to "
83 84
      "the 98th percentile. The rescale can be linear or log2. \n"
      "The choice of the output channels can be done with the extended filename, but "
85 86
      "less easy to handle. To do this, a 'channels' parameter allows you to "
      "select the desired bands at the output. There are 3 modes, the "
87 88
      "available choices are: \n\n"

89 90
      "* **All**: keep all bands.\n"
      "* **Grayscale**: to display mono image as standard color image \n"
91
      "* **RGB**: select 3 bands in the input image (multi-bands)\n"
92
    );
93
    SetDocLimitations("The application does not support complex pixel types as output.");
94
    SetDocAuthors("OTB-Team");
95
    SetDocSeeAlso("Rescale");
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111

    AddDocTag(Tags::Manip);
    AddDocTag("Conversion");
    AddDocTag("Image Dynamic");

    AddParameter(ParameterType_InputImage, "in", "Input image");
    SetParameterDescription("in", "Input image");

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

    AddParameter(ParameterType_Choice, "type", "Rescale type");
    SetParameterDescription("type", "Transfer function for the rescaling");
    AddChoice("type.linear", "Linear");
    AddChoice("type.log2", "Log2");
112
    SetParameterString("type", "linear");
113

114 115 116 117
    AddParameter(ParameterType_Float,"type.linear.gamma",
      "Gamma correction factor");
    SetParameterDescription("type.linear.gamma",
      "Gamma correction factor");
118 119 120 121
    SetDefaultParameterFloat("type.linear.gamma",1.0);
    MandatoryOff("type.linear.gamma");

    AddParameter(ParameterType_InputImage,  "mask",   "Input mask");
122
    SetParameterDescription("mask",
123 124
      "Optional mask to indicate which pixels are valid for computing the histogram quantiles. "
      "Pixels where the mask is zero will not contribute to the histogram. "
125
      "The mask must have the same dimensions as the input image.");
126 127 128 129
    MandatoryOff("mask");
    DisableParameter("mask");

    AddParameter(ParameterType_Group,"quantile","Histogram quantile cutting");
130 131
    SetParameterDescription("quantile",
      "Cut the histogram edges before rescaling");
132 133

    AddParameter(ParameterType_Float, "quantile.high", "High cut quantile");
134 135
    SetParameterDescription("quantile.high", 
      "Quantiles to cut from histogram high values "
136 137 138 139 140 141
      "before computing min/max rescaling (in percent, 2 by default)");
    MandatoryOff("quantile.high");
    SetDefaultParameterFloat("quantile.high", 2.0);
    DisableParameter("quantile.high");

    AddParameter(ParameterType_Float, "quantile.low", "Low cut quantile");
142 143
    SetParameterDescription("quantile.low", 
      "Quantiles to cut from histogram low values "
144 145 146 147 148 149 150 151 152 153
      "before computing min/max rescaling (in percent, 2 by default)");
    MandatoryOff("quantile.low");
    SetDefaultParameterFloat("quantile.low", 2.0);
    DisableParameter("quantile.low");

    AddParameter(ParameterType_Choice, "channels", "Channels selection");
    SetParameterDescription("channels", "It's possible to select the channels "
      "of the output image. There are 3 modes, the available choices are:");

    AddChoice("channels.all", "Default mode");
154 155
    SetParameterDescription("channels.all", 
      "Select all bands in the input image, (1,...,n).");
156 157

    AddChoice("channels.grayscale", "Grayscale mode");
158 159 160 161
    SetParameterDescription("channels.grayscale", 
      "Display single channel as standard color image.");
    AddParameter(ParameterType_Int, "channels.grayscale.channel", 
      "Grayscale channel");
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
    SetDefaultParameterInt("channels.grayscale.channel", 1);
    SetMinimumParameterIntValue("channels.grayscale.channel", 1);

    AddChoice("channels.rgb", "RGB composition");
    SetParameterDescription("channels.rgb", "Select 3 bands in the input image "
      "(multi-bands), by default (1,2,3).");

    AddParameter(ParameterType_Int, "channels.rgb.red", "Red Channel");
    SetParameterDescription("channels.rgb.red", "Red channel index.");
    SetMinimumParameterIntValue("channels.rgb.red", 1);
    AddParameter(ParameterType_Int, "channels.rgb.green", "Green Channel");
    SetParameterDescription("channels.rgb.green", "Green channel index.");
    SetMinimumParameterIntValue("channels.rgb.green", 1);
    AddParameter(ParameterType_Int, "channels.rgb.blue", "Blue Channel");
    SetParameterDescription("channels.rgb.blue", "Blue channel index.");
    SetMinimumParameterIntValue("channels.rgb.blue", 1);

    AddParameter(ParameterType_Float, "outmin", "Output min value");
    SetDefaultParameterFloat("outmin", 0.0);
    SetParameterDescription( "outmin", "Minimum value of the output image." );
    AddParameter(ParameterType_Float, "outmax", "Output max value");
    SetDefaultParameterFloat("outmax", 255.0);
    SetParameterDescription( "outmax", "Maximum value of the output image." );
    MandatoryOff("outmin");
    MandatoryOff("outmax");

    AddRAMParameter();

    // Doc example parameter settings
    SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_XS.tif");
    SetDocExampleParameterValue("out", "otbConvertWithScalingOutput.png");
    SetDocExampleParameterValue("type", "linear");
    SetDocExampleParameterValue("channels", "rgb");
    SetDocExampleParameterValue("outmin", "0");
    SetDocExampleParameterValue("outmax", "255");

    SetOfficialDocLink();
  }

201
  void DoUpdateParameters() override
202 203 204
  {
    // Read information
    if ( HasValue("in") )
205
    {
206 207
      typedef otb::ImageMetadataInterfaceBase ImageMetadataInterfaceType;
      ImageMetadataInterfaceType::Pointer metadataInterface =
208 209
        ImageMetadataInterfaceFactory::CreateIMI(
          GetParameterImage("in")->GetMetaDataDictionary());
210 211 212 213 214 215 216 217 218 219

      int nbBand = GetParameterImage("in")->GetNumberOfComponentsPerPixel();
      SetMaximumParameterIntValue("channels.grayscale.channel", nbBand);
      SetMaximumParameterIntValue("channels.rgb.red", nbBand);
      SetMaximumParameterIntValue("channels.rgb.green", nbBand);
      SetMaximumParameterIntValue("channels.rgb.blue", nbBand);

      if (nbBand > 1)
      {
        // get band index : Red/Green/Blue, in depending on the sensor
220 221 222 223
        auto const& display = metadataInterface->GetDefaultDisplay();
        SetDefaultParameterInt("channels.rgb.red", display[0] + 1);
        SetDefaultParameterInt("channels.rgb.green", display[1] + 1);
        SetDefaultParameterInt("channels.rgb.blue", display[2] + 1);
224
      }
225
    }
226
  }
227

228 229 230 231 232 233 234 235 236 237 238
  template<class TImageType>
  void GenericDoExecute()
  {
    // Clear previously registered filters
    m_Filters.clear();

    std::string rescaleType = this->GetParameterString("type");
    typedef otb::VectorRescaleIntensityImageFilter<FloatVectorImageType, TImageType> RescalerType;
    typename RescalerType::Pointer rescaler = RescalerType::New();

    // selected channel
239
    auto tempImage = GetSelectedChannels<FloatVectorImageType>();
240 241 242 243 244 245

    const unsigned int nbComp(tempImage->GetNumberOfComponentsPerPixel());

    // We need to subsample the input image in order to estimate its histogram
    // Shrink factor is computed so as to load a quicklook of 1000
    // pixels square at most
246
    auto imageSize = tempImage->GetLargestPossibleRegion().GetSize();
247 248
    unsigned int shrinkFactor = std::max({int(imageSize[0])/1000, 
      int(imageSize[1])/1000, 1});
249 250 251 252 253
    otbAppLogDEBUG( << "Shrink factor used to compute Min/Max: "<<shrinkFactor );

    otbAppLogDEBUG( << "Shrink starts..." );
    typename ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();
    shrinkFilter->SetShrinkFactor(shrinkFactor);
254 255 256 257
    shrinkFilter->GetStreamer()->
      SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
    AddProcess(shrinkFilter->GetStreamer(), 
      "Computing shrink Image for min/max estimation...");
258 259 260

    if ( rescaleType == "log2")
    {
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
      // define lambda function that applies a log to all bands of the input pixel
      auto logFunction = [](FloatVectorImageType::PixelType & vectorOut, const FloatVectorImageType::PixelType & vectorIn) {
	assert(vectorOut.Size() == vectorIn.Size() && "Input vector types don't have the same size");
	
	for (unsigned int i = 0; i < vectorIn.Size() ; i++) {
	  vectorOut[i] = std::log(vectorIn[i]);
	}
	
      };
// creates functor filter
      auto transferLogFilter = NewFunctorFilter(logFunction,tempImage->GetNumberOfComponentsPerPixel(),{{0,0}});
      
      // save a reference to the functor
      m_Filters.push_back(transferLogFilter.GetPointer());
      
      transferLogFilter->SetVariadicInputs(tempImage);
      transferLogFilter->UpdateOutputInformation();
      
      shrinkFilter->SetInput(transferLogFilter->GetOutput());
      rescaler->SetInput(transferLogFilter->GetOutput());
281 282 283
      shrinkFilter->Update();
    }
    else
284
    {
285 286 287
      shrinkFilter->SetInput(tempImage);
      rescaler->SetInput(tempImage);
      shrinkFilter->Update();
288 289
    }

290 291
    otbAppLogDEBUG( << "Evaluating input Min/Max..." );
    itk::ImageRegionConstIterator<FloatVectorImageType>
292 293
      it(shrinkFilter->GetOutput(), 
        shrinkFilter->GetOutput()->GetLargestPossibleRegion());
294 295

    typename ListSampleType::Pointer listSample = ListSampleType::New();
296 297
    listSample->SetMeasurementVectorSize( 
      tempImage->GetNumberOfComponentsPerPixel());
298 299 300

    // Now we generate the list of samples
    if (IsParameterEnabled("mask"))
301
    {
302 303
      UInt8ImageType::Pointer mask = this->GetParameterUInt8Image("mask");
      UInt8ShrinkFilterType::Pointer maskShrinkFilter = UInt8ShrinkFilterType::New();
304 305
      maskShrinkFilter->SetShrinkFactor(shrinkFactor);
      maskShrinkFilter->SetInput(mask);
306 307
      maskShrinkFilter->GetStreamer()->
        SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
308 309
      maskShrinkFilter->Update();

310
      auto itMask = itk::ImageRegionConstIterator<UInt8ImageType>(
311
        maskShrinkFilter->GetOutput(),
312
        maskShrinkFilter->GetOutput()->GetLargestPossibleRegion());
313 314 315 316

      // Remove masked pixels
      it.GoToBegin();
      itMask.GoToBegin();
317
      for(; !it.IsAtEnd(); ++it, ++itMask)
318
      {
319 320
        // valid pixels are non zero
        if (itMask.Get() != 0)
321 322 323 324
        {
          listSample->PushBack(it.Get());
        }
      }
325
      // if listSample is empty
326 327
      if (listSample->Size() == 0)
      {
328 329
        otbAppLogINFO( << "All pixels were masked, the application assume "
          "a wrong mask and include all the image");
330
      }
331 332
    }

333 334
    // get all pixels : if mask is disable or all pixels were masked
    if ((!IsParameterEnabled("mask")) || (listSample->Size() == 0))
335
    {
336 337 338 339 340 341 342
      for(it.GoToBegin(); !it.IsAtEnd(); ++it)
      {
        listSample->PushBack(it.Get());
      }
    }

    // And then the histogram
343 344
    typename HistogramsGeneratorType::Pointer histogramsGenerator = 
      HistogramsGeneratorType::New();
345 346 347 348 349 350 351 352 353 354 355
    histogramsGenerator->SetListSample(listSample);
    histogramsGenerator->SetNumberOfBins(255);
    // Samples with nodata values are ignored
    histogramsGenerator->NoDataFlagOn();
    histogramsGenerator->Update();
    auto histOutput = histogramsGenerator->GetOutput();
    assert(histOutput);

    // And extract the lower and upper quantile
    typename FloatVectorImageType::PixelType inputMin(nbComp), inputMax(nbComp);
    for(unsigned int i = 0; i < nbComp; ++i)
356
    {
357 358
      auto && elm = histOutput->GetNthElement(i);
      assert(elm);
359 360 361 362
      inputMin[i] = elm->Quantile(0, 
        0.01 * GetParameterFloat("quantile.low"));
      inputMax[i] = elm->Quantile(0, 
        1.0 - 0.01 * GetParameterFloat("quantile.high"));
363 364
    }

365 366 367
    otbAppLogDEBUG( << std::setprecision(5) 
                    << "Min/Max computation done : min=" 
                    << inputMin
368 369 370 371 372 373 374 375 376 377
                    << " max=" << inputMax );

    rescaler->AutomaticInputMinMaxComputationOff();
    rescaler->SetInputMinimum(inputMin);
    rescaler->SetInputMaximum(inputMax);

    if ( rescaleType == "linear")
    {
      rescaler->SetGamma(GetParameterFloat("type.linear.gamma"));
    }
378

379 380 381
    typename TImageType::PixelType minimum(nbComp);
    typename TImageType::PixelType maximum(nbComp);

382
    /*
383 384
    float outminvalue = std::numeric_limits<typename TImageType::InternalPixelType>::min();
    float outmaxvalue = std::numeric_limits<typename TImageType::InternalPixelType>::max();
385 386 387 388 389 390 391
    // TODO test outmin/outmax values
    if (outminvalue > GetParameterFloat("outmin"))
      itkExceptionMacro("The outmin value at " << GetParameterFloat("outmin") << 
                        " is too low, select a value in "<< outminvalue <<" min.");
    if ( outmaxvalue < GetParameterFloat("outmax") )
      itkExceptionMacro("The outmax value at " << GetParameterFloat("outmax") << 
                        " is too high, select a value in "<< outmaxvalue <<" max.");
392
    */
393

394 395 396 397 398 399 400 401 402 403 404
    maximum.Fill( GetParameterFloat("outmax") );
    minimum.Fill( GetParameterFloat("outmin") );

    rescaler->SetOutputMinimum(minimum);
    rescaler->SetOutputMaximum(maximum);

    m_Filters.push_back(rescaler.GetPointer());
    SetParameterOutputImage<TImageType>("out", rescaler->GetOutput());
  }

  // Get the bands order
405
  std::vector<int> const GetChannels()
406 407
  {
    std::vector<int> channels;
408

409 410
    int nbChan = GetParameterImage("in")->GetNumberOfComponentsPerPixel();
    std::string channelMode = GetParameterString("channels");
411

412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436
    if(channelMode == "grayscale")
    {
      if (GetParameterInt("channels.grayscale.channel") <= nbChan)
      {
        channels = {GetParameterInt("channels.grayscale.channel"),
        GetParameterInt("channels.grayscale.channel"),
        GetParameterInt("channels.grayscale.channel")};
      }
      else
      {
        itkExceptionMacro(<< "The channel has an invalid index");
      }
    }
    else if (channelMode == "rgb")
    {
      if ((GetParameterInt("channels.rgb.red") <= nbChan)
        && ( GetParameterInt("channels.rgb.green") <= nbChan)
        && ( GetParameterInt("channels.rgb.blue")   <= nbChan))
      {
        channels = {GetParameterInt("channels.rgb.red"),
        GetParameterInt("channels.rgb.green"),
        GetParameterInt("channels.rgb.blue")};
      }
      else
      {
437 438
        itkExceptionMacro(<< "At least one needed channel has an invalid "
          "index");
439 440 441 442 443 444 445 446 447 448
      }
    }
    else if (channelMode == "all")
    {
      // take all bands
      channels.resize(nbChan);
      std::iota(channels.begin(), channels.end(), 1);
    }
    return channels;
  }
449

450 451 452 453 454
  // return an image with the bands order modified of the input image
  template<class TImageType>
  typename TImageType::Pointer GetSelectedChannels()
  {
    typedef MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
455 456
      typename TImageType::InternalPixelType> ExtractROIFilterType;
    typedef otb::ImageList<otb::Image<typename TImageType::InternalPixelType> > ImageListType;
457
    typedef ImageListToVectorImageFilter<ImageListType,
458
                                         TImageType > ListConcatenerFilterType;
459

460 461 462
    typename ImageListType::Pointer imageList  = ImageListType::New();
    typename ListConcatenerFilterType::Pointer concatener = 
      ListConcatenerFilterType::New();
463

464 465
    //m_Filters.push_back(imageList.GetPointer());
    m_Filters.push_back(concatener.GetPointer());
466

467 468 469
    const bool monoChannel = IsParameterEnabled("channels.grayscale");

    // get band order
470
    const std::vector<int> channels = GetChannels();
471 472 473

    for (auto && channel : channels)
    {
474 475
      typename ExtractROIFilterType::Pointer extractROIFilter = 
        ExtractROIFilterType::New();
476 477
      m_Filters.push_back(extractROIFilter.GetPointer());
      extractROIFilter->SetInput(GetParameterImage("in"));
478 479 480
      if (!monoChannel) 
        extractROIFilter->SetChannel(channel);

481 482
      extractROIFilter->UpdateOutputInformation();
      imageList->PushBack(extractROIFilter->GetOutput());
483 484
    }

485 486 487 488 489
    concatener->SetInput(imageList);
    concatener->UpdateOutputInformation();

    return concatener->GetOutput();
  }
490

491

492
  void DoExecute() override
493 494
  {
    switch ( this->GetParameterOutputImagePixelType("out") )
495
    {
496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517
      case ImagePixelType_uint8:
        GenericDoExecute<UInt8VectorImageType>();
        break;
      case ImagePixelType_int16:
        GenericDoExecute<Int16VectorImageType>();
        break;
      case ImagePixelType_uint16:
        GenericDoExecute<UInt16VectorImageType>();
        break;
      case ImagePixelType_int32:
        GenericDoExecute<Int32VectorImageType>();
        break;
      case ImagePixelType_uint32:
        GenericDoExecute<UInt32VectorImageType>();
        break;
      case ImagePixelType_float:
        GenericDoExecute<FloatVectorImageType>();
        break;
      case ImagePixelType_double:
        GenericDoExecute<DoubleVectorImageType>();
        break;
      default:
518 519 520
        itkExceptionMacro("Unknown pixel type " << this->GetParameterOutputImagePixelType("out") <<"." << std::endl
                          << "The DynamicConvert application does not support complex pixel type as output." << std::endl
                          << "You can use instead the ExtractROI application to perform complex image conversion.");
521
        break;
522
    }
523
  }
524

525 526
  std::vector<itk::LightObject::Pointer> m_Filters;
};
527 528 529 530 531 532

}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::DynamicConvert)