otbContrastEnhancement.cxx 35.6 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3 4 5 6 7 8 9 10 11 12 13
 *
 * 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
 *
Yannick TANGUY's avatar
Yannick TANGUY committed
14
 * Unless required by applicable law or agreed to in writing, software
15 16 17 18 19
 * 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.
 */
Antoine Regimbeau's avatar
Antoine Regimbeau committed
20

21 22 23
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

24 25
#include "otbVectorImageToImageListFilter.h"
#include "otbImageListToVectorImageFilter.h"
Antoine Regimbeau's avatar
Antoine Regimbeau committed
26
#include "otbStreamingStatisticsVectorImageFilter.h"
Antoine Regimbeau's avatar
Antoine Regimbeau committed
27
#include "otbStreamingStatisticsImageFilter.h"
28
#include "otbFunctorImageFilter.h"
29
#include "itkStreamingImageFilter.h"
30
#include "otbInPlacePassFilter.h"
31

32
#include "otbComputeHistoFilter.h"
33
#include "otbComputeGainLutFilter.h"
34
#include "otbApplyGainFilter.h"
Antoine Regimbeau's avatar
Antoine Regimbeau committed
35
#include "otbImageFileWriter.h"
36
#include "itkImageRegionIterator.h"
37
#include <string>
38
#include "otbStreamingHistogramVectorImageFilter.h"
39

40 41
namespace otb
{
42

43 44 45
namespace Wrapper
{

46 47 48
namespace Functor
{

Cédric Traizet's avatar
Cédric Traizet committed
49
template <class TInput, class TOutput>
50 51 52
class LuminanceOperator
{
public:
53
  LuminanceOperator() = default;
Cédric Traizet's avatar
Cédric Traizet committed
54 55

  size_t OutputSize(const std::array<size_t, 1>&) const
56 57 58
  {
    return 1;
  }
59
  virtual ~LuminanceOperator() = default;
60

Cédric Traizet's avatar
Cédric Traizet committed
61 62 63 64 65
  TOutput operator()(TInput input)
  {
    TOutput out(1);
    out[0] = m_LumCoef[0] * input[m_Rgb[0]] + m_LumCoef[1] * input[m_Rgb[1]] + m_LumCoef[2] * input[m_Rgb[2]];
    return out;
66 67 68
  } // end operator ()


69
  void SetRgb( std::vector<unsigned int> rgb)
70 71 72 73
    {
    m_Rgb = rgb;
    }

Antoine Regimbeau's avatar
Antoine Regimbeau committed
74
  void SetLumCoef(std::vector<float> lumCoef)
75 76 77
    {
    m_LumCoef = lumCoef;
    }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
78 79 80 81
  std::vector<float> GetLumCoef()
    {
    return m_LumCoef;
    }
82
private:
83
  std::vector<unsigned int> m_Rgb;
Antoine Regimbeau's avatar
Antoine Regimbeau committed
84
  std::vector<float> m_LumCoef;
85
}; // end of functor class LuminanceOperator
86

Cédric Traizet's avatar
Cédric Traizet committed
87
} // namespace functor
88

89 90 91 92 93 94 95 96 97
class ContrastEnhancement : public Application
{
public:
  /** Standard class typedefs. */
  typedef ContrastEnhancement	              Self;
  typedef Application	                      Superclass;
  typedef itk::SmartPointer < Self >	      Pointer;
  typedef itk::SmartPointer < const Self >	ConstPointer;

98
  typedef otb::VectorImage < unsigned int , 2 > HistogramType;
99 100 101 102
  typedef otb::VectorImage < double , 2 > LutType;

  typedef FloatImageType::PixelType ImagePixelType;

103
  typedef otb::ComputeHistoFilter < FloatImageType , 
104
                                    HistogramType > 
105
          HistoFilterType;
106 107
  typedef otb::ComputeGainLutFilter < HistogramType , 
                                      LutType > 
108
          GainLutFilterType;
109
  typedef otb::ApplyGainFilter < FloatImageType , 
110 111 112 113 114 115
                                 LutType , FloatImageType > 
          ApplyFilterType;
  typedef otb::ImageList < FloatImageType > ImageListType;

  typedef otb::VectorImageToImageListFilter < FloatVectorImageType, 
                                              ImageListType > 
116
          VectorToImageListFilterType;
117 118 119

  typedef otb::ImageListToVectorImageFilter < ImageListType, 
                                              FloatVectorImageType > 
120
          ImageListToVectorFilterType;
121 122

  typedef otb::StreamingStatisticsVectorImageFilter < FloatVectorImageType >
Antoine Regimbeau's avatar
Antoine Regimbeau committed
123
          VectorStatsFilterType;
124 125

  typedef otb::StreamingStatisticsImageFilter < FloatImageType >
126
          StatsFilterType;
127

Cédric Traizet's avatar
Cédric Traizet committed
128
  typedef otb::FunctorImageFilter<Functor::LuminanceOperator<FloatVectorImageType::PixelType, FloatVectorImageType::PixelType>> LuminanceFunctorType;
129 130

  typedef itk::StreamingImageFilter < LutType , LutType > 
131
          StreamingImageFilterType;
132

133
  typedef otb::InPlacePassFilter < FloatImageType > BufferFilterType;
134

135 136
  typedef otb::StreamingHistogramVectorImageFilter < FloatVectorImageType > 
      HistoPersistentFilterType;
137

138 139 140 141 142 143
  /** Standard macro */
  itkNewMacro( Self );
 
