OpenEXRimages_HDR_Retina_toneMapping_video.cpp 20.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365

//============================================================================
// Name        : OpenEXRimages_HDR_Retina_toneMapping_video.cpp
// Author      : Alexandre Benoit (benoit.alexandre.vision@gmail.com)
// Version     : 0.2
// Copyright   : Alexandre Benoit, LISTIC Lab, december 2011
// Description : HighDynamicRange retina tone mapping for image sequences with the help of the Gipsa/Listic's retina in C++, Ansi-style
// Known issues: the input OpenEXR sequences can have bad computed pixels that should be removed
//               => a simple method consists of cutting histogram edges (a slider for this on the UI is provided)
//               => however, in image sequences, this histogramm cut must be done in an elegant way from frame to frame... still not done...
//============================================================================

#include <iostream>
#include <stdio.h>
#include <cstring>

#include "opencv2/bioinspired.hpp" // retina based algorithms
#include "opencv2/imgproc.hpp" // cvCvtcolor function
#include "opencv2/highgui.hpp" // display

#ifndef _CRT_SECURE_NO_WARNINGS
# define _CRT_SECURE_NO_WARNINGS
#endif

static void help(std::string errorMessage)
{
    std::cout<<"Program init error : "<<errorMessage<<std::endl;
    std::cout<<"\nProgram call procedure : ./OpenEXRimages_HDR_Retina_toneMapping [OpenEXR image sequence to process] [OPTIONNAL start frame] [OPTIONNAL end frame]"<<std::endl;
    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;
    std::cout<<"\t\t => WARNING : image index number of digits cannot exceed 10"<<std::endl;
    std::cout<<"\t[start frame] : the starting frame tat should be considered"<<std::endl;
    std::cout<<"\t[end frame] : the ending frame tat should be considered"<<std::endl;
    std::cout<<"\nExamples:"<<std::endl;
    std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping_video memorial%3d.exr 20 45"<<std::endl;
    std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping_video memorial%3d.exr 20 45 log"<<std::endl;
    std::cout<<"\t ==> to process images from memorial020d.exr to memorial045d.exr"<<std::endl;

}

// simple procedure for 1D curve tracing
static void drawPlot(const cv::Mat curve, const std::string figureTitle, const int lowerLimit, const int upperLimit)
{
    //std::cout<<"curve size(h,w) = "<<curve.size().height<<", "<<curve.size().width<<std::endl;
    cv::Mat displayedCurveImage = cv::Mat::ones(200, curve.size().height, CV_8U);

    cv::Mat windowNormalizedCurve;
    normalize(curve, windowNormalizedCurve, 0, 200, cv::NORM_MINMAX, CV_32F);

    displayedCurveImage = cv::Scalar::all(255); // set a white background
    int binW = cvRound((double)displayedCurveImage.cols/curve.size().height);

    for( int i = 0; i < curve.size().height; i++ )
        rectangle( displayedCurveImage, cv::Point(i*binW, displayedCurveImage.rows),
                cv::Point((i+1)*binW, displayedCurveImage.rows - cvRound(windowNormalizedCurve.at<float>(i))),
                cv::Scalar::all(0), -1, 8, 0 );
    rectangle( displayedCurveImage, cv::Point(0, 0),
            cv::Point((lowerLimit)*binW, 200),
            cv::Scalar::all(128), -1, 8, 0 );
    rectangle( displayedCurveImage, cv::Point(displayedCurveImage.cols, 0),
            cv::Point((upperLimit)*binW, 200),
            cv::Scalar::all(128), -1, 8, 0 );

    cv::imshow(figureTitle, displayedCurveImage);
}

/*
 * 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
 */
