2 //============================================================================
3 // Name : HighDynamicRange_RetinaCompression.cpp
4 // Author : Alexandre Benoit (benoit.alexandre.vision@gmail.com)
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 //============================================================================
13 #include "opencv2/bioinspired.hpp" // retina based algorithms
14 #include "opencv2/imgproc.hpp" // cvCvtcolor function
15 #include "opencv2/highgui.hpp" // display
17 static void help(std::string errorMessage)
19 std::cout<<"Program init error : "<<errorMessage<<std::endl;
20 std::cout<<"\nProgram call procedure : ./OpenEXRimages_HDR_Retina_toneMapping [OpenEXR image to process]"<<std::endl;
21 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;
22 std::cout<<"\nExamples:"<<std::endl;
23 std::cout<<"\t-Image processing : ./OpenEXRimages_HDR_Retina_toneMapping memorial.exr"<<std::endl;
26 // simple procedure for 1D curve tracing
27 static void drawPlot(const cv::Mat curve, const std::string figureTitle, const int lowerLimit, const int upperLimit)
29 //std::cout<<"curve size(h,w) = "<<curve.size().height<<", "<<curve.size().width<<std::endl;
30 cv::Mat displayedCurveImage = cv::Mat::ones(200, curve.size().height, CV_8U);
32 cv::Mat windowNormalizedCurve;
33 normalize(curve, windowNormalizedCurve, 0, 200, cv::NORM_MINMAX, CV_32F);
35 displayedCurveImage = cv::Scalar::all(255); // set a white background
36 int binW = cvRound((double)displayedCurveImage.cols/curve.size().height);
38 for( int i = 0; i < curve.size().height; i++ )
39 rectangle( displayedCurveImage, cv::Point(i*binW, displayedCurveImage.rows),
40 cv::Point((i+1)*binW, displayedCurveImage.rows - cvRound(windowNormalizedCurve.at<float>(i))),
41 cv::Scalar::all(0), -1, 8, 0 );
42 rectangle( displayedCurveImage, cv::Point(0, 0),
43 cv::Point((lowerLimit)*binW, 200),
44 cv::Scalar::all(128), -1, 8, 0 );
45 rectangle( displayedCurveImage, cv::Point(displayedCurveImage.cols, 0),
46 cv::Point((upperLimit)*binW, 200),
47 cv::Scalar::all(128), -1, 8, 0 );
49 cv::imshow(figureTitle, displayedCurveImage);
52 * objective : get the gray level map of the input image and rescale it to the range [0-255]
54 static void rescaleGrayLevelMat(const cv::Mat &inputMat, cv::Mat &outputMat, const float histogramClippingLimit)
57 // adjust output matrix wrt the input size but single channel
58 std::cout<<"Input image rescaling with histogram edges cutting (in order to eliminate bad pixels created during the HDR image creation) :"<<std::endl;
59 //std::cout<<"=> image size (h,w,channels) = "<<inputMat.size().height<<", "<<inputMat.size().width<<", "<<inputMat.channels()<<std::endl;
60 //std::cout<<"=> pixel coding (nbchannel, bytes per channel) = "<<inputMat.elemSize()/inputMat.elemSize1()<<", "<<inputMat.elemSize1()<<std::endl;
62 // rescale between 0-255, keeping floating point values
63 cv::normalize(inputMat, outputMat, 0.0, 255.0, cv::NORM_MINMAX);
65 // extract a 8bit image that will be used for histogram edge cut
67 if (inputMat.channels()==1)
69 outputMat.convertTo(intGrayImage, CV_8U);
73 outputMat.convertTo(rgbIntImg, CV_8UC3);
74 cv::cvtColor(rgbIntImg, intGrayImage, cv::COLOR_BGR2GRAY);
77 // get histogram density probability in order to cut values under above edges limits (here 5-95%)... usefull for HDR pixel errors cancellation
80 calcHist(&intGrayImage, 1, 0, cv::Mat(), hist, 1, &histSize, 0);
81 cv::Mat normalizedHist;
82 normalize(hist, normalizedHist, 1, 0, cv::NORM_L1, CV_32F); // normalize histogram so that its sum equals 1
84 double min_val, max_val;
85 minMaxLoc(normalizedHist, &min_val, &max_val);
86 //std::cout<<"Hist max,min = "<<max_val<<", "<<min_val<<std::endl;
88 // compute density probability
89 cv::Mat denseProb=cv::Mat::zeros(normalizedHist.size(), CV_32F);
90 denseProb.at<float>(0)=normalizedHist.at<float>(0);
91 int histLowerLimit=0, histUpperLimit=0;
92 for (int i=1;i<normalizedHist.size().height;++i)
94 denseProb.at<float>(i)=denseProb.at<float>(i-1)+normalizedHist.at<float>(i);
95 //std::cout<<normalizedHist.at<float>(i)<<", "<<denseProb.at<float>(i)<<std::endl;
96 if ( denseProb.at<float>(i)<histogramClippingLimit)
98 if ( denseProb.at<float>(i)<1-histogramClippingLimit)
101 // deduce min and max admitted gray levels
102 float minInputValue = (float)histLowerLimit/histSize*255;
103 float maxInputValue = (float)histUpperLimit/histSize*255;
105 std::cout<<"=> Histogram limits "
106 <<"\n\t"<<histogramClippingLimit*100<<"% index = "<<histLowerLimit<<" => normalizedHist value = "<<denseProb.at<float>(histLowerLimit)<<" => input gray level = "<<minInputValue
107 <<"\n\t"<<(1-histogramClippingLimit)*100<<"% index = "<<histUpperLimit<<" => normalizedHist value = "<<denseProb.at<float>(histUpperLimit)<<" => input gray level = "<<maxInputValue
109 //drawPlot(denseProb, "input histogram density probability", histLowerLimit, histUpperLimit);
110 drawPlot(normalizedHist, "input histogram", histLowerLimit, histUpperLimit);
112 // rescale image range [minInputValue-maxInputValue] to [0-255]
113 outputMat-=minInputValue;
114 outputMat*=255.0/(maxInputValue-minInputValue);
115 // cut original histogram and back project to original image
116 cv::threshold( outputMat, outputMat, 255.0, 255.0, 2 ); //THRESH_TRUNC, clips values above 255
117 cv::threshold( outputMat, outputMat, 0.0, 0.0, 3 ); //THRESH_TOZERO, clips values under 0
120 // basic callback method for interface management
122 cv::Mat imageInputRescaled;
123 int histogramClippingValue;
124 static void callBack_rescaleGrayLevelMat(int, void*)
126 std::cout<<"Histogram clipping value changed, current value = "<<histogramClippingValue<<std::endl;
127 rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)(histogramClippingValue/100.0));
128 normalize(imageInputRescaled, imageInputRescaled, 0.0, 255.0, cv::NORM_MINMAX);
131 cv::Ptr<cv::bioinspired::Retina> retina;
132 int retinaHcellsGain;
133 int localAdaptation_photoreceptors, localAdaptation_Gcells;
134 static void callBack_updateRetinaParams(int, void*)
136 retina->setupOPLandIPLParvoChannel(true, true, (float)(localAdaptation_photoreceptors/200.0), 0.5f, 0.43f, (float)retinaHcellsGain, 1.f, 7.f, (float)(localAdaptation_Gcells/200.0));
139 int colorSaturationFactor;
140 static void callback_saturateColors(int, void*)
142 retina->setColorSaturation(true, (float)colorSaturationFactor);
145 int main(int argc, char* argv[]) {
147 std::cout<<"*********************************************************************************"<<std::endl;
148 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;
149 std::cout<<"* This retina model allows spatio-temporal image processing (applied on still images, video sequences)."<<std::endl;
150 std::cout<<"* This demo focuses demonstration of the dynamic compression capabilities of the model"<<std::endl;
151 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;
152 std::cout<<"* The retina model still have the following properties:"<<std::endl;
153 std::cout<<"* => It applies a spectral whithening (mid-frequency details enhancement)"<<std::endl;
154 std::cout<<"* => high frequency spatio-temporal noise reduction"<<std::endl;
155 std::cout<<"* => low frequency luminance to be reduced (luminance range compression)"<<std::endl;
156 std::cout<<"* => local logarithmic luminance compression allows details to be enhanced in low light conditions\n"<<std::endl;
157 std::cout<<"* for more information, reer to the following papers :"<<std::endl;
158 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;
159 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;
160 std::cout<<"* => reports comments/remarks at benoit.alexandre.vision@gmail.com"<<std::endl;
161 std::cout<<"* => more informations and papers at : http://sites.google.com/site/benoitalexandrevision/"<<std::endl;
162 std::cout<<"*********************************************************************************"<<std::endl;
163 std::cout<<"** WARNING : this sample requires OpenCV to be configured with OpenEXR support **"<<std::endl;
164 std::cout<<"*********************************************************************************"<<std::endl;
165 std::cout<<"*** You can use free tools to generate OpenEXR images from images sets : ***"<<std::endl;
166 std::cout<<"*** => 1. take a set of photos from the same viewpoint using bracketing ***"<<std::endl;
167 std::cout<<"*** => 2. generate an OpenEXR image with tools like qtpfsgui.sourceforge.net ***"<<std::endl;
168 std::cout<<"*** => 3. apply tone mapping with this program ***"<<std::endl;
169 std::cout<<"*********************************************************************************"<<std::endl;
171 // basic input arguments checking
174 help("bad number of parameter");
178 bool useLogSampling = !strcmp(argv[argc-1], "log"); // check if user wants retina log sampling processing
180 if (!strcmp(argv[argc-1], "fast"))
183 std::cout<<"Using fast method (no spectral whithning), adaptation of Meylan&al 2008 method"<<std::endl;
186 std::string inputImageName=argv[1];
188 //////////////////////////////////////////////////////////////////////////////
189 // checking input media type (still image, video file, live video acquisition)
190 std::cout<<"RetinaDemo: processing image "<<inputImageName<<std::endl;
191 // image processing case
192 // declare the retina input buffer... that will be fed differently in regard of the input media
193 inputImage = cv::imread(inputImageName, -1); // load image in RGB mode
194 std::cout<<"=> image size (h,w) = "<<inputImage.size().height<<", "<<inputImage.size().width<<std::endl;
195 if (!inputImage.total())
197 help("could not load image, program end");
200 // rescale between 0 and 1
201 normalize(inputImage, inputImage, 0.0, 1.0, cv::NORM_MINMAX);
202 cv::Mat gammaTransformedImage;
203 cv::pow(inputImage, 1./5, gammaTransformedImage); // apply gamma curve: img = img ** (1./5)
204 imshow("EXR image original image, 16bits=>8bits linear rescaling ", inputImage);
205 imshow("EXR image with basic processing : 16bits=>8bits with gamma correction", gammaTransformedImage);
206 if (inputImage.empty())
208 help("Input image could not be loaded, aborting");
212 //////////////////////////////////////////////////////////////////////////////
213 // Program start in a try/catch safety context (Retina may throw errors)
216 /* create a retina instance with default parameters setup, uncomment the initialisation you wanna test
217 * -> if the last parameter is 'log', then activate log sampling (favour foveal vision and subsamples peripheral vision)
221 retina = cv::bioinspired::createRetina(inputImage.size(),true, cv::bioinspired::RETINA_COLOR_BAYER, true, 2.0, 10.0);
223 else// -> else allocate "classical" retina :
224 retina = cv::bioinspired::createRetina(inputImage.size());
226 // create a fast retina tone mapper (Meyla&al algorithm)
227 std::cout<<"Allocating fast tone mapper..."<<std::endl;
228 //cv::Ptr<cv::RetinaFastToneMapping> fastToneMapper=createRetinaFastToneMapping(inputImage.size());
229 std::cout<<"Fast tone mapper allocated"<<std::endl;
231 // save default retina parameters file in order to let you see this and maybe modify it and reload using method "setup"
232 retina->write("RetinaDefaultParameters.xml");
234 // desactivate Magnocellular pathway processing (motion information extraction) since it is not usefull here
235 retina->activateMovingContoursProcessing(false);
237 // declare retina output buffers
238 cv::Mat retinaOutput_parvo;
240 /////////////////////////////////////////////
241 // prepare displays and interactions
242 histogramClippingValue=0; // default value... updated with interface slider
243 //inputRescaleMat = inputImage;
244 //outputRescaleMat = imageInputRescaled;
245 cv::namedWindow("Processing configuration",1);
246 cv::createTrackbar("histogram edges clipping limit", "Processing configuration",&histogramClippingValue,50,callBack_rescaleGrayLevelMat);
248 colorSaturationFactor=3;
249 cv::createTrackbar("Color saturation", "Processing configuration", &colorSaturationFactor,5,callback_saturateColors);
252 cv::createTrackbar("Hcells gain", "Processing configuration",&retinaHcellsGain,100,callBack_updateRetinaParams);
254 localAdaptation_photoreceptors=197;
255 localAdaptation_Gcells=190;
256 cv::createTrackbar("Ph sensitivity", "Processing configuration", &localAdaptation_photoreceptors,199,callBack_updateRetinaParams);
257 cv::createTrackbar("Gcells sensitivity", "Processing configuration", &localAdaptation_Gcells,199,callBack_updateRetinaParams);
260 /////////////////////////////////////////////
261 // apply default parameters of user interaction variables
262 rescaleGrayLevelMat(inputImage, imageInputRescaled, (float)histogramClippingValue/100);
263 retina->setColorSaturation(true,(float)colorSaturationFactor);
264 callBack_updateRetinaParams(1,NULL); // first call for default parameters setup
266 // processing loop with stop condition
267 bool continueProcessing=true;
268 while(continueProcessing)
273 retina->run(imageInputRescaled);
274 // Retrieve and display retina output
275 retina->getParvo(retinaOutput_parvo);
276 cv::imshow("Retina input image (with cut edges histogram for basic pixels error avoidance)", imageInputRescaled/255.0);
277 cv::imshow("Retina Parvocellular pathway output : 16bit=>8bit image retina tonemapping", retinaOutput_parvo);
278 cv::imwrite("HDRinput.jpg",imageInputRescaled/255.0);
279 cv::imwrite("RetinaToneMapping.jpg",retinaOutput_parvo);
283 // apply the simplified hdr tone mapping method
284 cv::Mat fastToneMappingOutput;
285 retina->applyFastToneMapping(imageInputRescaled, fastToneMappingOutput);
286 cv::imshow("Retina fast tone mapping output : 16bit=>8bit image retina tonemapping", fastToneMappingOutput);
288 /*cv::Mat fastToneMappingOutput_specificObject;
289 fastToneMapper->setup(3.f, 1.5f, 1.f);
290 fastToneMapper->applyFastToneMapping(imageInputRescaled, fastToneMappingOutput_specificObject);
291 cv::imshow("### Retina fast tone mapping output : 16bit=>8bit image retina tonemapping", fastToneMappingOutput_specificObject);
295 }catch(cv::Exception e)
297 std::cerr<<"Error using Retina : "<<e.what()<<std::endl;
300 // Program end message
301 std::cout<<"Retina demo end"<<std::endl;