samples/cpp/OpenEXRimages_HDR_Retina_toneMapping.cpp

   1
   2 //============================================================================
   3 // Name        : HighDynamicRange_RetinaCompression.cpp
   4 // Author      : Alexandre Benoit (benoit.alexandre.vision@gmail.com)
   5 // Version     : 0.1
   6 // Copyright   : Alexandre Benoit, LISTIC Lab, july 2011
   7 // Description : HighDynamicRange compression (tone mapping) with the help of the Gipsa/Listic's retina in C++, Ansi-style
   8 //============================================================================
   9
  10 #include <iostream>
  11 #include <cstring>
  12
  13 #include "opencv2/opencv.hpp"
  14
  15 static void help(std::string errorMessage)
  16 {
  17     std::cout<<"Program init error : "<<errorMessage<<std::endl;
  18     std::cout<<"\nProgram call procedure : ./OpenEXRimages_HDR_Retina_toneMapping [OpenEXR image to process]"<<std::endl;
  19     std::cout<<"\t[OpenEXR image to process] : the input HDR image to process, must be an OpenEXR format, see http://www.openexr.com/ to get some samples or create your own using camera bracketing and Photoshop or equivalent software for OpenEXR image synthesis"<<std::endl;
  20     std::cout<<"\nExamples:"<<std::endl;
  21     std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping memorial.exr"<<std::endl;
  22 }
  23
  24 // simple procedure for 1D curve tracing
  25 static void drawPlot(const cv::Mat curve, const std::string figureTitle, const int lowerLimit, const int upperLimit)
  26 {
  27     //std::cout<<"curve size(h,w) = "<<curve.size().height<<", "<<curve.size().width<<std::endl;
  28     cv::Mat displayedCurveImage = cv::Mat::ones(200, curve.size().height, CV_8U);
  29
  30     cv::Mat windowNormalizedCurve;
  31     normalize(curve, windowNormalizedCurve, 0, 200, CV_MINMAX, CV_32F);
  32
  33     displayedCurveImage = cv::Scalar::all(255); // set a white background
  34     int binW = cvRound((double)displayedCurveImage.cols/curve.size().height);
  35
  36     for( int i = 0; i < curve.size().height; i++ )
  37         rectangle( displayedCurveImage, cv::Point(i*binW, displayedCurveImage.rows),
  38                 cv::Point((i+1)*binW, displayedCurveImage.rows - cvRound(windowNormalizedCurve.at<float>(i))),
  39                 cv::Scalar::all(0), -1, 8, 0 );
  40     rectangle( displayedCurveImage, cv::Point(0, 0),
  41             cv::Point((lowerLimit)*binW, 200),
  42             cv::Scalar::all(128), -1, 8, 0 );
  43     rectangle( displayedCurveImage, cv::Point(displayedCurveImage.cols, 0),
  44             cv::Point((upperLimit)*binW, 200),
  45             cv::Scalar::all(128), -1, 8, 0 );
  46
  47     cv::imshow(figureTitle, displayedCurveImage);
  48 }
  49 /*
  50  * objective : get the gray level map of the input image and rescale it to the range [0-255]
  51  */
  52  static void rescaleGrayLevelMat(const cv::Mat &inputMat, cv::Mat &outputMat, const float histogramClippingLimit)
  53  {
  54
  55      // adjust output matrix wrt the input size but single channel
  56      std::cout<<"Input image rescaling with histogram edges cutting (in order to eliminate bad pixels created during the HDR image creation) :"<<std::endl;
  57      //std::cout<<"=> image size (h,w,channels) = "<<inputMat.size().height<<", "<<inputMat.size().width<<", "<<inputMat.channels()<<std::endl;
  58      //std::cout<<"=> pixel coding (nbchannel, bytes per channel) = "<<inputMat.elemSize()/inputMat.elemSize1()<<", "<<inputMat.elemSize1()<<std::endl;
  59
  60      // rescale between 0-255, keeping floating point values
  61      cv::normalize(inputMat, outputMat, 0.0, 255.0, cv::NORM_MINMAX);
  62
  63      // extract a 8bit image that will be used for histogram edge cut
  64      cv::Mat intGrayImage;
  65      if (inputMat.channels()==1)
  66      {
  67          outputMat.convertTo(intGrayImage, CV_8U);
  68      }else
  69      {
  70          cv::Mat rgbIntImg;
  71          outputMat.convertTo(rgbIntImg, CV_8UC3);
  72          cvtColor(rgbIntImg, intGrayImage, CV_BGR2GRAY);
  73      }
  74
  75      // get histogram density probability in order to cut values under above edges limits (here 5-95%)... usefull for HDR pixel errors cancellation
  76      cv::Mat dst, hist;
  77      int histSize = 256;
  78      calcHist(&intGrayImage, 1, 0, cv::Mat(), hist, 1, &histSize, 0);
  79      cv::Mat normalizedHist;
  80      normalize(hist, normalizedHist, 1, 0, cv::NORM_L1, CV_32F); // normalize histogram so that its sum equals 1
  81
  82      double min_val, max_val;
  83      CvMat histArr(normalizedHist);
  84      cvMinMaxLoc(&histArr, &min_val, &max_val);
  85      //std::cout<<"Hist max,min = "<<max_val<<", "<<min_val<<std::endl;
  86
  87      // compute density probability
  88      cv::Mat denseProb=cv::Mat::zeros(normalizedHist.size(), CV_32F);
  89      denseProb.at<float>(0)=normalizedHist.at<float>(0);
  90      int histLowerLimit=0, histUpperLimit=0;
  91      for (int i=1;i<normalizedHist.size().height;++i)
  92      {
  93          denseProb.at<float>(i)=denseProb.at<float>(i-1)+normalizedHist.at<float>(i);
  94          //std::cout<<normalizedHist.at<float>(i)<<", "<<denseProb.at<float>(i)<<std::endl;
  95          if ( denseProb.at<float>(i)<histogramClippingLimit)
  96              histLowerLimit=i;
  97          if ( denseProb.at<float>(i)<1-histogramClippingLimit)
  98              histUpperLimit=i;
  99      }
 100      // deduce min and max admitted gray levels
 101      float minInputValue = (float)histLowerLimit/histSize*255;
 102      float maxInputValue = (float)histUpperLimit/histSize*255;
 103
 104      std::cout<<"=> Histogram limits "
 105              <<"\n\t"<<histogramClippingLimit*100<<"% index = "<<histLowerLimit<<" => normalizedHist value = "<<denseProb.at<float>(histLowerLimit)<<" => input gray level = "<<minInputValue
 106              <<"\n\t"<<(1-histogramClippingLimit)*100<<"% index = "<<histUpperLimit<<" => normalizedHist value = "<<denseProb.at<float>(histUpperLimit)<<" => input gray level = "<<maxInputValue
 107              <<std::endl;
 108      //drawPlot(denseProb, "input histogram density probability", histLowerLimit, histUpperLimit);
 109      drawPlot(normalizedHist, "input histogram", histLowerLimit, histUpperLimit);
 110
 111      // rescale image range [minInputValue-maxInputValue] to [0-255]
 112      outputMat-=minInputValue;
 113      outputMat*=255.0/(maxInputValue-minInputValue);
 114      // cut original histogram and back project to original image
 115      cv::threshold( outputMat, outputMat, 255.0, 255.0, 2 ); //THRESH_TRUNC, clips values above 255
 116      cv::threshold( outputMat, outputMat, 0.0, 0.0, 3 ); //THRESH_TOZERO, clips values under 0
 117
 118  }
 119  // basic callback method for interface management
 120  cv::Mat inputImage;
 121  cv::Mat imageInputRescaled;
 122  int histogramClippingValue;
 123  static void callBack_rescaleGrayLevelMat(int, void*)
 124  {
 125      std::cout<<"Histogram clipping value changed, current value = "<<histogramClippingValue<<std::endl;
 126      rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)(histogramClippingValue/100.0));
 127      normalize(imageInputRescaled, imageInputRescaled, 0.0, 255.0, cv::NORM_MINMAX);
 128  }
 129
 130  cv::Ptr<cv::Retina> retina;
 131  int retinaHcellsGain;
 132  int localAdaptation_photoreceptors, localAdaptation_Gcells;
 133  static void callBack_updateRetinaParams(int, void*)
 134  {
 135      retina->setupOPLandIPLParvoChannel(true, true, (float)(localAdaptation_photoreceptors/200.0), 0.5f, 0.43f, (float)retinaHcellsGain, 1.f, 7.f, (float)(localAdaptation_Gcells/200.0));
 136  }
 137
 138  int colorSaturationFactor;
 139  static void callback_saturateColors(int, void*)
 140  {
 141      retina->setColorSaturation(true, (float)colorSaturationFactor);
 142  }
 143
 144  int main(int argc, char* argv[]) {
 145      // welcome message
 146      std::cout<<"*********************************************************************************"<<std::endl;
 147      std::cout<<"* Retina demonstration for High Dynamic Range compression (tone-mapping) : demonstrates the use of a wrapper class of the Gipsa/Listic Labs retina model."<<std::endl;
 148      std::cout<<"* This retina model allows spatio-temporal image processing (applied on still images, video sequences)."<<std::endl;
 149      std::cout<<"* This demo focuses demonstration of the dynamic compression capabilities of the model"<<std::endl;
 150      std::cout<<"* => the main application is tone mapping of HDR images (i.e. see on a 8bit display a more than 8bits coded (up to 16bits) image with details in high and low luminance ranges"<<std::endl;
 151      std::cout<<"* The retina model still have the following properties:"<<std::endl;
 152      std::cout<<"* => It applies a spectral whithening (mid-frequency details enhancement)"<<std::endl;
 153      std::cout<<"* => high frequency spatio-temporal noise reduction"<<std::endl;
 154      std::cout<<"* => low frequency luminance to be reduced (luminance range compression)"<<std::endl;
 155      std::cout<<"* => local logarithmic luminance compression allows details to be enhanced in low light conditions\n"<<std::endl;
 156      std::cout<<"* for more information, reer to the following papers :"<<std::endl;
 157      std::cout<<"* Benoit A., Caplier A., Durette B., Herault, J., \"USING HUMAN VISUAL SYSTEM MODELING FOR BIO-INSPIRED LOW LEVEL IMAGE PROCESSING\", Elsevier, Computer Vision and Image Understanding 114 (2010), pp. 758-773, DOI: http://dx.doi.org/10.1016/j.cviu.2010.01.011"<<std::endl;
 158      std::cout<<"* Vision: Images, Signals and Neural Networks: Models of Neural Processing in Visual Perception (Progress in Neural Processing),By: Jeanny Herault, ISBN: 9814273686. WAPI (Tower ID): 113266891."<<std::endl;
 159      std::cout<<"* => reports comments/remarks at benoit.alexandre.vision@gmail.com"<<std::endl;
 160      std::cout<<"* => more informations and papers at : http://sites.google.com/site/benoitalexandrevision/"<<std::endl;
 161      std::cout<<"*********************************************************************************"<<std::endl;
 162      std::cout<<"** WARNING : this sample requires OpenCV to be configured with OpenEXR support **"<<std::endl;
 163      std::cout<<"*********************************************************************************"<<std::endl;
 164      std::cout<<"*** You can use free tools to generate OpenEXR images from images sets   :    ***"<<std::endl;
 165      std::cout<<"*** =>  1. take a set of photos from the same viewpoint using bracketing      ***"<<std::endl;
 166      std::cout<<"*** =>  2. generate an OpenEXR image with tools like qtpfsgui.sourceforge.net ***"<<std::endl;
 167      std::cout<<"*** =>  3. apply tone mapping with this program                               ***"<<std::endl;
 168      std::cout<<"*********************************************************************************"<<std::endl;
 169
 170      // basic input arguments checking
 171      if (argc<2)
 172      {
 173          help("bad number of parameter");
 174          return -1;
 175      }
 176
 177      bool useLogSampling = !strcmp(argv[argc-1], "log"); // check if user wants retina log sampling processing
 178
 179      std::string inputImageName=argv[1];
 180
 181      //////////////////////////////////////////////////////////////////////////////
 182      // checking input media type (still image, video file, live video acquisition)
 183      std::cout<<"RetinaDemo: processing image "<<inputImageName<<std::endl;
 184      // image processing case
 185      // declare the retina input buffer... that will be fed differently in regard of the input media
 186      inputImage = cv::imread(inputImageName, -1); // load image in RGB mode
 187      std::cout<<"=> image size (h,w) = "<<inputImage.size().height<<", "<<inputImage.size().width<<std::endl;
 188      if (!inputImage.total())
 189      {
 190         help("could not load image, program end");
 191             return -1;
 192          }
 193      // rescale between 0 and 1
 194      normalize(inputImage, inputImage, 0.0, 1.0, cv::NORM_MINMAX);
 195      cv::Mat gammaTransformedImage;
 196      cv::pow(inputImage, 1./5, gammaTransformedImage); // apply gamma curve: img = img ** (1./5)
 197      imshow("EXR image original image, 16bits=>8bits linear rescaling ", inputImage);
 198      imshow("EXR image with basic processing : 16bits=>8bits with gamma correction", gammaTransformedImage);
 199      if (inputImage.empty())
 200      {
 201          help("Input image could not be loaded, aborting");
 202          return -1;
 203      }
 204
 205      //////////////////////////////////////////////////////////////////////////////
 206      // Program start in a try/catch safety context (Retina may throw errors)
 207      try
 208      {
 209          /* create a retina instance with default parameters setup, uncomment the initialisation you wanna test
 210           * -> if the last parameter is 'log', then activate log sampling (favour foveal vision and subsamples peripheral vision)
 211           */
 212          if (useLogSampling)
 213                 {
 214                      retina = new cv::Retina(inputImage.size(),true, cv::RETINA_COLOR_BAYER, true, 2.0, 10.0);
 215                  }
 216          else// -> else allocate "classical" retina :
 217              retina = new cv::Retina(inputImage.size());
 218
 219         // save default retina parameters file in order to let you see this and maybe modify it and reload using method "setup"
 220         retina->write("RetinaDefaultParameters.xml");
 221
 222                  // desactivate Magnocellular pathway processing (motion information extraction) since it is not usefull here
 223                  retina->activateMovingContoursProcessing(false);
 224
 225          // declare retina output buffers
 226          cv::Mat retinaOutput_parvo;
 227
 228          /////////////////////////////////////////////
 229          // prepare displays and interactions
 230          histogramClippingValue=0; // default value... updated with interface slider
 231          //inputRescaleMat = inputImage;
 232          //outputRescaleMat = imageInputRescaled;
 233          cv::namedWindow("Retina input image (with cut edges histogram for basic pixels error avoidance)",1);
 234          cv::createTrackbar("histogram edges clipping limit", "Retina input image (with cut edges histogram for basic pixels error avoidance)",&histogramClippingValue,50,callBack_rescaleGrayLevelMat);
 235
 236          cv::namedWindow("Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", 1);
 237          colorSaturationFactor=3;
 238          cv::createTrackbar("Color saturation", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &colorSaturationFactor,5,callback_saturateColors);
 239
 240          retinaHcellsGain=40;
 241          cv::createTrackbar("Hcells gain", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping",&retinaHcellsGain,100,callBack_updateRetinaParams);
 242
 243          localAdaptation_photoreceptors=197;
 244          localAdaptation_Gcells=190;
 245          cv::createTrackbar("Ph sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_photoreceptors,199,callBack_updateRetinaParams);
 246          cv::createTrackbar("Gcells sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_Gcells,199,callBack_updateRetinaParams);
 247
 248
 249          /////////////////////////////////////////////
 250          // apply default parameters of user interaction variables
 251          rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)histogramClippingValue/100);
 252          retina->setColorSaturation(true,(float)colorSaturationFactor);
 253          callBack_updateRetinaParams(1,NULL); // first call for default parameters setup
 254
 255          // processing loop with stop condition
 256          bool continueProcessing=true;
 257          while(continueProcessing)
 258          {
 259              // run retina filter
 260              retina->run(imageInputRescaled);
 261              // Retrieve and display retina output
 262              retina->getParvo(retinaOutput_parvo);
 263              cv::imshow("Retina input image (with cut edges histogram for basic pixels error avoidance)", imageInputRescaled/255.0);
 264              cv::imshow("Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", retinaOutput_parvo);
 265              cv::waitKey(10);
 266          }
 267      }catch(cv::Exception e)
 268      {
 269          std::cerr<<"Error using Retina : "<<e.what()<<std::endl;
 270      }
 271
 272      // Program end message
 273      std::cout<<"Retina demo end"<<std::endl;
 274
 275      return 0;
 276  }