ITKeyword,专注技术干货聚合推荐

注册 | 登录

图像增强算法四种,图示与源码,包括retinex(ssr、msr、msrcr)和一种混合算法

android_asp 分享于 2013-07-24

2020腾讯云8月秒杀活动,优惠非常大!(领取2860元代金券),
地址https://cloud.tencent.com/act/cps/redirect?redirect=1040

2020阿里云最低价产品入口,含代金券(新老用户有优惠),
地址https://www.aliyun.com/minisite/goods

两组图像:左边较暗,右边较亮

第一行是原图像,他们下面是用四种算法处理的结果

依次为:

1.一种混合算法

2.msr,multi-scale retinex

3.msrcr,multi-scale retinex with color restoration

4.ssr,single scale retinex

 

 

 

 

 

 

 

源码,retinex算法的三种,其源码是国外一个研究生的毕设项目

头文件:

[cpp]  view plain copy
  1. /* 
  2.  * Copyright (c) 2006, Douglas Gray (dgray@soe.ucsc.edu, dr.de3ug@gmail.com) 
  3.  * All rights reserved. 
  4.  * 
  5.  * Redistribution and use in source and binary forms, with or without 
  6.  * modification, are permitted provided that the following conditions are met: 
  7.  *     * Redistributions of source code must retain the above copyright 
  8.  *       notice, this list of conditions and the following disclaimer. 
  9.  *     * Redistributions in binary form must reproduce the above copyright 
  10.  *       notice, this list of conditions and the following disclaimer in the 
  11.  *       documentation and/or other materials provided with the distribution. 
  12.  *     * Neither the name of the <organization> nor the 
  13.  *       names of its contributors may be used to endorse or promote products 
  14.  *       derived from this software without specific prior written permission. 
  15.  * 
  16.  * THIS SOFTWARE IS PROVIDED BY Douglas Gray ``AS IS'' AND ANY 
  17.  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
  18.  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
  19.  * DISCLAIMED. IN NO EVENT SHALL <copyright holder> BE LIABLE FOR ANY 
  20.  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
  21.  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 
  22.  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 
  23.  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
  24.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
  25.  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
  26.  */  
  27. #pragma once  
  28.   
  29. #include "cv.h"  
  30.   
  31. extern double* CreateKernel(double sigma);  
  32. extern int* CreateFastKernel(double sigma);  
  33.   
  34. extern void FilterGaussian(IplImage* img, double sigma);  
  35. extern void FastFilter(IplImage *img, double sigma);  
  36.   
  37. extern void Retinex  
  38. (IplImage *img, double sigma, int gain = 128, int offset = 128);  
  39.   
  40. extern void MultiScaleRetinex  
  41. (IplImage *img, int scales, double *weights, double *sigmas, int gain = 128, int offset = 128);  
  42.   
  43. extern void MultiScaleRetinexCR  
  44. (IplImage *img, int scales, double *weights, double *sigmas, int gain = 128, int offset = 128,  
  45.  double restoration_factor = 6, double color_gain = 2);  

实现:

[cpp]  view plain copy
  1. /* 
  2.  * Copyright (c) 2006, Douglas Gray (dgray@soe.ucsc.edu, dr.de3ug@gmail.com) 
  3.  * All rights reserved. 
  4.  * 
  5.  * Redistribution and use in source and binary forms, with or without 
  6.  * modification, are permitted provided that the following conditions are met: 
  7.  *     * Redistributions of source code must retain the above copyright 
  8.  *       notice, this list of conditions and the following disclaimer. 
  9.  *     * Redistributions in binary form must reproduce the above copyright 
  10.  *       notice, this list of conditions and the following disclaimer in the 
  11.  *       documentation and/or other materials provided with the distribution. 
  12.  *     * Neither the name of the <organization> nor the 
  13.  *       names of its contributors may be used to endorse or promote products 
  14.  *       derived from this software without specific prior written permission. 
  15.  * 
  16.  * THIS SOFTWARE IS PROVIDED BY Douglas Gray ``AS IS'' AND ANY 
  17.  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
  18.  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
  19.  * DISCLAIMED. IN NO EVENT SHALL <copyright holder> BE LIABLE FOR ANY 
  20.  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
  21.  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 
  22.  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 
  23.  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
  24.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
  25.  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
  26.  */  
  27. #include "retinex.h"  
  28.   
  29. #include <math.h>  
  30.   
  31. //#define USE_EXACT_SIGMA  
  32.   
  33.   
  34. #define pc(image, x, y, c) image->imageData[(image->widthStep * y) + (image->nChannels * x) + c]  
  35.   
  36. #define INT_PREC 1024.0  
  37. #define INT_PREC_BITS 10  
  38.   
  39. inline double int2double(int x) { return (double)x / INT_PREC; }  
  40. inline int double2int(double x) { return (int)(x * INT_PREC + 0.5); }  
  41.   
  42. inline int int2smallint(int x) { return (x >> INT_PREC_BITS); }  
  43. inline int int2bigint(int x) { return (x << INT_PREC_BITS); }  
  44.   
  45. //  
  46. // CreateKernel  
  47. //  
  48. // Summary:  
  49. // Creates a normalized 1 dimensional gaussian kernel.  
  50. //  
  51. // Arguments:  
  52. // sigma - the standard deviation of the gaussian kernel.  
  53. //  
  54. // Returns:  
  55. // double* - an array of values of length ((6*sigma)/2) * 2 + 1.  
  56. //  
  57. // Note:  
  58. // Caller is responsable for deleting the kernel.  
  59. //  
  60. double*  
  61. CreateKernel(double sigma)  
  62. {  
  63.     int i, x, filter_size;  
  64.     double* filter;  
  65.     double sum;  
  66.   
  67.     // Reject unreasonable demands  
  68.     if ( sigma > 200 ) sigma = 200;  
  69.   
  70.     // get needed filter size (enforce oddness)  
  71.     filter_size = (int)floor(sigma*6) / 2;  
  72.     filter_size = filter_size * 2 + 1;  
  73.   
  74.     // Allocate kernel space  
  75.     filter = new double[filter_size];  
  76.   
  77.     // Calculate exponential  
  78.     sum = 0;  
  79.     for (i = 0; i < filter_size; i++) {  
  80.         x = i - (filter_size / 2);  
  81.         filter[i] = exp( -(x*x) / (2*sigma*sigma) );  
  82.   
  83.         sum += filter[i];  
  84.     }  
  85.   
  86.     // Normalize  
  87.     for (i = 0, x; i < filter_size; i++)  
  88.         filter[i] /= sum;  
  89.   
  90.     return filter;  
  91. }  
  92.   
  93. //  
  94. // CreateFastKernel  
  95. //  
  96. // Summary:  
  97. // Creates a faster gaussian kernal using integers that  
  98. // approximate floating point (leftshifted by 8 bits)  
  99. //  
  100. // Arguments:  
  101. // sigma - the standard deviation of the gaussian kernel.  
  102. //  
  103. // Returns:  
  104. // int* - an array of values of length ((6*sigma)/2) * 2 + 1.  
  105. //  
  106. // Note:  
  107. // Caller is responsable for deleting the kernel.  
  108. //  
  109.   
  110. int*  
  111. CreateFastKernel(double sigma)  
  112. {  
  113.     double* fp_kernel;  
  114.     int* kernel;  
  115.     int i, filter_size;  
  116.       
  117.     // Reject unreasonable demands  
  118.     if ( sigma > 200 ) sigma = 200;  
  119.   
  120.     // get needed filter size (enforce oddness)  
  121.     filter_size = (int)floor(sigma*6) / 2;  
  122.     filter_size = filter_size * 2 + 1;  
  123.   
  124.     // Allocate kernel space  
  125.     kernel = new int[filter_size];  
  126.   
  127.     fp_kernel = CreateKernel(sigma);  
  128.   
  129.     for (i = 0; i < filter_size; i++)  
  130.         kernel[i] = double2int(fp_kernel[i]);  
  131.   
  132.     delete fp_kernel;  
  133.   
  134.     return kernel;  
  135. }  
  136.   
  137.   
  138. //  
  139. // FilterGaussian  
  140. //  
  141. // Summary:  
  142. // Performs a gaussian convolution for a value of sigma that is equal  
  143. // in both directions.  
  144. //  
  145. // Arguments:  
  146. // img - the image to be filtered in place.  
  147. // sigma - the standard deviation of the gaussian kernel to use.  
  148. //  
  149. void  
  150. FilterGaussian(IplImage* img, double sigma)  
  151. {  
  152.     int i, j, k, source, filter_size;  
  153.     int* kernel;  
  154.     IplImage* temp;  
  155.     int v1, v2, v3;  
  156.   
  157.     // Reject unreasonable demands  
  158.     if ( sigma > 200 ) sigma = 200;  
  159.   
  160.     // get needed filter size (enforce oddness)  
  161.     filter_size = (int)floor(sigma*6) / 2;  
  162.     filter_size = filter_size * 2 + 1;  
  163.   
  164.     kernel = CreateFastKernel(sigma);  
  165.   
  166.     temp = cvCreateImage(cvSize(img->width, img->height), img->depth, img->nChannels);  
  167.       
  168.     // filter x axis  
  169.     for (j = 0; j < temp->height; j++)  
  170.     for (i = 0; i < temp->width; i++) {  
  171.   
  172.         // inner loop has been unrolled  
  173.   
  174.         v1 = v2 = v3 = 0;  
  175.         for (k = 0; k < filter_size; k++) {  
  176.               
  177.             source = i + filter_size / 2 - k;  
  178.       
  179.             if (source < 0) source *= -1;  
  180.             if (source > img->width - 1) source = 2*(img->width - 1) - source;  
  181.   
  182.             v1 += kernel[k] * (unsigned char)pc(img, source, j, 0);  
  183.             if (img->nChannels == 1) continue;  
  184.             v2 += kernel[k] * (unsigned char)pc(img, source, j, 1);  
  185.             v3 += kernel[k] * (unsigned char)pc(img, source, j, 2);  
  186.   
  187.         }  
  188.   
  189.         // set value and move on  
  190.         pc(temp, i, j, 0) = (char)int2smallint(v1);  
  191.         if (img->nChannels == 1) continue;  
  192.         pc(temp, i, j, 1) = (char)int2smallint(v2);  
  193.         pc(temp, i, j, 2) = (char)int2smallint(v3);  
  194.   
  195.     }  
  196.       
  197.     // filter y axis  
  198.     for (j = 0; j < img->height; j++)  
  199.     for (i = 0; i < img->width; i++) {  
  200.   
  201.         v1 = v2 = v3 = 0;  
  202.         for (k = 0; k < filter_size; k++) {  
  203.               
  204.             source = j + filter_size / 2 - k;  
  205.       
  206.             if (source < 0) source *= -1;  
  207.             if (source > temp->height - 1) source = 2*(temp->height - 1) - source;  
  208.   
  209.             v1 += kernel[k] * (unsigned char)pc(temp, i, source, 0);  
  210.             if (img->nChannels == 1) continue;  
  211.             v2 += kernel[k] * (unsigned char)pc(temp, i, source, 1);  
  212.             v3 += kernel[k] * (unsigned char)pc(temp, i, source, 2);  
  213.   
  214.         }  
  215.   
  216.         // set value and move on  
  217.         pc(img, i, j, 0) = (char)int2smallint(v1);  
  218.         if (img->nChannels == 1) continue;  
  219.         pc(img, i, j, 1) = (char)int2smallint(v2);  
  220.         pc(img, i, j, 2) = (char)int2smallint(v3);  
  221.   
  222.     }  
  223.   
  224.   
  225.     cvReleaseImage( &temp );  
  226.   
  227.     delete kernel;  
  228.   
  229. }  
  230.   
  231. //  
  232. // FastFilter  
  233. //  
  234. // Summary:  
  235. // Performs gaussian convolution of any size sigma very fast by using  
  236. // both image pyramids and seperable filters.  Recursion is used.  
  237. //  
  238. // Arguments:  
  239. // img - an IplImage to be filtered in place.  
  240. //  
  241. void  
  242. FastFilter(IplImage *img, double sigma)  
  243. {  
  244.     int filter_size;  
  245.   
  246.     // Reject unreasonable demands  
  247.     if ( sigma > 200 ) sigma = 200;  
  248.   
  249.     // get needed filter size (enforce oddness)  
  250.     filter_size = (int)floor(sigma*6) / 2;  
  251.     filter_size = filter_size * 2 + 1;  
  252.   
  253.     // If 3 sigma is less than a pixel, why bother (ie sigma < 2/3)  
  254.     if(filter_size < 3) return;  
  255.   
  256.     // Filter, or downsample and recurse  
  257.     if (filter_size < 10) {  
  258.   
  259. #ifdef USE_EXACT_SIGMA  
  260.         FilterGaussian(img, sigma)  
  261. #else  
  262.         cvSmooth( img, img, CV_GAUSSIAN, filter_size, filter_size );  
  263. #endif  
  264.               
  265.     }  
  266.     else {  
  267.         if (img->width < 2 || img->height < 2) return;  
  268.   
  269.         IplImage* sub_img = cvCreateImage(cvSize(img->width / 2, img->height / 2), img->depth, img->nChannels);  
  270.   
  271.         cvPyrDown( img, sub_img );  
  272.   
  273.         FastFilter( sub_img, sigma / 2.0 );  
  274.   
  275.         cvResize( sub_img, img, CV_INTER_LINEAR );  
  276.   
  277.         cvReleaseImage( &sub_img );  
  278.     }  
  279.   
  280. }  
  281.   
  282. //  
  283. // Retinex  
  284. //  
  285. // Summary:  
  286. // Basic retinex restoration.  The image and a filtered image are converted  
  287. // to the log domain and subtracted.  
  288. //  
  289. // Arguments:  
  290. // img - an IplImage to be enhanced in place.  
  291. // sigma - the standard deviation of the gaussian kernal used to filter.  
  292. // gain - the factor by which to scale the image back into visable range.  
  293. // offset - an offset similar to the gain.  
  294. //  
  295. void  
  296. Retinex(IplImage *img, double sigma, int gain, int offset)  
  297. {  
  298.     IplImage *A, *fA, *fB, *fC;  
  299.   
  300.     // Initialize temp images  
  301.     fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  302.     fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  303.     fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  304.   
  305.     // Compute log image  
  306.     cvConvert( img, fA );  
  307.     cvLog( fA, fB );  
  308.       
  309.     // Compute log of blured image  
  310.     A = cvCloneImage( img );  
  311.     FastFilter( A, sigma );  
  312.     cvConvert( A, fA );  
  313.     cvLog( fA, fC );  
  314.   
  315.     // Compute difference  
  316.     cvSub( fB, fC, fA );  
  317.   
  318.     // Restore  
  319.     cvConvertScale( fA, img, gain, offset);  
  320.   
  321.     // Release temp images  
  322.     cvReleaseImage( &A );  
  323.     cvReleaseImage( &fA );  
  324.     cvReleaseImage( &fB );  
  325.     cvReleaseImage( &fC );  
  326.   
  327. }  
  328.   
  329. //  
  330. // MultiScaleRetinex  
  331. //  
  332. // Summary:  
  333. // Multiscale retinex restoration.  The image and a set of filtered images are  
  334. // converted to the log domain and subtracted from the original with some set  
  335. // of weights. Typicaly called with three equaly weighted scales of fine,  
  336. // medium and wide standard deviations.  
  337. //  
  338. // Arguments:  
  339. // img - an IplImage to be enhanced in place.  
  340. // sigma - the standard deviation of the gaussian kernal used to filter.  
  341. // gain - the factor by which to scale the image back into visable range.  
  342. // offset - an offset similar to the gain.  
  343. //  
  344. void  
  345. MultiScaleRetinex(IplImage *img, int scales, double *weights, double *sigmas, int gain, int offset)  
  346. {  
  347.     int i;  
  348.     double weight;  
  349.     IplImage *A, *fA, *fB, *fC;  
  350.   
  351.     // Initialize temp images  
  352.     fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  353.     fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  354.     fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  355.   
  356.   
  357.     // Compute log image  
  358.     cvConvert( img, fA );  
  359.     cvLog( fA, fB );  
  360.   
  361.     // Normalize according to given weights  
  362.     for (i = 0, weight = 0; i < scales; i++)  
  363.         weight += weights[i];  
  364.   
  365.     if (weight != 1.0) cvScale( fB, fB, weight );  
  366.   
  367.     // Filter at each scale  
  368.     for (i = 0; i < scales; i++) {  
  369.         A = cvCloneImage( img );  
  370.         FastFilter( A, sigmas[i] );  
  371.       
  372.         cvConvert( A, fA );  
  373.         cvLog( fA, fC );  
  374.         cvReleaseImage( &A );  
  375.   
  376.         // Compute weighted difference  
  377.         cvScale( fC, fC, weights[i] );  
  378.         cvSub( fB, fC, fB );  
  379.     }  
  380.   
  381.     // Restore  
  382.     cvConvertScale( fB, img, gain, offset);  
  383.   
  384.     // Release temp images  
  385.     cvReleaseImage( &fA );  
  386.     cvReleaseImage( &fB );  
  387.     cvReleaseImage( &fC );  
  388. }  
  389.   
  390. //  
  391. // MultiScaleRetinexCR  
  392. //  
  393. // Summary:  
  394. // Multiscale retinex restoration with color restoration.  The image and a set of  
  395. // filtered images are converted to the log domain and subtracted from the  
  396. // original with some set of weights. Typicaly called with three equaly weighted  
  397. // scales of fine, medium and wide standard deviations. A color restoration weight  
  398. // is then applied to each color channel.  
  399. //  
  400. // Arguments:  
  401. // img - an IplImage to be enhanced in place.  
  402. // sigma - the standard deviation of the gaussian kernal used to filter.  
  403. // gain - the factor by which to scale the image back into visable range.  
  404. // offset - an offset similar to the gain.  
  405. // restoration_factor - controls the non-linearaty of the color restoration.  
  406. // color_gain - controls the color restoration gain.  
  407. //  
  408. void  
  409. MultiScaleRetinexCR(IplImage *img, int scales, double *weights, double *sigmas,  
  410.                     int gain, int offset, double restoration_factor, double color_gain)  
  411. {  
  412.     int i;  
  413.     double weight;  
  414.     IplImage *A, *B, *C, *fA, *fB, *fC, *fsA, *fsB, *fsC, *fsD, *fsE, *fsF;  
  415.   
  416.     // Initialize temp images  
  417.     fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  418.     fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  419.     fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  420.     fsA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  421.     fsB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  422.     fsC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  423.     fsD = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  424.     fsE = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  425.     fsF = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  426.   
  427.     // Compute log image  
  428.     cvConvert( img, fB );  
  429.     cvLog( fB, fA );  
  430.   
  431.     // Normalize according to given weights  
  432.     for (i = 0, weight = 0; i < scales; i++)  
  433.         weight += weights[i];  
  434.   
  435.     if (weight != 1.0) cvScale( fA, fA, weight );  
  436.   
  437.     // Filter at each scale  
  438.     for (i = 0; i < scales; i++) {  
  439.         A = cvCloneImage( img );  
  440.         FastFilter( A, sigmas[i] );  
  441.       
  442.         cvConvert( A, fB );  
  443.         cvLog( fB, fC );  
  444.         cvReleaseImage( &A );  
  445.   
  446.         // Compute weighted difference  
  447.         cvScale( fC, fC, weights[i] );  
  448.         cvSub( fA, fC, fA );  
  449.     }  
  450.   
  451.     // Color restoration  
  452.     if (img->nChannels > 1) {  
  453.         A = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
  454.         B = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
  455.         C = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
  456.           
  457.         // Divide image into channels, convert and store sum  
  458.         cvCvtPixToPlane( img, A, B, C, NULL );  
  459.   
  460.         cvConvert( A, fsA );  
  461.         cvConvert( B, fsB );  
  462.         cvConvert( C, fsC );  
  463.   
  464.         cvReleaseImage( &A );  
  465.         cvReleaseImage( &B );  
  466.         cvReleaseImage( &C );  
  467.   
  468.         // Sum components  
  469.         cvAdd( fsA, fsB, fsD );  
  470.         cvAdd( fsD, fsC, fsD );  
  471.   
  472.         // Normalize weights  
  473.         cvDiv( fsA, fsD, fsA, restoration_factor);  
  474.         cvDiv( fsB, fsD, fsB, restoration_factor);  
  475.         cvDiv( fsC, fsD, fsC, restoration_factor);  
  476.   
  477.         cvConvertScale( fsA, fsA, 1, 1 );  
  478.         cvConvertScale( fsB, fsB, 1, 1 );  
  479.         cvConvertScale( fsC, fsC, 1, 1 );  
  480.   
  481.         // Log weights  
  482.         cvLog( fsA, fsA );  
  483.         cvLog( fsB, fsB );  
  484.         cvLog( fsC, fsC );  
  485.   
  486.         // Divide retinex image, weight accordingly and recombine  
  487.         cvCvtPixToPlane( fA, fsD, fsE, fsF, NULL );  
  488.   
  489.         cvMul( fsD, fsA, fsD, color_gain);  
  490.         cvMul( fsE, fsB, fsE, color_gain );  
  491.         cvMul( fsF, fsC, fsF, color_gain );  
  492.   
  493.         cvCvtPlaneToPix( fsD, fsE, fsF, NULL, fA );  
  494.     }  
  495.   
  496.     // Restore  
  497.     cvConvertScale( fA, img, gain, offset);  
  498.   
  499.     // Release temp images  
  500.     cvReleaseImage( &fA );  
  501.     cvReleaseImage( &fB );  
  502.     cvReleaseImage( &fC );  
  503.     cvReleaseImage( &fsA );  
  504.     cvReleaseImage( &fsB );  
  505.     cvReleaseImage( &fsC );  
  506.     cvReleaseImage( &fsD );  
  507.     cvReleaseImage( &fsE );  
  508.     cvReleaseImage( &fsF );  
  509. }  
 

 