  itkTypeMacro( ContrastEnhancement , otb::Application );

private:
144

145
	void DoInit() override
146
	{
147 148 149
		SetName("ContrastEnhancement");
    SetDescription("This application is the implementation of the histogram "
      "equalization algorithm. It can be used to enhance contrast in an image "
Guillaume Pasero's avatar
Guillaume Pasero committed
150
      "or to reduce the dynamic of the image without losing too much contrast. "
151
      "It offers several options as a nodata value, "
152
      "a contrast limitation factor, a local version of the algorithm and "
153
      "also a mode to equalize the luminance of the image.");
154 155

    // Documentation
156 157 158
    SetDocLongDescription("This application is the implementation of the "
      "histogram equalization algorithm. The idea of the algorithm is to use "
      "the whole available dynamic. In order to do so it computes a histogram "
Antoine Regimbeau's avatar
Antoine Regimbeau committed
159
      "over the image and then use the whole dynamic: meaning flattening the "
160 161
      "histogram. That gives us gain for each bin that transform the original "
      "histogram into the flat one. This gain is then apply on the original "
162 163 164
      "image.\n\n"
      "The application proposes several options to allow a finer result:\n\n"
      "* There is an option to limit contrast. We choose to limit the contrast "
165
      "by modifying the original histogram. To do so, we clip the histogram at a "
166
      "given height and redistribute equally among the bins the clipped population. "
167 168
      "Then we add a local version of the algorithm.\n"
      "* It is possible to apply the algorithm on tiles of the image, instead "
169
      "of on the whole image. That gives us gain depending on "
170 171
      "the value of the pixel and its position in the image. In order to "
      "smoothen the result we interpolate the gain between tiles.");
172 173
    SetDocLimitations("None");
    SetDocAuthors("OTB-Team");
174
    SetDocSeeAlso(" ");
175 176 177 178 179

    AddDocTag(Tags::Filter);

    AddParameter(ParameterType_InputImage,  "in",   "Input Image");
    SetParameterDescription("in", "Input image.");
180
    
181 182
    AddParameter(ParameterType_OutputImage, "out",  "Output Image");
    SetParameterDescription("out", "Output image.");
183
    
184
    AddParameter(ParameterType_Int , "bins" , "Number of bins");
185 186
    SetDefaultParameterInt("bins", 256);
    SetParameterDescription("bins",
187
      "Number of bins in the histogram");
188

189
    AddParameter(ParameterType_Float , "hfact" , "Contrast Limitation");  
190
    SetParameterDescription("hfact","This parameter will set the maximum "
191
      "height accepted in a bin on the input image histogram. "
192
      "The maximum height will be computed as hfact*eqHeight where eqHeight "
193
      "is the height of the theoretical flat histogram. The higher hfact, the "
194
      "higher the contrast."
Antoine Regimbeau's avatar
Antoine Regimbeau committed
195
      "\nWhen using 'luminance mode', it is recommended to limit this factor to a small value (ex: 4)");
196 197
    MandatoryOff("hfact");

198 199
    AddParameter(ParameterType_Float , "nodata" , "Nodata Value");
    SetParameterDescription("nodata","If there is a value in the "
200 201
      "image that has no visualization meaning, it can be ignored by the "
      "algorithm.");
202
    MandatoryOff("nodata");
203

204 205 206 207
    AddParameter(ParameterType_Choice , "spatial" , "Spatial parameters "
      "for the histogram computation");
    AddChoice( "spatial.local" , "Local" );
    SetParameterDescription("spatial.local" , "The histograms will be "
208
      "computed on each thumbnail. Each of the histogram will be "
209
      "equalized and the corresponding gain will be interpolated.");
210 211
    AddChoice( "spatial.global" , "Global" );
    SetParameterDescription("spatial.global" , "The histogram will be "
212 213
      "computed on the whole image. The equalization will be computed on "
      "this histogram.");
214

215
 
216
    AddParameter(ParameterType_Int,"spatial.local.h" , 
217 218
      "Thumbnail height");
    SetParameterDescription("spatial.local.h","Height of the thumbnail "
Antoine Regimbeau's avatar
Antoine Regimbeau committed
219
      "over which the histogram will be computed. The value is in pixels.");
220
    AddParameter(ParameterType_Int,"spatial.local.w" , 
221 222
      "Thumbnail width");
    SetParameterDescription("spatial.local.w","Width of the thumbnail "
Antoine Regimbeau's avatar
Antoine Regimbeau committed
223
      "over which the histogram will be computed. The value is in pixels.");
224

225
    AddParameter(ParameterType_Choice , "minmax" , "Minimum and maximum "
226
      "settings");
227
    SetParameterDescription("minmax","Minimum and maximum value that will "
228 229
      "bound the histogram and thus the dynamic of the resulting image. Values "
      "over those boundaries will be clipped.");
230
    AddChoice( "minmax.auto" , "Automatic" );
231
    SetParameterDescription("minmax.auto" , "Minimum and maximum value will "
232 233
      "be computed on the image (nodata value won't be taken "
      "into account) . Each band will have a minimum and a maximum.");
234
    AddParameter(ParameterType_Bool, "minmax.auto.global", "Global");
235
    SetParameterDescription("minmax.auto.global" , 
236
      "Min/max computation will result in the same minimum and maximum for "
237 238 239
      "all the bands.");
    AddChoice( "minmax.manual" , "Manual settings for min/max values" );
    SetParameterDescription("minmax.manual","Minimum and maximum value will be "
240
      "set by the user");
241 242
    AddParameter(ParameterType_Float , "minmax.manual.min" , "Minimum value");
    AddParameter(ParameterType_Float , "minmax.manual.max" , "Maximum value");
243 244
    // SetDefaultParameterFloat("minmax.manual.min", 0 );
    // SetDefaultParameterFloat("minmax.manual.max", 255 );
245 246
    MandatoryOff("minmax.manual.min");
    MandatoryOff("minmax.manual.max");
247

248
    AddParameter(ParameterType_Choice , "mode" , "What to equalized");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
249
    AddChoice( "mode.each" , "Channels" );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
250
    SetParameterDescription( "mode.each" ,
251
      "Each channel is equalized independently" );
252
    AddChoice( "mode.lum" , "Luminance" );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
253
    SetParameterDescription( "mode.lum" ,
Antoine Regimbeau's avatar
Antoine Regimbeau committed
254 255 256
      "The relative luminance is computed according to the coefficients."
      "Then the histogram is equalized and the gain is applied to each of the "
      "channels. The channel gain will depend on "
257 258 259 260 261
      "the weight (coef) of the channel in the luminance." 
      "\nNote that default values come from color space theories "
      "on how human eyes perceive colors)" 

);
262 263
    AddParameter(ParameterType_Group , "mode.lum.red" , "Red channel" );
    AddParameter(ParameterType_Int , "mode.lum.red.ch" , "Red channel" );
264
    SetDefaultParameterInt("mode.lum.red.ch", 0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
265
    AddParameter(ParameterType_Float , "mode.lum.red.coef" ,
266
      "Value for luminance computation for the red channel" );
267 268
    SetDefaultParameterFloat("mode.lum.red.coef", 0.21 );

269 270
    AddParameter(ParameterType_Group , "mode.lum.green" , "Green channel" );
    AddParameter(ParameterType_Int , "mode.lum.green.ch" , "Green channel" );
271 272
    SetDefaultParameterInt("mode.lum.green.ch", 1 );
    AddParameter(ParameterType_Float , "mode.lum.green.coef" ,
273
      "Value for luminance computation of the green channel" );
274
    SetDefaultParameterFloat("mode.lum.green.coef", 0.71 );
275

276 277
    AddParameter(ParameterType_Group , "mode.lum.blue" , "Blue channel" );
    AddParameter(ParameterType_Int , "mode.lum.blue.ch" , "Blue channel" );
278 279
    SetDefaultParameterInt("mode.lum.blue.ch", 2 );
    AddParameter(ParameterType_Float , "mode.lum.blue.coef" ,
280
      "Value for luminance computation of the blue channel" );
281
    SetDefaultParameterFloat("mode.lum.blue.coef", 0.08 );
282

283 284 285
    SetDefaultParameterInt( "spatial.local.w" , 256 );
    SetDefaultParameterInt( "spatial.local.h" , 256 );

Antoine Regimbeau's avatar
Antoine Regimbeau committed
286
    SetMinimumParameterIntValue("mode.lum.red.ch", 0);
287 288
    SetMinimumParameterIntValue("mode.lum.green.ch", 0);
    SetMinimumParameterIntValue("mode.lum.blue.ch", 0);
289 290 291
    SetMinimumParameterIntValue("bins", 2);
    SetMinimumParameterIntValue("spatial.local.h", 1);
    SetMinimumParameterIntValue("spatial.local.w", 1);
Antoine Regimbeau's avatar
Antoine Regimbeau committed
292

293
    SetExampleComment( "Local contrast enhancement by luminance" , 0 );
294 295
    SetDocExampleParameterValue( "in" , "colours.tif" );
    SetDocExampleParameterValue( "out" , "equalizedcolors.tif float" );
296 297 298 299
    SetDocExampleParameterValue( "bins" , "256" );
    SetDocExampleParameterValue( "spatial.local.w" , "500" );
    SetDocExampleParameterValue( "spatial.local.h" , "500");
    SetDocExampleParameterValue( "mode" , "lum" );
300

301
    AddRAMParameter(); 
302
  }
303

304
  void DoUpdateParameters() override
305
  {
306
    if ( HasValue("in") )
307
      {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
308
      FloatVectorImageType * inImage = GetParameterImage("in");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
309 310 311
      FloatVectorImageType::RegionType::SizeType size;
      size = inImage->GetLargestPossibleRegion().GetSize() ;

312 313 314
      if ( GetParameterString("spatial") == "local" &&
           HasValue("spatial.local.h") && HasValue("spatial.local.w") &&
           HasValue("bins") )
315
        CheckValidity();
Antoine Regimbeau's avatar
Antoine Regimbeau committed
316 317
      
      
318
      if ( !HasUserValue("nodata") && IsParameterEnabled("nodata"))
Antoine Regimbeau's avatar
Antoine Regimbeau committed
319 320 321 322
        SetDefaultValue( inImage , "NODATA" );

      if ( GetParameterString( "mode" ) == "lum" && 
           !HasUserValue("mode.lum.red.ch") &&
323 324
           !HasUserValue("mode.lum.green.ch") &&
           !HasUserValue("mode.lum.blue.ch") )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
325 326 327
        SetDefaultValue( inImage , "RGB" );
      }

328
    if ( GetParameterString("minmax") == "manual" )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
329
      {
330 331
      MandatoryOn("minmax.manual.min");
      MandatoryOn("minmax.manual.max");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
332 333 334
      }
    else if ( GetParameterString("minmax") == "auto" )
      {
335 336
      MandatoryOff("minmax.manual.min");
      MandatoryOff("minmax.manual.max");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
337 338 339
      }
  }

340
  void DoExecute() override
Antoine Regimbeau's avatar
Antoine Regimbeau committed
341
  {
342 343 344
    m_MinMaxMode = GetParameterString("minmax");
    m_EqMode = GetParameterString("mode");
    m_SpatialMode = GetParameterString("spatial");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
345
    FloatVectorImageType * inImage = GetParameterImage("in");
346 347 348
    WarningGlobalOrNot( inImage );
    WarningMinMax();
    LogInfo();
Antoine Regimbeau's avatar
Antoine Regimbeau committed
349
    ImageListType::Pointer outputImageList ( ImageListType::New() ); 
350 351 352
    m_VectorToImageListFilter = VectorToImageListFilterType::New() ;
    m_VectorToImageListFilter->SetInput( inImage );
    m_VectorToImageListFilter->UpdateOutputInformation();
Antoine Regimbeau's avatar
Antoine Regimbeau committed
353
    ImageListType::Pointer inputImageList = 
354
                  m_VectorToImageListFilter->GetOutput();
355
    unsigned int nbChannel = inImage->GetVectorLength ();
Antoine Regimbeau's avatar
Antoine Regimbeau committed
356

357
    if ( m_EqMode == "each")
Antoine Regimbeau's avatar
Antoine Regimbeau committed
358
      {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
359
      // Each channel will be equalized
360 361 362 363
      m_GainLutFilter.resize(nbChannel);
      m_ApplyFilter.resize(nbChannel);
      m_BufferFilter.resize(nbChannel);
      m_StreamingFilter.resize(nbChannel);
Antoine Regimbeau's avatar
Antoine Regimbeau committed
364
      PerBandEqualization( inImage , inputImageList , 
Antoine Regimbeau's avatar
Antoine Regimbeau committed
365
                           nbChannel , outputImageList );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
366
      }
367
    else if ( m_EqMode == "lum")
Antoine Regimbeau's avatar
Antoine Regimbeau committed
368
      {
369
      std::vector< unsigned int > rgb( 3 , 0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
370
      rgb[0] = GetParameterInt("mode.lum.red.ch");
371 372
      rgb[1] = GetParameterInt("mode.lum.green.ch");
      rgb[2] = GetParameterInt("mode.lum.blue.ch");
373 374 375 376 377 378 379
      if ( !( nbChannel > std::max( rgb[0] , std::max( rgb[1] , rgb[2] ) ) ) )
        {
        std::ostringstream oss;
        oss<<"One of the selected channel needed for luminance computation "
        "exceed the number of component of the image.";
        otbAppLogFATAL( << oss.str() )
        }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
380 381 382 383
      ComputeLuminance( inImage , rgb );
      LuminanceEqualization( inputImageList , rgb , outputImageList );
      }

384 385
    m_ImageListToVectorFilterOut = ImageListToVectorFilterType::New() ;
    m_ImageListToVectorFilterOut->SetInput(outputImageList);
386 387
    SetParameterOutputImage( "out" , 
        m_ImageListToVectorFilterOut->GetOutput() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
  }

  // Look for default values in the image metadata
  void SetDefaultValue( const FloatVectorImageType * inImage ,
                        std::string what)
  {
    typedef ImageMetadataInterfaceBase ImageMetadataInterfaceType;
    ImageMetadataInterfaceType::Pointer metadataInterface = 
          ImageMetadataInterfaceFactory
                ::CreateIMI(inImage->GetMetaDataDictionary());
    if ( what == "NODATA" )
      {
      std::vector<double> values;
      std::vector<bool> flags;

      bool ret = metadataInterface->GetNoDataFlags(flags,values);

      if(ret && !values.empty() && !flags.empty() && flags[0])
406
        {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
407
        SetParameterFloat( "nodata" , static_cast<float>( values[0] ) );
408
        }
409 410 411 412
      else
        {
        SetParameterFloat( "nodata" , 0 );
        }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
413 414
      }
    else if ( what == "RGB" )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
415
      {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
416 417 418 419
      std::vector<unsigned int> rgb = 
                    metadataInterface->GetDefaultDisplay() ;
      unsigned int m = inImage->GetVectorLength ();
      SetParameterInt( "mode.lum.red.ch" , rgb[0] );
420 421
      SetParameterInt( "mode.lum.green.ch" , rgb[1] );
      SetParameterInt( "mode.lum.blue.ch" , rgb[2] );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
422
      if( m < rgb[ 0 ] )
423
        {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
424 425
        SetParameterFloat ("mode.lum.red.coef" , 0.0 );
        SetParameterInt( "mode.lum.red.ch" , 0 );
426
        }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
427 428
      if( m < rgb[ 1 ] )
        {
429
        SetParameterFloat ("mode.lum.green.coef" , 0.0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
430
        SetParameterInt( "mode.lum.gre.ch" , 0 );
431
        }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
432
      if( m < rgb[ 2 ] )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
433
        {
434 435
        SetParameterFloat ("mode.lum.blue.coef" , 0.0 );
        SetParameterInt( "mode.lum.blue.ch" , 0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
436 437
        }
      }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
438
  }
439

440 441 442 443
  // Log info for the user
  void LogInfo()
  {
    std::ostringstream oss;
444
    oss << "The application has been launched with the following parameters :"
445
    <<std::endl;
446 447 448 449 450 451 452 453 454 455
    oss << "- number of bins : " << GetParameterInt("bins") << std::endl;
    if ( HasValue("hfact") )
      {
      oss << "- contrast limtaition factor : "
      << GetParameterFloat("hfact") << std::endl;
      }
    else
      {
      oss << "- no contrast limitation factor" << std::endl;
      }
456 457 458
    oss << "- spatial parameters : " << m_SpatialMode ;
    if ( m_SpatialMode == "local" )
      {
459
      oss<< " with a thumbnail of " << m_ThumbSize[0] <<" X "<< m_ThumbSize[1] ;
460 461 462 463 464 465 466 467 468 469 470 471 472 473
      }
    oss << std::endl << "- equalisation of ";
    if ( m_EqMode == "each" )
      {
      oss << "each channel";
      }
    else
      {
      oss << "the luminance";
      }
    oss << std::endl << "- Min/Max parameters : ";
    if ( m_MinMaxMode == "auto" )
      {
      oss << "automatic";
474
      if ( GetParameterInt( "minmax.auto.global" ) )
475 476 477 478 479 480
        { 
        oss << " and global";
        }
      }
    else 
      {
481 482
      oss << GetParameterFloat("minmax.manual.min") << "/" << 
        GetParameterFloat("minmax.manual.max");
483 484 485 486 487
      }

    otbAppLogINFO( << oss.str() );
  }

488
  // Force global computation if the thumbsize is equal to the largest region
489
  void WarningGlobalOrNot( FloatVectorImageType * input)
490 491
  {
    auto size = input->GetLargestPossibleRegion().GetSize();
492 493 494 495 496 497 498 499 500 501
    if ( m_SpatialMode == "global" )
      {
      m_ThumbSize[0] = size[0];
      m_ThumbSize[1] = size[1];
      }  
    else 
      {
      m_ThumbSize[0] = GetParameterInt("spatial.local.w");
      m_ThumbSize[1] = GetParameterInt("spatial.local.h");
      if ( size[0] == m_ThumbSize[0] && size[1] == m_ThumbSize[1] )
502 503 504 505 506 507 508 509
        {
        std::ostringstream oss;
        oss<<"Warning you choose to compute the histogram with a local "
        "method whereas you have selected the whole image for the thumbnail "
        "size. In order to use less memory consider using the global option for "
        "histogram computation.";
        otbAppLogWARNING( << oss.str() );
        }
510 511 512 513 514 515
      }
  }

  // Check for min max validity 
  void WarningMinMax()
  {
516 517 518
    if ( m_MinMaxMode == "manual" &&  
         GetParameterFloat( "minmax.manual.min" ) > 
         GetParameterFloat( "minmax.manual.max" ) )
519 520
      {
      std::ostringstream oss;
521
      oss<<"The minimum (" << GetParameterFloat( "minmax.manual.min" ) <<
522
      ") is superior to the maximum (" 
523
      << GetParameterFloat( "minmax.manual.max" )
524 525 526 527
      << ") please correct this error or allow the application to compute "
      "those parameters";
      otbAppLogFATAL( << oss.str() )
      }
528 529
  }

Antoine Regimbeau's avatar
Antoine Regimbeau committed
530
  // Check if the image size is a multiple of the thumbnail size
531
  void CheckValidity()
Antoine Regimbeau's avatar
Antoine Regimbeau committed
532 533
  {
    std::ostringstream oss;
534 535 536
    long nbPixel = GetParameterInt("spatial.local.w") * 
                   GetParameterInt("spatial.local.h");
    if ( nbPixel < 10 * GetParameterInt("bins"))
537 538 539
      {   
      oss<<"Warning in parameters selection the thumbnail size is small "
      "in comparison with the number of bin. Histogram may not have much sens. "
540
      "For better result enlarge thumbnail size or reduce number of bin.";
541
      otbAppLogINFO( << oss.str() );
542
      }
543 544
  }

Antoine Regimbeau's avatar
Antoine Regimbeau committed
545 546 547 548
  // Compute min max from a vector image
  void ComputeVectorMinMax( const FloatVectorImageType::Pointer inImage ,
                            FloatVectorImageType::PixelType & max ,
                            FloatVectorImageType::PixelType & min )
549
  {
550
    if ( m_MinMaxMode == "manual" )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
551
      {
552 553
      min.Fill( GetParameterFloat("minmax.manual.min") );
      max.Fill( GetParameterFloat("minmax.manual.max") );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
554
      }
555
    else
Antoine Regimbeau's avatar
Antoine Regimbeau committed
556
      {
557 558
      VectorStatsFilterType::Pointer 
                statFilter ( VectorStatsFilterType::New() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
559
      statFilter->SetIgnoreInfiniteValues(true);
560
      if( IsParameterEnabled("nodata") )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
561 562 563 564 565
        {
        statFilter->SetIgnoreUserDefinedValue(true);
        statFilter->SetUserIgnoredValue( GetParameterFloat("nodata") );
        }
      statFilter->SetInput( inImage );
566
      AddProcess(statFilter->GetStreamer(), "Computing statistics");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
567 568 569
      statFilter->Update();
      min = statFilter->GetMinimum();
      max = statFilter->GetMaximum();
570
      if ( GetParameterInt("minmax.auto.global") )
571 572
        {
        float temp(min[0]);
573
        for ( unsigned int i = 1 ; i < min.GetSize() ; i++ )
574 575 576 577 578
          {
          temp = std::min(temp , min[i]);
          }
        min.Fill(temp);
        temp = max[0];
579
        for ( unsigned int i = 1 ; i < max.GetSize() ; i++ )
580 581 582 583 584
          {
          temp = std::max(temp , max[i]);
          }
        max.Fill(temp);
        }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
585
      }
586 587
    std::ostringstream oss;
    oss<<"Minimum and maximum are for each channel : ";
588
    if ( GetParameterInt("minmax.auto.global") || 
589
          m_MinMaxMode == "manual" )
590 591 592 593 594 595 596 597 598 599 600
      {
      oss<<std::endl<<min[0]<<" and "<<max[0];
      }
    else 
      {
      for ( unsigned int i = 0 ; i < min.GetSize() ; i++ )
        {
          oss<<std::endl<<min[i]<<" and "<<max[i];
        }
      }
    otbAppLogINFO( << oss.str() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
601
  }
Antoine Regimbeau's avatar
Antoine Regimbeau committed
602

603 604 605 606 607 608 609 610 611 612 613
  // Prepare the first half of the pipe that is common to every methode of 
  // equalization
  void SetUpPipeline( unsigned int channel ,
                      const FloatImageType::Pointer input )
  {
    m_GainLutFilter[channel] = GainLutFilterType::New();
    m_ApplyFilter[channel] = ApplyFilterType::New();
    m_StreamingFilter[channel] = StreamingImageFilterType::New();
    m_BufferFilter[channel] = BufferFilterType::New();
    m_BufferFilter[channel]->SetInput( input );
    m_GainLutFilter[channel]->SetInput ( m_Histogram[channel] );
614 615
    m_StreamingFilter[channel]->SetInput( 
      m_GainLutFilter[channel]->GetOutput() );
616 617 618 619 620 621
    m_ApplyFilter[channel]->SetInputImage ( 
      m_BufferFilter[channel]->GetOutput() );
    m_ApplyFilter[channel]->SetInputLut( 
      m_StreamingFilter[channel]->GetOutput() );
  }

Antoine Regimbeau's avatar
Antoine Regimbeau committed
622
  // Function corresponding to the "each" mode
623
  void PerBandEqualization( const FloatVectorImageType::Pointer inImage ,
Antoine Regimbeau's avatar
Antoine Regimbeau committed
624
                            const ImageListType::Pointer inputImageList ,
625
                            const unsigned int nbChannel,
Antoine Regimbeau's avatar
Antoine Regimbeau committed
626 627
                            ImageListType::Pointer outputImageList )
  {
Antoine Regimbeau's avatar
Antoine Regimbeau committed
628
    FloatVectorImageType::PixelType min(nbChannel) , max(nbChannel);
Antoine Regimbeau's avatar
Antoine Regimbeau committed
629 630 631
    min.Fill(0);
    max.Fill(0);
    ComputeVectorMinMax( inImage , max , min );
632

633
    if ( m_SpatialMode == "global" )
634
      PersistentComputation( inImage , nbChannel , max , min );
635 636
    else
      {
637 638 639 640 641 642
      float thresh (-1);
      if ( HasValue("hfact") )
        {
        thresh = GetParameterInt("hfact");
        }

643 644 645 646 647 648 649
      m_HistoFilter.resize(nbChannel);
      m_Histogram.resize(nbChannel);
      for (unsigned int channel = 0 ; channel < nbChannel ; channel++)
        {
        m_HistoFilter[channel] = HistoFilterType::New();
        m_HistoFilter[channel]->SetInput( 
            inputImageList->GetNthElement(channel) );
650 651 652
        SetHistoFilterParameter( m_HistoFilter[channel] ,
                                 min[channel] ,
                                 max[channel] ,
653
                                 GetParameterInt("bins") ,
654
                                 thresh );
655 656 657 658
        m_Histogram[channel] = m_HistoFilter[channel]->GetHistoOutput();
        }
      }

659
    for ( unsigned int channel = 0 ; channel < nbChannel ; channel++ ) 
660
      {
661
      SetUpPipeline( channel , inputImageList->GetNthElement(channel) );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
662

Antoine Regimbeau's avatar
Antoine Regimbeau committed
663
      if ( min[channel] == max[channel] )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
664
        {
665 666 667 668
          std::ostringstream oss;
          oss<< "Channel " <<channel << " is constant : " << "min = " <<
          min[channel] << " and max = " << max[channel] ;
          otbAppLogINFO( << oss.str() );
669 670
          m_BufferFilter[channel]->SetInput( 
                inputImageList->GetNthElement(channel) );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
671
          outputImageList->PushBack( m_BufferFilter[channel]->GetOutput() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
672
          continue;
Antoine Regimbeau's avatar
Antoine Regimbeau committed
673
        }
674

675 676 677
      SetGainLutFilterParameter( m_GainLutFilter[channel] ,
                                 min[channel] ,
                                 max[channel]);
678 679 680
      SetApplyFilterParameter( m_ApplyFilter[channel] ,
                               min[channel] ,
                               max[channel]);
681 682

      outputImageList->PushBack( m_ApplyFilter[channel]->GetOutput() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
683 684 685 686 687
      }
  }

  // Compute the luminance with user parameters
  void ComputeLuminance( const FloatVectorImageType::Pointer inImage ,
688
                         std::vector < unsigned int > rgb )
Antoine Regimbeau's avatar
Antoine Regimbeau committed
689
  {
Guillaume Pasero's avatar
Guillaume Pasero committed
690
    // Retrieve coeffs for each channel
Antoine Regimbeau's avatar
Antoine Regimbeau committed
691
    std::vector < float > lumCoef( 3 , 0.0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
692
    lumCoef[0] = GetParameterFloat("mode.lum.red.coef");
693 694
    lumCoef[1] = GetParameterFloat("mode.lum.green.coef");
    lumCoef[2] = GetParameterFloat("mode.lum.blue.coef");
Antoine Regimbeau's avatar
Antoine Regimbeau committed
695
    // Normalize those coeffs
696
    float sum = std::accumulate( lumCoef.begin() , lumCoef.end() , 0.0 );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
697 698 699 700 701
    assert(sum>0);
    for (int i = 0 ; i<3 ; i++ )
      {
      lumCoef[i] /= sum;
      }
702
    m_LuminanceFunctor =  LuminanceFunctorType::New() ;
Cédric Traizet's avatar
Cédric Traizet committed
703 704
    m_LuminanceFunctor->GetModifiableFunctor().SetRgb(rgb);
    m_LuminanceFunctor->GetModifiableFunctor().SetLumCoef(lumCoef);
705
    m_LuminanceFunctor->SetInput( inImage );
706
    m_LuminanceFunctor->UpdateOutputInformation();
Antoine Regimbeau's avatar
Antoine Regimbeau committed
707 708
  }

709
  // Equalize the luminance and apply the corresponding gain on each channel
Antoine Regimbeau's avatar
Antoine Regimbeau committed
710 711
  // used to compute this luminance
  void LuminanceEqualization( const ImageListType::Pointer inputImageList ,
712
                              const std::vector < unsigned int > rgb ,
Antoine Regimbeau's avatar
Antoine Regimbeau committed
713 714
                              ImageListType::Pointer outputImageList )
  {
715
    m_GainLutFilter.resize( 1 , GainLutFilterType::New() );
716 717
    m_HistoFilter.resize( 1 , HistoFilterType::New() );
    m_StreamingFilter.resize( 1 , StreamingImageFilterType::New() );
718
    m_ApplyFilter.resize(3);
719
    m_BufferFilter.resize(3);
720 721 722
    FloatVectorImageType::PixelType min(1) , max(1);
    ComputeVectorMinMax( m_LuminanceFunctor->GetOutput() , max , min );

723
    if ( m_SpatialMode == "global" )
724 725 726
      PersistentComputation( m_LuminanceFunctor->GetOutput() , 1 , max , min );
    else
      {
727 728 729 730 731
        float thresh (-1);
        if ( HasValue("hfact") )
          {
          thresh = GetParameterInt("hfact");
          }
732
        m_Histogram.resize(1);
733
        m_HistoFilter[0] = HistoFilterType::New();
734 735 736
        m_LuminanceToImageListFilter = VectorToImageListFilterType::New();
        m_LuminanceToImageListFilter->SetInput( 
          m_LuminanceFunctor->GetOutput() );
737 738 739
        SetHistoFilterParameter( m_HistoFilter[0] ,
                                 min[0] ,
                                 max[0] ,
740
                                 GetParameterInt("bins") ,
741
                                 thresh ); 
742 743 744 745 746
        m_LuminanceToImageListFilter->UpdateOutputInformation();
        m_HistoFilter[0]->SetInput( 
          m_LuminanceToImageListFilter->GetOutput()->GetNthElement(0) ) ;
        m_Histogram[0] = m_HistoFilter[0]->GetHistoOutput();
      }
747

748 749
    m_GainLutFilter[0]->SetInput ( m_Histogram[0] );
    m_StreamingFilter[0]->SetInput( m_GainLutFilter[0]->GetOutput() );
Antoine Regimbeau's avatar
Antoine Regimbeau committed
750

751 752 753
    SetGainLutFilterParameter( m_GainLutFilter[0] ,
                               min[0] ,
                               max[0] );
754

Antoine Regimbeau's avatar
Antoine Regimbeau committed
755
    for ( int channel = 0 ; channel < 3 ; channel++ ) 
Antoine Regimbeau's avatar
Antoine Regimbeau committed
756
      {
757

Antoine Regimbeau's avatar
Antoine Regimbeau committed
758
      m_BufferFilter[channel] = BufferFilterType::New();
759
      m_ApplyFilter[channel] = ApplyFilterType::New();
760 761 762
      SetApplyFilterParameter( m_ApplyFilter[channel] ,
                               min[0] ,
                               max[0] );
763 764 765 766 767 768 769 770 771 772 773

      m_BufferFilter[channel]->SetInput( 
          inputImageList->GetNthElement( rgb[channel] ) );
      m_ApplyFilter[channel]->SetInputImage ( 
          m_BufferFilter[channel]->GetOutput() );
      m_ApplyFilter[channel]->SetInputLut( 
          m_StreamingFilter[0]->GetOutput() );
      

      outputImageList->PushBack( m_ApplyFilter[channel]->GetOutput() );
    }
774
  }
775

776
  // Function that compute histograms with HistoPersistentFilterType
777 778 779 780 781 782 783 784
  void PersistentComputation( const FloatVectorImageType::Pointer inImage ,
                              const unsigned int nbChannel ,
                              const FloatVectorImageType::PixelType & max ,
                              const FloatVectorImageType::PixelType & min)
  { 

    HistoPersistentFilterType::Pointer histoPersistent (
        HistoPersistentFilterType::New());
785
    unsigned int nbBin ( GetParameterInt( "bins" ) );
786
    histoPersistent->SetInput( inImage );
787
    FloatVectorImageType::PixelType pixel( nbChannel ) , step( nbChannel );
788 789 790 791 792 793 794 795 796 797 798 799 800
    pixel.Fill( nbBin );
    histoPersistent->GetFilter()->SetNumberOfBins( pixel );
    if ( IsParameterEnabled("nodata") )
      {
      histoPersistent->GetFilter()->NoDataFlagOn();
      histoPersistent->GetFilter()->SetNoDataValue( 
                GetParameterFloat( "nodata" ) );
      }
    step = ( max - min ) / nbBin;
    pixel = min - 0.5 * step;
    histoPersistent->GetFilter()->SetHistogramMin(pixel);
    pixel = max + 0.5 * step;
    histoPersistent->GetFilter()->SetHistogramMax(pixel);
801
    AddProcess(histoPersistent->GetStreamer(), "Computing histogram");
802
    histoPersistent->Update();
803
    HistoPersistentFilterType::HistogramListType * histoList = 
804
            histoPersistent->GetHistogramList();
805

806
    Transfer( histoList );
807

808 809
    if ( HasValue("hfact") )
      {
810 811 812 813
      Threshold( histoList , nbBin );
      }
  }

Guillaume Pasero's avatar
Guillaume Pasero committed
814
  // Threshold function that is normally done in ComputeHistoFilter  and here is
815
  // used on the output of HistoPersistentFilterType.
816 817 818
  void Threshold( HistoPersistentFilterType::HistogramListType * histoList ,
                  unsigned int nbBin )
  {
819
    for ( unsigned int j = 0 ; j < histoList->Size() ; j++ )
820 821
      {
      unsigned int rest(0) , height ( static_cast<unsigned int>( 
822
          GetParameterFloat( "hfact" ) * 
823
          ( histoList->GetNthElement(j)->GetTotalFrequency() / nbBin ) ) );
824

825 826
      HistogramType::IndexType zero;
      HistogramType::Pointer & histoToThresh = m_Histogram[j];
827
      zero.Fill(0);
828
      for ( unsigned int i = 0 ; i < nbBin ; i++ )
829 830 831 832 833 834 835 836 837
        {
        if ( histoToThresh->GetPixel(zero)[i] > height )
          {
          rest += histoToThresh->GetPixel(zero)[i] - height ;
          histoToThresh->GetPixel(zero)[i] = height ;
          }
        }
      height = rest / nbBin;
      rest = rest % nbBin;
838
      for ( unsigned int i = 0 ; i < nbBin ; i++ )
839 840 841 842 843 844 845 846 847 848
        {
        histoToThresh->GetPixel(zero)[i] += height ;
        if ( i > (nbBin - rest)/2 && i <= (nbBin - rest)/2 + rest )
          {
          ++histoToThresh->GetPixel(zero)[i];
          }
        }
      }
  }

849 850 851
  // Transfer the output of the HistoPersistentFilterType to the input type
  // of ComputeGainLutFilter
  void Transfer( HistoPersistentFilterType::HistogramListType * histoList )
852
  {
853
    unsigned int nbBin( GetParameterInt( "bins" ) );
854 855

    HistoPersistentFilterType::HistogramType::Pointer histo;
856
    FloatImageType::SpacingType inputSpacing ( 
857
      GetParameterImage("in")->GetSignedSpacing() );
858 859
    FloatImageType::PointType inputOrigin ( 
      GetParameterImage("in")->GetOrigin() );
860

861
    HistogramType::SpacingType histoSpacing ;
862 863
    histoSpacing[0] = inputSpacing[0] * m_ThumbSize[0] ;
    histoSpacing[1] = inputSpacing[1] * m_ThumbSize[1] ;
864
    HistogramType::PointType histoOrigin ;
865 866 867 868
    histoOrigin[0] = histoSpacing[0] / 2 + 
      inputOrigin[0] - inputSpacing[0] / 2 ;
    histoOrigin[1] = histoSpacing[1] / 2 + 
      inputOrigin[1] - inputSpacing[1] / 2 ;
869 870 871

    for ( unsigned int i = 0 ; i < histoList->Size() ; i++ )
      {
872 873 874 875 876 877
      HistogramType::Pointer histoVectorImage =
        CreateHistogramFrom( histoList->GetNthElement( i ) ,
                             nbBin ,
                             histoSpacing ,
                             histoOrigin );
     
878 879 880 881
      m_Histogram.push_back( histoVectorImage );
      }
  }

882
  // Creating a histogram (VectorImage) from an itkHistogram
883 884 885 886 887
  HistogramType::Pointer CreateHistogramFrom( 
      HistoPersistentFilterType::HistogramType::Pointer histo ,
      unsigned int nbBin ,
      HistogramType::SpacingType histoSpacing ,
      HistogramType::PointType histoOrigin )
888
  {
889 890 891 892 893 894 895 896 897 898 899 900 901 902
    HistogramType::SizeType sizeOne;
    sizeOne.Fill( 1 );
    HistogramType::IndexType index;
    index.Fill(0);

    HistogramType::PixelType histoPixel( nbBin );

    HistogramType::Pointer histoVectorImage ( 
            HistogramType::New() );

    histoVectorImage->SetVectorLength( nbBin );
    histoVectorImage->SetBufferedRegion( sizeOne );
    histoVectorImage->SetRequestedRegion( sizeOne );
    histoVectorImage->SetLargestPossibleRegion( sizeOne );
903
    histoVectorImage->SetSignedSpacing( histoSpacing ) ;
904 905 906
    histoVectorImage->SetOrigin( histoOrigin );
    histoVectorImage->Allocate();
    for (unsigned int j = 0 ; j < nbBin ; j++ )
907
      {
908
      histoPixel[j] = histo->GetFrequency( j );
909
      }
910 911
    histoVectorImage->SetPixel( index , histoPixel );
    return histoVectorImage;
912 913
  }

914
  // Set correct parameters for the ComputeHistoFilter
915 916 917
  void SetHistoFilterParameter( HistoFilterType::Pointer histoFilter ,
                                float min ,
                                float max ,
918 919
                                unsigned int nbBin ,
                                float thresh = -1 )
920 921 922 923 924 925
  {
    histoFilter->SetMin( min );
    histoFilter->SetMax( max );
    histoFilter->SetNbBin( nbBin );
    histoFilter->SetThumbSize( m_ThumbSize );
    histoFilter->SetThreshold( thresh );
926 927 928 929 930
    if ( IsParameterEnabled("nodata") )
      {
      histoFilter->SetNoData( GetParameterFloat( "nodata" ) );
      histoFilter->SetNoDataFlag( true );
      }
931 932
  }

933
  // Set correct parameters for the ComputeGainLutFilter
934 935 936 937 938 939 940 941 942
  void SetGainLutFilterParameter( GainLutFilterType::Pointer gainLutFilter ,
                                  ImagePixelType min ,
                                  ImagePixelType max )
  {
    gainLutFilter->SetMin( min );
    gainLutFilter->SetMax( max );
    gainLutFilter->SetNbPixel( m_ThumbSize[0]*m_ThumbSize[1] );
  }

943
  // Set correct parameters for the ApplyGainFilter
944 945
  void SetApplyFilterParameter( ApplyFilterType::Pointer applyFilter ,
                                ImagePixelType min ,
946
                                ImagePixelType max )
947 948 949 950
  {
    applyFilter->SetMin( min );
    applyFilter->SetMax( max );
    applyFilter->SetThumbSize( m_ThumbSize );
Antoine Regimbeau's avatar