static void rescaleGrayLevelMat(const cv::Mat &inputMat, cv::Mat &outputMat, const float histogramClippingLimit, const bool rescale0_255)
 {
     // adjust output matrix wrt the input size but single channel
     std::cout<<"Input image rescaling with histogram edges cutting (in order to eliminate bad pixels created during the HDR image creation) :"<<std::endl;
     //std::cout<<"=> image size (h,w,channels) = "<<inputMat.size().height<<", "<<inputMat.size().width<<", "<<inputMat.channels()<<std::endl;
     //std::cout<<"=> pixel coding (nbchannel, bytes per channel) = "<<inputMat.elemSize()/inputMat.elemSize1()<<", "<<inputMat.elemSize1()<<std::endl;

     // get min and max values to use afterwards if no 0-255 rescaling is used
     double maxInput, minInput, histNormRescalefactor=1.f;
     double histNormOffset=0.f;
     minMaxLoc(inputMat, &minInput, &maxInput);
     histNormRescalefactor=255.f/(maxInput-minInput);
     histNormOffset=minInput;
     std::cout<<"Hist max,min = "<<maxInput<<", "<<minInput<<" => scale, offset = "<<histNormRescalefactor<<", "<<histNormOffset<<std::endl;
     // rescale between 0-255, keeping floating point values
     cv::Mat normalisedImage;
     cv::normalize(inputMat, normalisedImage, 0.f, 255.f, cv::NORM_MINMAX);
     if (rescale0_255)
        normalisedImage.copyTo(outputMat);
     // extract a 8bit image that will be used for histogram edge cut
     cv::Mat intGrayImage;
     if (inputMat.channels()==1)
     {
         normalisedImage.convertTo(intGrayImage, CV_8U);
     }else
     {
         cv::Mat rgbIntImg;
         normalisedImage.convertTo(rgbIntImg, CV_8UC3);
         cvtColor(rgbIntImg, intGrayImage, cv::COLOR_BGR2GRAY);
     }

     // get histogram density probability in order to cut values under above edges limits (here 5-95%)... usefull for HDR pixel errors cancellation
     cv::Mat dst, hist;
     int histSize = 256;
     calcHist(&intGrayImage, 1, 0, cv::Mat(), hist, 1, &histSize, 0);
     cv::Mat normalizedHist;

     normalize(hist, normalizedHist, 1.f, 0.f, cv::NORM_L1, CV_32F); // normalize histogram so that its sum equals 1

     // compute density probability
     cv::Mat denseProb=cv::Mat::zeros(normalizedHist.size(), CV_32F);
     denseProb.at<float>(0)=normalizedHist.at<float>(0);
     int histLowerLimit=0, histUpperLimit=0;
     for (int i=1;i<normalizedHist.size().height;++i)
     {
         denseProb.at<float>(i)=denseProb.at<float>(i-1)+normalizedHist.at<float>(i);
         //std::cout<<normalizedHist.at<float>(i)<<", "<<denseProb.at<float>(i)<<std::endl;
         if ( denseProb.at<float>(i)<histogramClippingLimit)
             histLowerLimit=i;
         if ( denseProb.at<float>(i)<1.f-histogramClippingLimit)
             histUpperLimit=i;
     }
     // deduce min and max admitted gray levels
     float minInputValue = (float)histLowerLimit/histSize*255.f;
     float maxInputValue = (float)histUpperLimit/histSize*255.f;

     std::cout<<"=> Histogram limits "
             <<"\n\t"<<histogramClippingLimit*100.f<<"% index = "<<histLowerLimit<<" => normalizedHist value = "<<denseProb.at<float>(histLowerLimit)<<" => input gray level = "<<minInputValue
             <<"\n\t"<<(1.f-histogramClippingLimit)*100.f<<"% index = "<<histUpperLimit<<" => normalizedHist value = "<<denseProb.at<float>(histUpperLimit)<<" => input gray level = "<<maxInputValue
             <<std::endl;
     //drawPlot(denseProb, "input histogram density probability", histLowerLimit, histUpperLimit);
     drawPlot(normalizedHist, "input histogram", histLowerLimit, histUpperLimit);

    if(rescale0_255) // rescale between 0-255 if asked to
    {
        cv::threshold( outputMat, outputMat, maxInputValue, maxInputValue, 2 ); //THRESH_TRUNC, clips values above maxInputValue
        cv::threshold( outputMat, outputMat, minInputValue, minInputValue, 3 ); //THRESH_TOZERO, clips values under minInputValue
        // rescale image range [minInputValue-maxInputValue] to [0-255]
        outputMat-=minInputValue;
        outputMat*=255.f/(maxInputValue-minInputValue);
    }else
    {
        inputMat.copyTo(outputMat);
        // update threshold in the initial input image range
        maxInputValue=(float)((maxInputValue-255.f)/histNormRescalefactor+maxInput);
        minInputValue=(float)(minInputValue/histNormRescalefactor+minInput);
        std::cout<<"===> Input Hist clipping values (max,min) = "<<maxInputValue<<", "<<minInputValue<<std::endl;
        cv::threshold( outputMat, outputMat, maxInputValue, maxInputValue, 2 ); //THRESH_TRUNC, clips values above maxInputValue
        cv::threshold( outputMat, outputMat, minInputValue, minInputValue, 3 ); //
    }
 }

 // basic callback method for interface management
 cv::Mat inputImage;
 cv::Mat imageInputRescaled;
 float globalRescalefactor=1;
 cv::Scalar globalOffset=0;
 int histogramClippingValue;
 static void callBack_rescaleGrayLevelMat(int, void*)
 {
     std::cout<<"Histogram clipping value changed, current value = "<<histogramClippingValue<<std::endl;
    // rescale and process
    inputImage+=globalOffset;
    inputImage*=globalRescalefactor;
    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
    rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)histogramClippingValue/100.f, true);

 }

 cv::Ptr<cv::bioinspired::Retina> retina;
 int retinaHcellsGain;
 int localAdaptation_photoreceptors, localAdaptation_Gcells;
 static void callBack_updateRetinaParams(int, void*)
 {
     retina->setupOPLandIPLParvoChannel(true, true, (float)(localAdaptation_photoreceptors/200.0), 0.5f, 0.43f, (float)retinaHcellsGain, 1.f, 7.f, (float)(localAdaptation_Gcells/200.0));
 }

 int colorSaturationFactor;
 static void callback_saturateColors(int, void*)
 {
     retina->setColorSaturation(true, (float)colorSaturationFactor);
 }