这种混合算法代码:

我参考一个matlab代码写的 

c版

[cpp]  view plain copy
  1. // author : onezeros.lee@gmail.com  
  2. // data   : 4/20/2011   
  3. /* 
  4. %%%%%%%%%%%%%%%RGB normalization%%%%%%%%%%%%%%%%%%%%%% 
  5. %its cascaded implementation of 1 section of paper "A FAST SKIN REGION DETECTOR" by 
  6. %Phil Chen, Dr.Christos Greecos 
  7. %and 
  8. %section 2.1 of paper "Simple and accurate face detection in color images" by 
  9. %YUI TING PAI et al 
  10. % Coding by Madhava.S.Bhat, Dept of Electronics an communication 
  11. %Dr.Ambedkar Institute of Technology, Bangalore 
  12. %madhava.s@dr-ait.org 
  13.  */  
  14.   
  15. #include <cv.h>  
  16. #include <cxcore.h>  
  17. #include <highgui.h>  
  18.   
  19. #include <iostream>  
  20. #include <algorithm>  
  21. using namespace std;  
  22.   
  23. void cvCompensationGlobal1(IplImage* _src,IplImage* dst)  
  24. {  
  25.     static const int B=0;  
  26.     static const int G=1;  
  27.     static const int R=2;  
  28.   
  29.     int width=_src->width;  
  30.     int height=_src->height;  
  31.     //double  
  32.     IplImage* src=cvCreateImage(cvGetSize(_src),IPL_DEPTH_64F,3);  
  33.     CvScalar minV=cvScalar(255,255,255,0);  
  34.     for (int h=0;h<height;h++) {  
  35.         unsigned char* p=(unsigned char*)_src->imageData+h*_src->widthStep;  
  36.         double* pc=(double*)(src->imageData+h*src->widthStep);  
  37.         for (int w=0;w<width;w++) {  
  38.             for (int i=0;i<3;i++) {  
  39.                 *pc=*p;  
  40.                 if (minV.val[i]>*pc) {  
  41.                     minV.val[i]=*pc;  
  42.                 }  
  43.                 pc++;  
  44.                 p++;  
  45.             }  
  46.         }  
  47.     }  
  48.     cvSubS(src,minV,src);  
  49.   
  50.     int blackNum=0;  
  51.     double total=0;  
  52.     double acc[3]={0};  
  53.     for (int h=0;h<height;h++) {  
  54.         double* p=(double*)(src->imageData+h*src->widthStep);  
  55.         for (int w=0;w<width;w++) {  
  56.             if (p[0]<0.001&&p[1]<0.001&&p[2]<0.001) {  
  57.                 blackNum++;  
  58.                 p+=3;  
  59.                 continue;  
  60.             }  
  61.             double a=p[R];  
  62.             double b=p[R];  
  63.             if (p[B]<a) {  
  64.                 a=p[B];  
  65.             }else {  
  66.                 b=p[B];  
  67.             }  
  68.             if (p[G]<a) {  
  69.                 a=p[G];  
  70.             }else {  
  71.                 b=p[G];  
  72.             }  
  73.             total+=a+b;  
  74.   
  75.             for (int i=0;i<3;i++) {  
  76.                 acc[i]+=p[i];  
  77.             }  
  78.             p+=3;  
  79.         }  
  80.     }  
  81.     double avgT=total/(2*(width*height-blackNum));  
  82.     CvScalar avg;  
  83.     IplImage* rgb[3];  
  84.     for (int i=0;i<3;i++) {  
  85.         rgb[i]=cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
  86.         avg.val[i]=(double)acc[i]/(width*height);  
  87.         avg.val[i]=avgT/avg.val[i];  
  88.     }  
  89.     cvSplit(src,rgb[0],rgb[1],rgb[2],0);  
  90.     for (int i=0;i<3;i++) {  
  91.         cvScale(rgb[i],rgb[i],avg.val[i]);//bigger than 255  
  92.     }  
  93.     cvMerge(rgb[0],rgb[1],rgb[2],0,src);  
  94.   
  95.     //y component only  
  96.     IplImage* y=cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
  97.     for (int h=0;h<height;h++) {  
  98.         double* psrc=(double*)(src->imageData+h*src->widthStep);  
  99.         double* py=(double*)(y->imageData+h*y->widthStep);  
  100.         for (int w=0;w<width;w++) {  
  101.             *py++=psrc[R]*0.299+psrc[G]*0.587+psrc[B]*0.114;  
  102.             psrc+=3;  
  103.         }  
  104.     }  
  105.     double maxY=0;  
  106.     double minY=0;  
  107.     cvMinMaxLoc(y,&minY,&maxY);  
  108.       
  109.     double sumY=0;  
  110.     //scale   
  111.     for (int h=0;h<height;h++) {  
  112.         double* p=(double*)(y->imageData+h*y->widthStep);  
  113.         for (int w=0;w<width;w++) {  
  114.             *p=(*p-minY)/(maxY-minY);  
  115.             sumY+=*p;  
  116.             p++;  
  117.         }  
  118.     }  
  119.     sumY*=255;  
  120.     sumY/=width*height;  
  121.       
  122.     double scale=1.;  
  123.     if (sumY<64) {  
  124.         scale=1.4;  
  125.     }else if (sumY>192) {  
  126.         scale=0.6;  
  127.     }  
  128.   
  129.     if (abs(scale-1.)>0.001) {  
  130.         for (int h=0;h<height;h++) {  
  131.             double* psrc=(double*)(src->imageData+h*src->widthStep);  
  132.             unsigned char* pdst=(unsigned char*)dst->imageData+h*dst->widthStep;  
  133.             double t[3];  
  134.             for (int w=0;w<width;w++) {  
  135.                 t[0]=pow(psrc[R],scale);  
  136.                 t[1]=pow(psrc[G],scale);  
  137.                 t[2]=psrc[B];  
  138.                 for (int i=0;i<3;i++) {  
  139.                     if (t[i]>255) {  
  140.                         pdst[i]=255;  
  141.                     }else {  
  142.                         pdst[i]=(unsigned char)t[i];  
  143.                     }  
  144.                 }  
  145.                 psrc+=3;  
  146.                 pdst+=3;  
  147.             }  
  148.         }  
  149.     }else{  
  150.         for (int h=0;h<height;h++) {  
  151.             double* psrc=(double*)(src->imageData+h*src->widthStep);  
  152.             unsigned char* pdst=(unsigned char*)dst->imageData+h*dst->widthStep;  
  153.             for (int w=0;w<width;w++) {  
  154.                 for (int i=0;i<3;i++) {  
  155.                     double t=*psrc++;  
  156.                     if (t>255) {  
  157.                         *pdst++=255;  
  158.                     }else {  
  159.                         *pdst++=(unsigned char)t;  
  160.                     }                     
  161.                 }  
  162.             }  
  163.         }  
  164.     }  
  165.       
  166.     //free memory  
  167.     cvReleaseImage(&src);  
  168.     cvReleaseImage(&y);  
  169.     for (int i=0;i<3;i++) {  
  170.         cvReleaseImage(&rgb[i]);  
  171.     }  
  172. }  

