First draw of the histogram equalisation working on black&white pictures and color ones

Modules/Feature/Contrast/CMakeLists.txt

+cmake_minimum_required (VERSION 2.6)
+project (projchambolle)
+
+find_package(OTB REQUIRED)
+include(${OTB_USE_FILE})
+
+
+set(SOURCES otbContrastEnhancementFilter.cxx)
+
+add_executable(Contrast otbContrastEnhancementFilter.cxx)
+
+target_link_libraries(Contrast ${OTB_LIBRARIES})
+
+
+# message(WARNING "Include dirs : ${OTB_INCLUDE_DIRS}")
\ No newline at end of file

Modules/Feature/Contrast/otbContrastEnhancementFilter.cxx

+#include <iostream>
+#include <math.h>
+#include "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkAdditiveGaussianNoiseImageFilter.h"
+#include "otbStreamingStatisticsImageFilter.h"
+#include "itkImageLinearIteratorWithIndex.h"
+#include "itkMultiplyImageFilter.h"
+#include "itkDivideImageFilter.h"
+#include "itkAbsImageFilter.h"
+#include "itkSubtractImageFilter.h"
+#include "itkAddImageFilter.h"
+#include "itkSquareImageFilter.h"
+#include "itkSqrtImageFilter.h"
+#include "itkMaximumImageFilter.h"
+
+typedef float                     PixelType;
+typedef itk::Image< PixelType , 2>  ImageType;
+typedef itk::AddImageFilter<ImageType, ImageType, ImageType>  AddFilterType;
+typedef  otb::StreamingStatisticsImageFilter<ImageType>      StatFilter;
+
+ImageType::Pointer 
+createImagefrom( ImageType::Pointer input )
+{
+  ImageType::Pointer cloneImage( ImageType::New() );
+  cloneImage->SetRegions( input->GetRequestedRegion() ); 
+  cloneImage->Allocate();
+  cloneImage->SetOrigin( input->GetOrigin() );
+  return cloneImage ;
+
+}
+void
+copyAtoB( const ImageType::Pointer input ,
+          ImageType::Pointer output )
+{
+  itk::ImageRegionConstIterator< ImageType > itinp( input , input->GetRequestedRegion() );
+  itk::ImageRegionIterator< ImageType >      itout( output , output->GetRequestedRegion() );
+  itinp.GoToBegin();
+  itout.GoToBegin();
+  while (!itinp.IsAtEnd())
+  {
+    itout.Set( itinp.Get() );
+    ++itout;
+    ++itinp;
+  }
+}
+
+void
+normL( ImageType::Pointer inputX ,
+       ImageType::Pointer inputY ,
+       ImageType::Pointer output )
+{
+  typedef itk::SquareImageFilter< ImageType , ImageType > SquareFilterType;
+  typedef itk::SqrtImageFilter< ImageType , ImageType > SqrtFitlerType;
+  SquareFilterType::Pointer square( SquareFilterType::New() );
+  SqrtFitlerType::Pointer sqrt( SqrtFitlerType::New() );
+  AddFilterType::Pointer add( AddFilterType::New() );
+  square->SetInput( inputX );
+  square->Update();
+  copyAtoB( square->GetOutput() , output );
+  add->SetInput1( output );
+  square->SetInput( inputY );
+  square->Update();
+  add->SetInput2( square->GetOutput() );
+  sqrt->SetInput( add->GetOutput() );
+  sqrt->Update();
+  copyAtoB( sqrt->GetOutput() , output );
+}
+
+int
+main( int argc, 
+      char const *argv[] )
+{
+	const char * inputFilename  = argv[ 1 ];
+  const char * outputFilename = argv[ 2 ];
+
+
+  typedef otb::ImageFileReader< ImageType > ReaderType;
+
+ 	ReaderType::Pointer reader( ReaderType::New() );
+ 	reader->SetFileName( argv[ 1 ] );
+
+  ImageType::Pointer noisyImage( ImageType::New() );
+
+  typedef otb::ImageFileWriter<ImageType>  WriterType;
+ 	WriterType::Pointer writer( WriterType::New());
+  return 0;
+}
\ No newline at end of file

Modules/Filtering/Contrast/CMakeLists.txt

+#
+# 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.
+#
+
+project(OTBContrast)
+otb_module_impl()
\ No newline at end of file

Modules/Filtering/Contrast/otb-module.cmake

+#
+# 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.
+#
+
+set(DOCUMENTATION "work in progress")
+
+otb_module(OTBContrast
+  DEPENDS
+    OTBITK
+  	OTBCommon
+    OTBImageBase
+    
+
+  TEST_DEPENDS
+    OTBTestKernel
+    OTBImageIO
+    OTBImageBase
+  	
+
+  DESCRIPTION
+    "${DOCUMENTATION}"
+)

Modules/Filtering/Contrast/test/CMakeLists.txt

+#
+# 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.
+#
+
+otb_module_test()
+
+set(OTBContrastTests
+otbContrastEnhancementFilter.cxx
+)
+
+add_executable(otbContrastTestDriver ${OTBContrastTests})
+target_link_libraries(otbContrastTestDriver ${OTBContrast-Test_LIBRARIES})
+otb_module_target_label(otbContrastTestDriver)

Modules/Filtering/Contrast/test/otbContrastEnhancementFilter.cxx