// 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)
static void loadNewFrame(const std::string filenamePrototype, const int currentFileIndex, const bool firstTimeread)
{
     char *currentImageName=NULL;
    currentImageName = (char*)malloc(sizeof(char)*filenamePrototype.size()+10);

    // grab the first frame
    sprintf(currentImageName, filenamePrototype.c_str(), currentFileIndex);

     //////////////////////////////////////////////////////////////////////////////
     // checking input media type (still image, video file, live video acquisition)
     std::cout<<"RetinaDemo: reading image : "<<currentImageName<<std::endl;
     // image processing case
     // declare the retina input buffer... that will be fed differently in regard of the input media
     inputImage = cv::imread(currentImageName, -1); // load image in RGB mode
     std::cout<<"=> image size (h,w) = "<<inputImage.size().height<<", "<<inputImage.size().width<<std::endl;
     if (inputImage.empty())
     {
        help("could not load image, program end");
            return;;
         }

    // rescaling/histogram clipping stage
    // rescale between 0 and 1
    // 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
    double maxInput, minInput;
    minMaxLoc(inputImage, &minInput, &maxInput);
    std::cout<<"ORIGINAL IMAGE pixels values range (max,min) : "<<maxInput<<", "<<minInput<<std::endl;

    if (firstTimeread)
    {
        /* the first time, get the pixel values range and rougthly update scaling value
        in order to center values around 128 and getting a range close to [0-255],
        => actually using a little less in order to let some more flexibility in range evolves...
        */
        double maxInput1, minInput1;
        minMaxLoc(inputImage, &minInput1, &maxInput1);
        std::cout<<"FIRST IMAGE pixels values range (max,min) : "<<maxInput1<<", "<<minInput1<<std::endl;
        globalRescalefactor=(float)(50.0/(maxInput1-minInput1)); // less than 255 for flexibility... experimental value to be carefull about
        double channelOffset = -1.5*minInput;
        globalOffset= cv::Scalar(channelOffset, channelOffset, channelOffset, channelOffset);
    }
    // call the generic input image rescaling callback
    callBack_rescaleGrayLevelMat(1,NULL);
}

 int main(int argc, char* argv[]) {
     // welcome message
     std::cout<<"*********************************************************************************"<<std::endl;
     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;
     std::cout<<"* This retina model allows spatio-temporal image processing (applied on still images, video sequences)."<<std::endl;
     std::cout<<"* This demo focuses demonstration of the dynamic compression capabilities of the model"<<std::endl;
     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;
     std::cout<<"* The retina model still have the following properties:"<<std::endl;
     std::cout<<"* => It applies a spectral whithening (mid-frequency details enhancement)"<<std::endl;
     std::cout<<"* => high frequency spatio-temporal noise reduction"<<std::endl;
     std::cout<<"* => low frequency luminance to be reduced (luminance range compression)"<<std::endl;
     std::cout<<"* => local logarithmic luminance compression allows details to be enhanced in low light conditions\n"<<std::endl;
     std::cout<<"* for more information, reer to the following papers :"<<std::endl;
     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;
     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;
     std::cout<<"* => reports comments/remarks at benoit.alexandre.vision@gmail.com"<<std::endl;
     std::cout<<"* => more informations and papers at : http://sites.google.com/site/benoitalexandrevision/"<<std::endl;
     std::cout<<"*********************************************************************************"<<std::endl;
     std::cout<<"** WARNING : this sample requires OpenCV to be configured with OpenEXR support **"<<std::endl;
     std::cout<<"*********************************************************************************"<<std::endl;
     std::cout<<"*** You can use free tools to generate OpenEXR images from images sets   :    ***"<<std::endl;
     std::cout<<"*** =>  1. take a set of photos from the same viewpoint using bracketing      ***"<<std::endl;
     std::cout<<"*** =>  2. generate an OpenEXR image with tools like qtpfsgui.sourceforge.net ***"<<std::endl;
     std::cout<<"*** =>  3. apply tone mapping with this program                               ***"<<std::endl;
     std::cout<<"*********************************************************************************"<<std::endl;

     // basic input arguments checking
     if (argc<4)
     {
         help("bad number of parameter");
         return -1;
     }

     bool useLogSampling = !strcmp(argv[argc-1], "log"); // check if user wants retina log sampling processing

     int startFrameIndex=0, endFrameIndex=0, currentFrameIndex=0;
     sscanf(argv[2], "%d", &startFrameIndex);
     sscanf(argv[3], "%d", &endFrameIndex);
     std::string inputImageNamePrototype(argv[1]);

     //////////////////////////////////////////////////////////////////////////////
     // checking input media type (still image, video file, live video acquisition)
     std::cout<<"RetinaDemo: setting up system with first image..."<<std::endl;
     loadNewFrame(inputImageNamePrototype, startFrameIndex, true);

     if (inputImage.empty())
     {
        help("could not load image, program end");
            return -1;
         }

     //////////////////////////////////////////////////////////////////////////////
     // Program start in a try/catch safety context (Retina may throw errors)
     try
     {
         /* create a retina instance with default parameters setup, uncomment the initialisation you wanna test
          * -> if the last parameter is 'log', then activate log sampling (favour foveal vision and subsamples peripheral vision)
          */
         if (useLogSampling)
                {
                     retina = cv::bioinspired::createRetina(inputImage.size(),true, cv::bioinspired::RETINA_COLOR_BAYER, true, 2.0, 10.0);
                 }
         else// -> else allocate "classical" retina :
             retina = cv::bioinspired::createRetina(inputImage.size());

        // save default retina parameters file in order to let you see this and maybe modify it and reload using method "setup"
        retina->write("RetinaDefaultParameters.xml");

                 // desactivate Magnocellular pathway processing (motion information extraction) since it is not usefull here
                 retina->activateMovingContoursProcessing(false);

         // declare retina output buffers
         cv::Mat retinaOutput_parvo;

         /////////////////////////////////////////////
         // prepare displays and interactions
         histogramClippingValue=0; // default value... updated with interface slider

         std::string retinaInputCorrected("Retina input image (with cut edges histogram for basic pixels error avoidance)");
         cv::namedWindow(retinaInputCorrected,1);
         cv::createTrackbar("histogram edges clipping limit", "Retina input image (with cut edges histogram for basic pixels error avoidance)",&histogramClippingValue,50,callBack_rescaleGrayLevelMat);

         std::string RetinaParvoWindow("Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping");
         cv::namedWindow(RetinaParvoWindow, 1);
         colorSaturationFactor=3;
         cv::createTrackbar("Color saturation", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &colorSaturationFactor,5,callback_saturateColors);

         retinaHcellsGain=40;
         cv::createTrackbar("Hcells gain", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping",&retinaHcellsGain,100,callBack_updateRetinaParams);

         localAdaptation_photoreceptors=197;
         localAdaptation_Gcells=190;
         cv::createTrackbar("Ph sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_photoreceptors,199,callBack_updateRetinaParams);
         cv::createTrackbar("Gcells sensitivity", "Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", &localAdaptation_Gcells,199,callBack_updateRetinaParams);

        std::string powerTransformedInput("EXR image with basic processing : 16bits=>8bits with gamma correction");

         /////////////////////////////////////////////
         // apply default parameters of user interaction variables
         callBack_updateRetinaParams(1,NULL); // first call for default parameters setup
         callback_saturateColors(1, NULL);

         // processing loop with stop condition
         currentFrameIndex=startFrameIndex;
         while(currentFrameIndex <= endFrameIndex)
         {
             loadNewFrame(inputImageNamePrototype, currentFrameIndex, false);

             if (inputImage.empty())
             {
                std::cout<<"Could not load new image (index = "<<currentFrameIndex<<"), program end"<<std::endl;
                return -1;
             }
            // display input & process standard power transformation
            imshow("EXR image original image, 16bits=>8bits linear rescaling ", imageInputRescaled);
            cv::Mat gammaTransformedImage;
            cv::pow(imageInputRescaled, 1./5, gammaTransformedImage); // apply gamma curve: img = img ** (1./5)
            imshow(powerTransformedInput, gammaTransformedImage);
             // run retina filter
             retina->run(imageInputRescaled);
             // Retrieve and display retina output
             retina->getParvo(retinaOutput_parvo);
             cv::imshow(retinaInputCorrected, imageInputRescaled/255.f);
             cv::imshow(RetinaParvoWindow, retinaOutput_parvo);
             cv::waitKey(4);
            // jump to next frame
            ++currentFrameIndex;
         }
     }catch(cv::Exception e)
     {
         std::cerr<<"Error using Retina : "<<e.what()<<std::endl;
     }

     // Program end message
     std::cout<<"Retina demo end"<<std::endl;

     return 0;
 }