matlab 版,算法没问题,但代码写的不太好,我做了点修改

[c-sharp]  view plain copy
  1. %%%%%%%%%%%%%%%RGB normalisation%%%%%%%%%%%%%%%%%%%%%%  
  2. %its cascaded implementain of 1 section of paper "A FAST SKIN REGION DETECTOR" by  
  3. %Phil Chen, Dr.Christos Greecos  
  4. %and  
  5. %section 2.1 of paper "Simple and accurate face detection in color images" by  
  6. %YUI TING PAI et al  
  7. % Coding by Madhava.S.Bhat, Dept of Electronics an communication  
  8. %Dr.Ambedkar Institute of Technology, Bangalore  
  9. %madhava.s@dr-ait.org  
  10. function[]= imag_improve_rgb(IMG)  
  11. figure,imshow(IMG)  
  12. title('original')  
  13. R=double(IMG(:,:,1));  
  14. G=double(IMG(:,:,2));  
  15. B=double(IMG(:,:,3));  
  16. [H,W]=size(R);  
  17.   
  18. % minR=0;  
  19. % minG=0;  
  20. % minB=0;  
  21. % [srow,scol]=find(R==0 & G==0 & B==0);  
  22. % if(isempty(srow) && isempty(scol))  
  23.     minR=min(min(R))  
  24.     minG=min(min(G))  
  25.     minB=min(min(B))  
  26. % end  
  27. R=R-minR;  
  28. G=G-minG;  
  29. B=B-minB;  
  30.   
  31. S=zeros(H,W);  
  32. [srow,scol]=find(R==0 & G==0 & B==0);  
  33. [sm,sn]=size(srow);  
  34.   
  35. for i=1:sm  
  36.      S(srow(i),scol(i))=1;  
  37. end  
  38. mstd=sum(sum(S))  
  39. Nstd=(H*W)-mstd;  
  40.   
  41. Cst=0;  
  42. Cst=double(Cst);  
  43. for i=1:H  
  44.     for j=1:W  
  45.          a=R(i,j);  
  46.          b=R(i,j);  
  47.            
  48.             if(B(i,j)<a)  
  49.                a=B(i,j);  
  50.             else  
  51.                b=B(i,j);  
  52.             end  
  53.            
  54.              if(G(i,j)<a)  
  55.                 a=G(i,j);  
  56.              else  
  57.                 b=G(i,j);  
  58.              end  
  59.            
  60.          Cst=a+b+Cst;  
  61.            
  62.     end  
  63. end  
  64. %%%%sum of black pixels%%%%%%%%%%%  
  65. blacksumR=0;  
  66. blacksumG=0;  
  67. blacksumB=0;  
  68. for i=1:sm  
  69.     blacksumR=blacksumR+R(srow(i),scol(i));  
  70.     blacksumG=blacksumG+G(srow(i),scol(i));  
  71.     blacksumB=blacksumB+B(srow(i),scol(i));  
  72. end  
  73. Cstd = Cst/(2*Nstd)  
  74. CavgR=sum(sum(R))./(H*W)  
  75. CavgB=sum(sum(B))./(H*W)  
  76. CavgG=sum(sum(G))./(H*W)  
  77. Rsc=Cstd./CavgR  
  78. Gsc=Cstd./CavgG  
  79. Bsc=Cstd/CavgB  
  80.   
  81. R=R.*Rsc;  
  82. G=G.*Gsc;  
  83. B=B.*Bsc;  
  84.   
  85. C(:,:,1)=R;  
  86. C(:,:,2)=G;  
  87. C(:,:,3)=B;  
  88. C=C/255;  
  89. YCbCr=rgb2ycbcr(C);  
  90. Y=YCbCr(:,:,1);  
  91.   
  92. figure,imshow(C)  
  93. title('aft 1st stage of compensation')  
  94.   
  95. %normalize Y  
  96. minY=min(min(Y));  
  97. maxY=max(max(Y));  
  98. Y=255.0*(Y-minY)./(maxY-minY);  
  99. YEye=Y;  
  100. Yavg=sum(sum(Y))/(W*H)  
  101.   
  102. T=1;  
  103. if (Yavg<64)  
  104.     T=1.4  
  105. elseif (Yavg>192)  
  106.     T=0.6  
  107. end  
  108. T  
  109.   
  110. if (T~=1)  
  111.     RI=R.^T;  
  112.     GI=G.^T;  
  113. else  
  114.     RI=R;  
  115.     GI=G;  
  116. end  
  117.   
  118.   
  119. Cfinal(:,:,1)=uint8(RI);  
  120. Cfinal(:,:,2)=uint8(GI);  
  121. Cfinal(:,:,3)=uint8(B);  
  122. figure,imshow(Cfinal)  
  123. title('Light intensity compensated')  
  124.   
  125. % YCbCr=rgb2ycbcr(Cnew);  

 转自:http://blog.csdn.net/onezeros/article/details/6342661

两组图像:左边较暗,右边较亮 第一行是原图像,他们下面是用四种算法处理的结果 依次为: 1.一种混合算法 2.msr,multi-scale retinex 3.msrcr,multi-scale retinex with color restoration 4

相关阅读排行


相关内容推荐

最新文章

×

×

请激活账号

为了能正常使用评论、编辑功能及以后陆续为用户提供的其他产品,请激活账号。

您的注册邮箱: 修改

重新发送激活邮件 进入我的邮箱

如果您没有收到激活邮件,请注意检查垃圾箱。