+#include <iostream>
+#include <math.h>
+#include <vector>
+#include "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkImageLinearIteratorWithIndex.h"
+#include "itkImageToHistogramFilter.h"
+#include "itkVectorIndexSelectionCastImageFilter.h"
+#include "itkMultiplyImageFilter.h"
+
+
+typedef int                   PixelType;
+typedef itk::Image< PixelType , 2 >  ImageType;
+typedef itk::VectorImage< PixelType , 2 >  VectorImageType;
+
+void
+computeLuminance( ImageType::Pointer luminance ,
+                  VectorImageType::Pointer input )
+{
+  itk::ImageRegionIterator< VectorImageType >  itinp( input , input->GetRequestedRegion() );
+  itk::ImageRegionIterator< ImageType >  itlum( luminance , luminance->GetRequestedRegion() );
+  itinp.GoToBegin();
+  itlum.GoToBegin();
+  while ( !itinp.IsAtEnd() )
+  {
+    // 0.3R 0.6G 0.1B
+    itlum.Set( round(0.3 * itinp.Get()[0] + 0.6 * itinp.Get()[1] + 0.1  * itinp.Get()[2]) );
+    ++itinp;
+    ++itlum;
+  }
+}
+
+int
+main( int argc, 
+      char const *argv[] )
+{
+  // const char * inputfilename  = argv[ 1 ];
+  // const char * outputfilename = argv[ 2 ];
+
+
+  typedef otb::ImageFileReader< VectorImageType > ReaderType;
+  typedef itk::Statistics::ImageToHistogramFilter< ImageType > HistogramFilterType;
+  typedef HistogramFilterType::HistogramType HistogramType;
+  typedef itk::VectorIndexSelectionCastImageFilter< VectorImageType , ImageType > VectorToImageFilterType;
+
+
+  ReaderType::Pointer reader( ReaderType::New() );
+  reader->SetFileName( argv[ 1 ] );
+  reader->Update();
+  VectorImageType::Pointer input = reader->GetOutput();
+  ImageType::Pointer lum ( ImageType::New() ); 
+  lum->SetRegions( input->GetRequestedRegion() );
+  lum->Allocate();
+  lum->SetOrigin( input->GetOrigin() );
+  VectorToImageFilterType::Pointer vectorToImageFilter ( VectorToImageFilterType::New() );
+  if ( input->GetVectorLength() == 3 )
+  {
+    computeLuminance(lum, input);
+  }
+  else
+  {
+    // default is a black and white image, if length>3 and is not equal to one
+    // other case :  how to apply the equalization :
+    // - channel per channel?
+    // - with a general luminance?
+    vectorToImageFilter->SetInput( input );
+    vectorToImageFilter->SetIndex( 0 );
+    vectorToImageFilter->Update();
+    lum = vectorToImageFilter->GetOutput();
+  }
+
+
+  HistogramFilterType::Pointer histofilter( HistogramFilterType::New() );
+  histofilter->SetInput( lum );
+  histofilter->Update();
+  // Compute the histogram of the bw image
+  HistogramType * histo = histofilter->GetOutput();
+
+  // Define propertise of the target histogram : nomber of bin, height per bin
+  ImageType::SizeType lumSize = lum->GetLargestPossibleRegion().GetSize(); 
+  int npx = lumSize[0] * lumSize[1];
+  int nbin = histo->GetSize()[0];
+  float height = npx / nbin;
+
+  // The transfer function will be a look up table
+  std::vector< int > lut;
+  float nbOri, nbNew ; 
+  int countnew, countold;
+  countnew = 0;
+  countold = 0;
+  nbNew = height;
+  HistogramType::Iterator it = histo->Begin();
+  nbOri = it.GetFrequency();
+  while ( it != histo->End() )
+  {
+    if (nbOri> nbNew)
+    {
+      nbNew += height;
+      // index[0] += 1; 
+      ++countnew;
+    }
+    else
+    {
+      ++it;
+      nbOri += it.GetFrequency();
+      ++countold;
+      lut.push_back(countnew);
+    }
+  }
+
+  // Here we compute a gain image corresponding to the associated gain of the equalization
+  typedef itk::Image< float , 2 > ImageGainType;
+  ImageGainType::Pointer gainImage ( ImageGainType::New() );
+  gainImage->SetRegions( lum->GetRequestedRegion() );
+  gainImage->Allocate();
+  gainImage->SetOrigin( lum->GetOrigin() );
+  itk::ImageRegionIterator< ImageType > itlum( lum , lum->GetRequestedRegion() );
+  itk::ImageRegionIterator< ImageGainType > itgain( gainImage , gainImage->GetRequestedRegion() );
+  itlum.GoToBegin();
+  itgain.GoToBegin();
+  while( !itlum.IsAtEnd() )
+  {
+    itgain.Set( static_cast<float>(lut[itlum.Get()]) / static_cast<float>( itlum.Get() ) ) ;
+    ++itlum;
+    ++itgain;
+  }
+
+  typedef  itk::MultiplyImageFilter< VectorImageType, ImageGainType, VectorImageType > MultiplyImageFilterType;
+  MultiplyImageFilterType::Pointer gainMultiplyer ( MultiplyImageFilterType::New() );
+  gainMultiplyer->SetInput1( input );
+  gainMultiplyer->SetInput2( gainImage );
+  gainMultiplyer->Update();
+  typedef otb::ImageFileWriter<VectorImageType>  writertype;
+  writertype::Pointer writer( writertype::New());
+  writer->SetFileName( argv[2] );
+  writer->SetInput( gainMultiplyer->GetOutput() );
+  writer->Update();
+  return 0;
+}