samples/cpp/OpenEXRimages_HDR_Retina_toneMapping_video.cpp

   1
   2 //============================================================================
   3 // Name        : OpenEXRimages_HDR_Retina_toneMapping_video.cpp
   4 // Author      : Alexandre Benoit (benoit.alexandre.vision@gmail.com)
   5 // Version     : 0.2
   6 // Copyright   : Alexandre Benoit, LISTIC Lab, december 2011
   7 // Description : HighDynamicRange retina tone mapping for image sequences with the help of the Gipsa/Listic's retina in C++, Ansi-style
   8 // Known issues: the input OpenEXR sequences can have bad computed pixels that should be removed
   9 //               => a simple method consists of cutting histogram edges (a slider for this on the UI is provided)
  10 //               => however, in image sequences, this histogramm cut must be done in an elegant way from frame to frame... still not done...
  11 //============================================================================
  12
  13 #include <iostream>
  14 #include <stdio.h>
  15 #include <cstring>
  16
  17 #include "opencv2/bioinspired.hpp" // retina based algorithms
  18 #include "opencv2/imgproc.hpp" // cvCvtcolor function
  19 #include "opencv2/highgui.hpp" // display
  20
  21 #ifndef _CRT_SECURE_NO_WARNINGS
  22 # define _CRT_SECURE_NO_WARNINGS
  23 #endif
  24
  25 static void help(std::string errorMessage)
  26 {
  27     std::cout<<"Program init error : "<<errorMessage<<std::endl;
  28     std::cout<<"\nProgram call procedure : ./OpenEXRimages_HDR_Retina_toneMapping [OpenEXR image sequence to process] [OPTIONNAL start frame] [OPTIONNAL end frame]"<<std::endl;
  29     std::cout<<"\t[OpenEXR image sequence to process] : std::sprintf style ready prototype filename of the input HDR images 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;
  30     std::cout<<"\t\t => WARNING : image index number of digits cannot exceed 10"<<std::endl;
  31     std::cout<<"\t[start frame] : the starting frame tat should be considered"<<std::endl;
  32     std::cout<<"\t[end frame] : the ending frame tat should be considered"<<std::endl;
  33     std::cout<<"\nExamples:"<<std::endl;
  34     std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping_video memorial%3d.exr 20 45"<<std::endl;
  35     std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping_video memorial%3d.exr 20 45 log"<<std::endl;
  36     std::cout<<"\t ==> to process images from memorial020d.exr to memorial045d.exr"<<std::endl;
  37
  38 }
  39
  40 // simple procedure for 1D curve tracing
  41 static void drawPlot(const cv::Mat curve, const std::string figureTitle, const int lowerLimit, const int upperLimit)
  42 {
  43     //std::cout<<"curve size(h,w) = "<<curve.size().height<<", "<<curve.size().width<<std::endl;
  44     cv::Mat displayedCurveImage = cv::Mat::ones(200, curve.size().height, CV_8U);
  45
  46     cv::Mat windowNormalizedCurve;
  47     normalize(curve, windowNormalizedCurve, 0, 200, cv::NORM_MINMAX, CV_32F);
  48
  49     displayedCurveImage = cv::Scalar::all(255); // set a white background
  50     int binW = cvRound((double)displayedCurveImage.cols/curve.size().height);
  51
  52     for( int i = 0; i < curve.size().height; i++ )
  53         rectangle( displayedCurveImage, cv::Point(i*binW, displayedCurveImage.rows),
  54                 cv::Point((i+1)*binW, displayedCurveImage.rows - cvRound(windowNormalizedCurve.at<float>(i))),
  55                 cv::Scalar::all(0), -1, 8, 0 );
  56     rectangle( displayedCurveImage, cv::Point(0, 0),
  57             cv::Point((lowerLimit)*binW, 200),
  58             cv::Scalar::all(128), -1, 8, 0 );
  59     rectangle( displayedCurveImage, cv::Point(displayedCurveImage.cols, 0),
  60             cv::Point((upperLimit)*binW, 200),
  61             cv::Scalar::all(128), -1, 8, 0 );
  62
  63     cv::imshow(figureTitle, displayedCurveImage);
  64 }
  65
  66 /*
  67  * objective : get the gray level map of the input image and rescale it to the range [0-255] if rescale0_255=TRUE, simply trunks else
  68  */
  69 static void rescaleGrayLevelMat(const cv::Mat &inputMat, cv::Mat &outputMat, const float histogramClippingLimit, const bool rescale0_255)
  70  {
  71      // adjust output matrix wrt the input size but single channel
  72      std::cout<<"Input image rescaling with histogram edges cutting (in order to eliminate bad pixels created during the HDR image creation) :"<<std::endl;
  73      //std::cout<<"=> image size (h,w,channels) = "<<inputMat.size().height<<", "<<inputMat.size().width<<", "<<inputMat.channels()<<std::endl;
  74      //std::cout<<"=> pixel coding (nbchannel, bytes per channel) = "<<inputMat.elemSize()/inputMat.elemSize1()<<", "<<inputMat.elemSize1()<<std::endl;
  75
  76      // get min and max values to use afterwards if no 0-255 rescaling is used
  77      double maxInput, minInput, histNormRescalefactor=1.f;
  78      double histNormOffset=0.f;
  79      minMaxLoc(inputMat, &minInput, &maxInput);
  80      histNormRescalefactor=255.f/(maxInput-minInput);
  81      histNormOffset=minInput;
  82      std::cout<<"Hist max,min = "<<maxInput<<", "<<minInput<<" => scale, offset = "<<histNormRescalefactor<<", "<<histNormOffset<<std::endl;
  83      // rescale between 0-255, keeping floating point values
  84      cv::Mat normalisedImage;
  85      cv::normalize(inputMat, normalisedImage, 0.f, 255.f, cv::NORM_MINMAX);
  86      if (rescale0_255)
  87         normalisedImage.copyTo(outputMat);
  88      // extract a 8bit image that will be used for histogram edge cut
  89      cv::Mat intGrayImage;
  90      if (inputMat.channels()==1)
  91      {
  92          normalisedImage.convertTo(intGrayImage, CV_8U);
  93      }else
  94      {
  95          cv::Mat rgbIntImg;
  96          normalisedImage.convertTo(rgbIntImg, CV_8UC3);
  97          cvtColor(rgbIntImg, intGrayImage, cv::COLOR_BGR2GRAY);
  98      }
  99
 100      // get histogram density probability in order to cut values under above edges limits (here 5-95%)... usefull for HDR pixel errors cancellation
 101      cv::Mat dst, hist;
 102      int histSize = 256;
 103      calcHist(&intGrayImage, 1, 0, cv::Mat(), hist, 1, &histSize, 0);
 104      cv::Mat normalizedHist;
 105
 106      normalize(hist, normalizedHist, 1.f, 0.f, cv::NORM_L1, CV_32F); // normalize histogram so that its sum equals 1
 107
 108      // compute density probability
 109      cv::Mat denseProb=cv::Mat::zeros(normalizedHist.size(), CV_32F);
 110      denseProb.at<float>(0)=normalizedHist.at<float>(0);
 111      int histLowerLimit=0, histUpperLimit=0;
 112      for (int i=1;i<normalizedHist.size().height;++i)
 113      {
 114          denseProb.at<float>(i)=denseProb.at<float>(i-1)+normalizedHist.at<float>(i);
 115          //std::cout<<normalizedHist.at<float>(i)<<", "<<denseProb.at<float>(i)<<std::endl;
 116          if ( denseProb.at<float>(i)<histogramClippingLimit)
 117              histLowerLimit=i;
 118          if ( denseProb.at<float>(i)<1.f-histogramClippingLimit)
 119              histUpperLimit=i;
 120      }
 121      // deduce min and max admitted gray levels
 122      float minInputValue = (float)histLowerLimit/histSize*255.f;
 123      float maxInputValue = (float)histUpperLimit/histSize*255.f;
 124
 125      std::cout<<"=> Histogram limits "
 126              <<"\n\t"<<histogramClippingLimit*100.f<<"% index = "<<histLowerLimit<<" => normalizedHist value = "<<denseProb.at<float>(histLowerLimit)<<" => input gray level = "<<minInputValue
 127              <<"\n\t"<<(1.f-histogramClippingLimit)*100.f<<"% index = "<<histUpperLimit<<" => normalizedHist value = "<<denseProb.at<float>(histUpperLimit)<<" => input gray level = "<<maxInputValue
 128              <<std::endl;
 129      //drawPlot(denseProb, "input histogram density probability", histLowerLimit, histUpperLimit);
 130      drawPlot(normalizedHist, "input histogram", histLowerLimit, histUpperLimit);
 131
 132     if(rescale0_255) // rescale between 0-255 if asked to
 133     {
 134         cv::threshold( outputMat, outputMat, maxInputValue, maxInputValue, 2 ); //THRESH_TRUNC, clips values above maxInputValue
 135         cv::threshold( outputMat, outputMat, minInputValue, minInputValue, 3 ); //THRESH_TOZERO, clips values under minInputValue
 136         // rescale image range [minInputValue-maxInputValue] to [0-255]
 137         outputMat-=minInputValue;
 138         outputMat*=255.f/(maxInputValue-minInputValue);
 139     }else
 140     {
 141         inputMat.copyTo(outputMat);
 142         // update threshold in the initial input image range
 143         maxInputValue=(float)((maxInputValue-255.f)/histNormRescalefactor+maxInput);
 144         minInputValue=(float)(minInputValue/histNormRescalefactor+minInput);
 145         std::cout<<"===> Input Hist clipping values (max,min) = "<<maxInputValue<<", "<<minInputValue<<std::endl;
 146         cv::threshold( outputMat, outputMat, maxInputValue, maxInputValue, 2 ); //THRESH_TRUNC, clips values above maxInputValue
 147         cv::threshold( outputMat, outputMat, minInputValue, minInputValue, 3 ); //
 148     }
 149  }
 150
 151  // basic callback method for interface management
 152  cv::Mat inputImage;
 153  cv::Mat imageInputRescaled;
 154  float globalRescalefactor=1;
 155  cv::Scalar globalOffset=0;
 156  int histogramClippingValue;
 157  static void callBack_rescaleGrayLevelMat(int, void*)
 158  {
 159      std::cout<<"Histogram clipping value changed, current value = "<<histogramClippingValue<<std::endl;
 160     // rescale and process
 161     inputImage+=globalOffset;
 162     inputImage*=globalRescalefactor;
 163     inputImage+=cv::Scalar(50, 50, 50, 50); // WARNING value linked to the hardcoded value (200.0) used in the globalRescalefactor in order to center on the 128 mean value... experimental but... basic compromise
 164     rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)histogramClippingValue/100.f, true);
 165
 166  }
 167
 168  cv::Ptr<cv::bioinspired::Retina> retina;
 169  int retinaHcellsGain;
 170  int localAdaptation_photoreceptors, localAdaptation_Gcells;
 171  static void callBack_updateRetinaParams(int, void*)
 172  {
 173      retina->setupOPLandIPLParvoChannel(true, true, (float)(localAdaptation_photoreceptors/200.0), 0.5f, 0.43f, (float)retinaHcellsGain, 1.f, 7.f, (float)(localAdaptation_Gcells/200.0));
 174  }
 175
 176  int colorSaturationFactor;
 177  static void callback_saturateColors(int, void*)
 178  {
 179      retina->setColorSaturation(true, (float)colorSaturationFactor);
 180  }
 181
 182 // loadNewFrame : loads a n image wrt filename parameters. it also manages image rescaling/histogram edges cutting (acts differently at first image i.e. if firstTimeread=true)
 183 static void loadNewFrame(const std::string filenamePrototype, const int currentFileIndex, const bool firstTimeread)
 184 {
 185      char *currentImageName=NULL;
 186     currentImageName = (char*)malloc(sizeof(char)*filenamePrototype.size()+10);
 187
 188     // grab the first frame
 189     sprintf(currentImageName, filenamePrototype.c_str(), currentFileIndex);
 190
 191      //////////////////////////////////////////////////////////////////////////////
 192      // checking input media type (still image, video file, live video acquisition)
 193      std::cout<<"RetinaDemo: reading image : "<<currentImageName<<std::endl;
 194      // image processing case
 195      // declare the retina input buffer... that will be fed differently in regard of the input media
 196      inputImage = cv::imread(currentImageName, -1); // load image in RGB mode
 197      std::cout<<"=> image size (h,w) = "<<inputImage.size().height<<", "<<inputImage.size().width<<std::endl;
 198      if (inputImage.empty())
 199      {
 200         help("could not load image, program end");
 201             return;;
 202          }
 203
 204     // rescaling/histogram clipping stage
 205     // rescale between 0 and 1
 206     // TODO : take care of this step !!! maybe disable of do this in a nicer way ... each successive image should get the same transformation... but it depends on the initial image format
 207     double maxInput, minInput;
 208     minMaxLoc(inputImage, &minInput, &maxInput);
 209     std::cout<<"ORIGINAL IMAGE pixels values range (max,min) : "<<maxInput<<", "<<minInput<<std::endl;
 210
 211     if (firstTimeread)
 212     {
 213         /* the first time, get the pixel values range and rougthly update scaling value
 214         in order to center values around 128 and getting a range close to [0-255],
 215         => actually using a little less in order to let some more flexibility in range evolves...
 216         */
 217         double maxInput1, minInput1;
 218         minMaxLoc(inputImage, &minInput1, &maxInput1);
 219         std::cout<<"FIRST IMAGE pixels values range (max,min) : "<<maxInput1<<", "<<minInput1<<std::endl;
 220         globalRescalefactor=(float)(50.0/(maxInput1-minInput1)); // less than 255 for flexibility... experimental value to be carefull about
 221         double channelOffset = -1.5*minInput;
 222         globalOffset= cv::Scalar(channelOffset, channelOffset, channelOffset, channelOffset);
 223     }
 224     // call the generic input image rescaling callback
 225     callBack_rescaleGrayLevelMat(1,NULL);
 226 }
 227
 228  int main(int argc, char* argv[]) {
 229      // welcome message
 230      std::cout<<"*********************************************************************************"<<std::endl;
 231      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;
 232      std::cout<<"* This retina model allows spatio-temporal image processing (applied on still images, video sequences)."<<std::endl;
 233      std::cout<<"* This demo focuses demonstration of the dynamic compression capabilities of the model"<<std::endl;
 234      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;
 235      std::cout<<"* The retina model still have the following properties:"<<std::endl;
 236      std::cout<<"* => It applies a spectral whithening (mid-frequency details enhancement)"<<std::endl;
 237      std::cout<<"* => high frequency spatio-temporal noise reduction"<<std::endl;
 238      std::cout<<"* => low frequency luminance to be reduced (luminance range compression)"<<std::endl;
 239      std::cout<<"* => local logarithmic luminance compression allows details to be enhanced in low light conditions\n"<<std::endl;
 240      std::cout<<"* for more information, reer to the following papers :"<<std::endl;
 241      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;
 242      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;
 243      std::cout<<"* => reports comments/remarks at benoit.alexandre.vision@gmail.com"<<std::endl;
 244      std::cout<<"* => more informations and papers at : http://sites.google.com/site/benoitalexandrevision/"<<std::endl;
 245      std::cout<<"*********************************************************************************"<<std::endl;
 246      std::cout<<"** WARNING : this sample requires OpenCV to be configured with OpenEXR support **"<<std::endl;
 247      std::cout<<"*********************************************************************************"<<std::endl;
 248      std::cout<<"*** You can use free tools to generate OpenEXR images from images sets   :    ***"<<std::endl;
 249      std::cout<<"*** =>  1. take a set of photos from the same viewpoint using bracketing      ***"<<std::endl;
 250      std::cout<<"*** =>  2. generate an OpenEXR image with tools like qtpfsgui.sourceforge.net ***"<<std::endl;
 251      std::cout<<"*** =>  3. apply tone mapping with this program                               ***"<<std::endl;
 252      std::cout<<"*********************************************************************************"<<std::endl;
 253
 254      // basic input arguments checking
 255      if (argc<4)
 256      {
 257          help("bad number of parameter");
 258          return -1;
 259      }
 260
 261      bool useLogSampling = !strcmp(argv[argc-1], "log"); // check if user wants retina log sampling processing
 262
 263      int startFrameIndex=0, endFrameIndex=0, currentFrameIndex=0;
 264      sscanf(argv[2], "%d", &startFrameIndex);
 265      sscanf(argv[3], "%d", &endFrameIndex);
 266      std::string inputImageNamePrototype(argv[1]);
 267
 268      //////////////////////////////////////////////////////////////////////////////
 269      // checking input media type (still image, video file, live video acquisition)
 270      std::cout<<"RetinaDemo: setting up system with first image..."<<std::endl;
 271      loadNewFrame(inputImageNamePrototype, startFrameIndex, true);
 272
 273      if (inputImage.empty())
 274      {
 275         help("could not load image, program end");
 276             return -1;
 277          }
 278
 279      //////////////////////////////////////////////////////////////////////////////
 280      // Program start in a try/catch safety context (Retina may throw errors)
 281      try
 282      {
 283          /* create a retina instance with default parameters setup, uncomment the initialisation you wanna test
 284           * -> if the last parameter is 'log', then activate log sampling (favour foveal vision and subsamples peripheral vision)
 285           */
 286          if (useLogSampling)
 287                 {
 288                      retina = cv::bioinspired::createRetina(inputImage.size(),true, cv::bioinspired::RETINA_COLOR_BAYER, true, 2.0, 10.0);
 289                  }
 290          else// -> else allocate "classical" retina :
 291              retina = cv::bioinspired::createRetina(inputImage.size());
 292
 293         // save default retina parameters file in order to let you see this and maybe modify it and reload using method "setup"
 294         retina->write("RetinaDefaultParameters.xml");
 295
 296                  // desactivate Magnocellular pathway processing (motion information extraction) since it is not usefull here
 297                  retina->activateMovingContoursProcessing(false);
 298
 299          // declare retina output buffers
 300          cv::Mat retinaOutput_parvo;
 301
 302          /////////////////////////////////////////////
 303          // prepare displays and interactions
 304          histogramClippingValue=0; // default value... updated with interface slider
 305
 306          std::string retinaInputCorrected("Retina input image (with cut edges histogram for basic pixels error avoidance)");
 307          cv::namedWindow(retinaInputCorrected,1);
 308          cv::createTrackbar("histogram edges clipping limit", "Retina input image (with cut edges histogram for basic pixels error avoidance)",&histogramClippingValue,50,callBack_rescaleGrayLevelMat);
 309
 310          std::string RetinaParvoWindow("Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping");
 311          cv::namedWindow(RetinaParvoWindow, 1);
 312          colorSaturationFactor=3;
 313          cv::createTrackbar("Color saturation", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &colorSaturationFactor,5,callback_saturateColors);
 314
 315          retinaHcellsGain=40;
 316          cv::createTrackbar("Hcells gain", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping",&retinaHcellsGain,100,callBack_updateRetinaParams);
 317
 318          localAdaptation_photoreceptors=197;
 319          localAdaptation_Gcells=190;
 320          cv::createTrackbar("Ph sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_photoreceptors,199,callBack_updateRetinaParams);
 321          cv::createTrackbar("Gcells sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_Gcells,199,callBack_updateRetinaParams);
 322
 323         std::string powerTransformedInput("EXR image with basic processing : 16bits=>8bits with gamma correction");
 324
 325          /////////////////////////////////////////////
 326          // apply default parameters of user interaction variables
 327          callBack_updateRetinaParams(1,NULL); // first call for default parameters setup
 328          callback_saturateColors(1, NULL);
 329
 330          // processing loop with stop condition
 331          currentFrameIndex=startFrameIndex;
 332          while(currentFrameIndex <= endFrameIndex)
 333          {
 334              loadNewFrame(inputImageNamePrototype, currentFrameIndex, false);
 335
 336              if (inputImage.empty())
 337              {
 338                 std::cout<<"Could not load new image (index = "<<currentFrameIndex<<"), program end"<<std::endl;
 339                 return -1;
 340              }
 341             // display input & process standard power transformation
 342             imshow("EXR image original image, 16bits=>8bits linear rescaling ", imageInputRescaled);
 343             cv::Mat gammaTransformedImage;
 344             cv::pow(imageInputRescaled, 1./5, gammaTransformedImage); // apply gamma curve: img = img ** (1./5)
 345             imshow(powerTransformedInput, gammaTransformedImage);
 346              // run retina filter
 347              retina->run(imageInputRescaled);
 348              // Retrieve and display retina output
 349              retina->getParvo(retinaOutput_parvo);
 350              cv::imshow(retinaInputCorrected, imageInputRescaled/255.f);
 351              cv::imshow(RetinaParvoWindow, retinaOutput_parvo);
 352              cv::waitKey(4);
 353             // jump to next frame
 354             ++currentFrameIndex;
 355          }
 356      }catch(cv::Exception e)
 357      {
 358          std::cerr<<"Error using Retina : "<<e.what()<<std::endl;
 359      }
 360
 361      // Program end message
 362      std::cout<<"Retina demo end"<<std::endl;
 363
 364      return 0;
 365  }