图像增强算法四种，图示与源码，包括retinex（ssr、msr、msrcr）和一种混合算法

来源：互联网发布：java list sort 排序编辑：程序博客网时间：2024/05/16 23:44

申明：本文非笔者原创，原文转载自：http://blog.csdn.net/onezeros/article/details/6342661

两组图像：左边较暗，右边较亮

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

依次为：

1.一种混合算法

2.msr，multi-scale retinex

3.msrcr，multi-scale retinex with color restoration

4.ssr，single scale retinex

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

头文件：

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

实现：

[cpp] view plaincopy
/* 
 * Copyright (c) 2006, Douglas Gray (dgray@soe.ucsc.edu, dr.de3ug@gmail.com) 
 * All rights reserved. 
 * 
 * Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions are met: 
 *     * Redistributions of source code must retain the above copyright 
 *       notice, this list of conditions and the following disclaimer. 
 *     * Redistributions in binary form must reproduce the above copyright 
 *       notice, this list of conditions and the following disclaimer in the 
 *       documentation and/or other materials provided with the distribution. 
 *     * Neither the name of the <organization> nor the 
 *       names of its contributors may be used to endorse or promote products 
 *       derived from this software without specific prior written permission. 
 * 
 * THIS SOFTWARE IS PROVIDED BY Douglas Gray ``AS IS'' AND ANY 
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL <copyright holder> BE LIABLE FOR ANY 
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
 */  
#include "retinex.h"  
  
#include <math.h>  
  
//#define USE_EXACT_SIGMA  
  
  
#define pc(image, x, y, c) image->imageData[(image->widthStep * y) + (image->nChannels * x) + c]  
  
#define INT_PREC 1024.0  
#define INT_PREC_BITS 10  
  
inline double int2double(int x) { return (double)x / INT_PREC; }  
inline int double2int(double x) { return (int)(x * INT_PREC + 0.5); }  
  
inline int int2smallint(int x) { return (x >> INT_PREC_BITS); }  
inline int int2bigint(int x) { return (x << INT_PREC_BITS); }  
  
//  
// CreateKernel  
//  
// Summary:  
// Creates a normalized 1 dimensional gaussian kernel.  
//  
// Arguments:  
// sigma - the standard deviation of the gaussian kernel.  
//  
// Returns:  
// double* - an array of values of length ((6*sigma)/2) * 2 + 1.  
//  
// Note:  
// Caller is responsable for deleting the kernel.  
//  
double*  
CreateKernel(double sigma)  
{  
    int i, x, filter_size;  
    double* filter;  
    double sum;  
  
    // Reject unreasonable demands  
    if ( sigma > 200 ) sigma = 200;  
  
    // get needed filter size (enforce oddness)  
    filter_size = (int)floor(sigma*6) / 2;  
    filter_size = filter_size * 2 + 1;  
  
    // Allocate kernel space  
    filter = new double[filter_size];  
  
    // Calculate exponential  
    sum = 0;  
    for (i = 0; i < filter_size; i++) {  
        x = i - (filter_size / 2);  
        filter[i] = exp( -(x*x) / (2*sigma*sigma) );  
  
        sum += filter[i];  
    }  
  
    // Normalize  
    for (i = 0, x; i < filter_size; i++)  
        filter[i] /= sum;  
  
    return filter;  
}  
  
//  
// CreateFastKernel  
//  
// Summary:  
// Creates a faster gaussian kernal using integers that  
// approximate floating point (leftshifted by 8 bits)  
//  
// Arguments:  
// sigma - the standard deviation of the gaussian kernel.  
//  
// Returns:  
// int* - an array of values of length ((6*sigma)/2) * 2 + 1.  
//  
// Note:  
// Caller is responsable for deleting the kernel.  
//  
  
int*  
CreateFastKernel(double sigma)  
{  
    double* fp_kernel;  
    int* kernel;  
    int i, filter_size;  
      
    // Reject unreasonable demands  
    if ( sigma > 200 ) sigma = 200;  
  
    // get needed filter size (enforce oddness)  
    filter_size = (int)floor(sigma*6) / 2;  
    filter_size = filter_size * 2 + 1;  
  
    // Allocate kernel space  
    kernel = new int[filter_size];  
  
    fp_kernel = CreateKernel(sigma);  
  
    for (i = 0; i < filter_size; i++)  
        kernel[i] = double2int(fp_kernel[i]);  
  
    delete fp_kernel;  
  
    return kernel;  
}  
  
  
//  
// FilterGaussian  
//  
// Summary:  
// Performs a gaussian convolution for a value of sigma that is equal  
// in both directions.  
//  
// Arguments:  
// img - the image to be filtered in place.  
// sigma - the standard deviation of the gaussian kernel to use.  
//  
void  
FilterGaussian(IplImage* img, double sigma)  
{  
    int i, j, k, source, filter_size;  
    int* kernel;  
    IplImage* temp;  
    int v1, v2, v3;  
  
    // Reject unreasonable demands  
    if ( sigma > 200 ) sigma = 200;  
  
    // get needed filter size (enforce oddness)  
    filter_size = (int)floor(sigma*6) / 2;  
    filter_size = filter_size * 2 + 1;  
  
    kernel = CreateFastKernel(sigma);  
  
    temp = cvCreateImage(cvSize(img->width, img->height), img->depth, img->nChannels);  
      
    // filter x axis  
    for (j = 0; j < temp->height; j++)  
    for (i = 0; i < temp->width; i++) {  
  
        // inner loop has been unrolled  
  
        v1 = v2 = v3 = 0;  
        for (k = 0; k < filter_size; k++) {  
              
            source = i + filter_size / 2 - k;  
      
            if (source < 0) source *= -1;  
            if (source > img->width - 1) source = 2*(img->width - 1) - source;  
  
            v1 += kernel[k] * (unsigned char)pc(img, source, j, 0);  
            if (img->nChannels == 1) continue;  
            v2 += kernel[k] * (unsigned char)pc(img, source, j, 1);  
            v3 += kernel[k] * (unsigned char)pc(img, source, j, 2);  
  
        }  
  
        // set value and move on  
        pc(temp, i, j, 0) = (char)int2smallint(v1);  
        if (img->nChannels == 1) continue;  
        pc(temp, i, j, 1) = (char)int2smallint(v2);  
        pc(temp, i, j, 2) = (char)int2smallint(v3);  
  
    }  
      
    // filter y axis  
    for (j = 0; j < img->height; j++)  
    for (i = 0; i < img->width; i++) {  
  
        v1 = v2 = v3 = 0;  
        for (k = 0; k < filter_size; k++) {  
              
            source = j + filter_size / 2 - k;  
      
            if (source < 0) source *= -1;  
            if (source > temp->height - 1) source = 2*(temp->height - 1) - source;  
  
            v1 += kernel[k] * (unsigned char)pc(temp, i, source, 0);  
            if (img->nChannels == 1) continue;  
            v2 += kernel[k] * (unsigned char)pc(temp, i, source, 1);  
            v3 += kernel[k] * (unsigned char)pc(temp, i, source, 2);  
  
        }  
  
        // set value and move on  
        pc(img, i, j, 0) = (char)int2smallint(v1);  
        if (img->nChannels == 1) continue;  
        pc(img, i, j, 1) = (char)int2smallint(v2);  
        pc(img, i, j, 2) = (char)int2smallint(v3);  
  
    }  
  
  
    cvReleaseImage( &temp );  
  
    delete kernel;  
  
}  
  
//  
// FastFilter  
//  
// Summary:  
// Performs gaussian convolution of any size sigma very fast by using  
// both image pyramids and seperable filters.  Recursion is used.  
//  
// Arguments:  
// img - an IplImage to be filtered in place.  
//  
void  
FastFilter(IplImage *img, double sigma)  
{  
    int filter_size;  
  
    // Reject unreasonable demands  
    if ( sigma > 200 ) sigma = 200;  
  
    // get needed filter size (enforce oddness)  
    filter_size = (int)floor(sigma*6) / 2;  
    filter_size = filter_size * 2 + 1;  
  
    // If 3 sigma is less than a pixel, why bother (ie sigma < 2/3)  
    if(filter_size < 3) return;  
  
    // Filter, or downsample and recurse  
    if (filter_size < 10) {  
  
#ifdef USE_EXACT_SIGMA  
        FilterGaussian(img, sigma)  
#else  
        cvSmooth( img, img, CV_GAUSSIAN, filter_size, filter_size );  
#endif  
              
    }  
    else {  
        if (img->width < 2 || img->height < 2) return;  
  
        IplImage* sub_img = cvCreateImage(cvSize(img->width / 2, img->height / 2), img->depth, img->nChannels);  
  
        cvPyrDown( img, sub_img );  
  
        FastFilter( sub_img, sigma / 2.0 );  
  
        cvResize( sub_img, img, CV_INTER_LINEAR );  
  
        cvReleaseImage( &sub_img );  
    }  
  
}  
  
//  
// Retinex  
//  
// Summary:  
// Basic retinex restoration.  The image and a filtered image are converted  
// to the log domain and subtracted.  
//  
// Arguments:  
// img - an IplImage to be enhanced in place.  
// sigma - the standard deviation of the gaussian kernal used to filter.  
// gain - the factor by which to scale the image back into visable range.  
// offset - an offset similar to the gain.  
//  
void  
Retinex(IplImage *img, double sigma, int gain, int offset)  
{  
    IplImage *A, *fA, *fB, *fC;  
  
    // Initialize temp images  
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  
    // Compute log image  
    cvConvert( img, fA );  
    cvLog( fA, fB );  
      
    // Compute log of blured image  
    A = cvCloneImage( img );  
    FastFilter( A, sigma );  
    cvConvert( A, fA );  
    cvLog( fA, fC );  
  
    // Compute difference  
    cvSub( fB, fC, fA );  
  
    // Restore  
    cvConvertScale( fA, img, gain, offset);  
  
    // Release temp images  
    cvReleaseImage( &A );  
    cvReleaseImage( &fA );  
    cvReleaseImage( &fB );  
    cvReleaseImage( &fC );  
  
}  
  
//  
// MultiScaleRetinex  
//  
// Summary:  
// Multiscale retinex restoration.  The image and a set of filtered images are  
// converted to the log domain and subtracted from the original with some set  
// of weights. Typicaly called with three equaly weighted scales of fine,  
// medium and wide standard deviations.  
//  
// Arguments:  
// img - an IplImage to be enhanced in place.  
// sigma - the standard deviation of the gaussian kernal used to filter.  
// gain - the factor by which to scale the image back into visable range.  
// offset - an offset similar to the gain.  
//  
void  
MultiScaleRetinex(IplImage *img, int scales, double *weights, double *sigmas, int gain, int offset)  
{  
    int i;  
    double weight;  
    IplImage *A, *fA, *fB, *fC;  
  
    // Initialize temp images  
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
  
  
    // Compute log image  
    cvConvert( img, fA );  
    cvLog( fA, fB );  
  
    // Normalize according to given weights  
    for (i = 0, weight = 0; i < scales; i++)  
        weight += weights[i];  
  
    if (weight != 1.0) cvScale( fB, fB, weight );  
  
    // Filter at each scale  
    for (i = 0; i < scales; i++) {  
        A = cvCloneImage( img );  
        FastFilter( A, sigmas[i] );  
      
        cvConvert( A, fA );  
        cvLog( fA, fC );  
        cvReleaseImage( &A );  
  
        // Compute weighted difference  
        cvScale( fC, fC, weights[i] );  
        cvSub( fB, fC, fB );  
    }  
  
    // Restore  
    cvConvertScale( fB, img, gain, offset);  
  
    // Release temp images  
    cvReleaseImage( &fA );  
    cvReleaseImage( &fB );  
    cvReleaseImage( &fC );  
}  
  
//  
// MultiScaleRetinexCR  
//  
// Summary:  
// Multiscale retinex restoration with color restoration.  The image and a set of  
// filtered images are converted to the log domain and subtracted from the  
// original with some set of weights. Typicaly called with three equaly weighted  
// scales of fine, medium and wide standard deviations. A color restoration weight  
// is then applied to each color channel.  
//  
// Arguments:  
// img - an IplImage to be enhanced in place.  
// sigma - the standard deviation of the gaussian kernal used to filter.  
// gain - the factor by which to scale the image back into visable range.  
// offset - an offset similar to the gain.  
// restoration_factor - controls the non-linearaty of the color restoration.  
// color_gain - controls the color restoration gain.  
//  
void  
MultiScaleRetinexCR(IplImage *img, int scales, double *weights, double *sigmas,  
                    int gain, int offset, double restoration_factor, double color_gain)  
{  
    int i;  
    double weight;  
    IplImage *A, *B, *C, *fA, *fB, *fC, *fsA, *fsB, *fsC, *fsD, *fsE, *fsF;  
  
    // Initialize temp images  
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);  
    fsA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
    fsB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
    fsC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
    fsD = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
    fsE = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
    fsF = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);  
  
    // Compute log image  
    cvConvert( img, fB );  
    cvLog( fB, fA );  
  
    // Normalize according to given weights  
    for (i = 0, weight = 0; i < scales; i++)  
        weight += weights[i];  
  
    if (weight != 1.0) cvScale( fA, fA, weight );  
  
    // Filter at each scale  
    for (i = 0; i < scales; i++) {  
        A = cvCloneImage( img );  
        FastFilter( A, sigmas[i] );  
      
        cvConvert( A, fB );  
        cvLog( fB, fC );  
        cvReleaseImage( &A );  
  
        // Compute weighted difference  
        cvScale( fC, fC, weights[i] );  
        cvSub( fA, fC, fA );  
    }  
  
    // Color restoration  
    if (img->nChannels > 1) {  
        A = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
        B = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
        C = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);  
          
        // Divide image into channels, convert and store sum  
        cvCvtPixToPlane( img, A, B, C, NULL );  
  
        cvConvert( A, fsA );  
        cvConvert( B, fsB );  
        cvConvert( C, fsC );  
  
        cvReleaseImage( &A );  
        cvReleaseImage( &B );  
        cvReleaseImage( &C );  
  
        // Sum components  
        cvAdd( fsA, fsB, fsD );  
        cvAdd( fsD, fsC, fsD );  
  
        // Normalize weights  
        cvDiv( fsA, fsD, fsA, restoration_factor);  
        cvDiv( fsB, fsD, fsB, restoration_factor);  
        cvDiv( fsC, fsD, fsC, restoration_factor);  
  
        cvConvertScale( fsA, fsA, 1, 1 );  
        cvConvertScale( fsB, fsB, 1, 1 );  
        cvConvertScale( fsC, fsC, 1, 1 );  
  
        // Log weights  
        cvLog( fsA, fsA );  
        cvLog( fsB, fsB );  
        cvLog( fsC, fsC );  
  
        // Divide retinex image, weight accordingly and recombine  
        cvCvtPixToPlane( fA, fsD, fsE, fsF, NULL );  
  
        cvMul( fsD, fsA, fsD, color_gain);  
        cvMul( fsE, fsB, fsE, color_gain );  
        cvMul( fsF, fsC, fsF, color_gain );  
  
        cvCvtPlaneToPix( fsD, fsE, fsF, NULL, fA );  
    }  
  
    // Restore  
    cvConvertScale( fA, img, gain, offset);  
  
    // Release temp images  
    cvReleaseImage( &fA );  
    cvReleaseImage( &fB );  
    cvReleaseImage( &fC );  
    cvReleaseImage( &fsA );  
    cvReleaseImage( &fsB );  
    cvReleaseImage( &fsC );  
    cvReleaseImage( &fsD );  
    cvReleaseImage( &fsE );  
    cvReleaseImage( &fsF );  
}  

这种混合算法代码：

我参考一个matlab代码写的

c版

[cpp] view plaincopy
// author : onezeros.lee@gmail.com  
// data   : 4/20/2011   
/* 
%%%%%%%%%%%%%%%RGB normalization%%%%%%%%%%%%%%%%%%%%%% 
%its cascaded implementation of 1 section of paper "A FAST SKIN REGION DETECTOR" by 
%Phil Chen, Dr.Christos Greecos 
%and 
%section 2.1 of paper "Simple and accurate face detection in color images" by 
%YUI TING PAI et al 
% Coding by Madhava.S.Bhat, Dept of Electronics an communication 
%Dr.Ambedkar Institute of Technology, Bangalore 
%madhava.s@dr-ait.org 
 */  
  
#include <cv.h>  
#include <cxcore.h>  
#include <highgui.h>  
  
#include <iostream>  
#include <algorithm>  
using namespace std;  
  
void cvCompensationGlobal1(IplImage* _src,IplImage* dst)  
{  
    static const int B=0;  
    static const int G=1;  
    static const int R=2;  
  
    int width=_src->width;  
    int height=_src->height;  
    //double  
    IplImage* src=cvCreateImage(cvGetSize(_src),IPL_DEPTH_64F,3);  
    CvScalar minV=cvScalar(255,255,255,0);  
    for (int h=0;h<height;h++) {  
        unsigned char* p=(unsigned char*)_src->imageData+h*_src->widthStep;  
        double* pc=(double*)(src->imageData+h*src->widthStep);  
        for (int w=0;w<width;w++) {  
            for (int i=0;i<3;i++) {  
                *pc=*p;  
                if (minV.val[i]>*pc) {  
                    minV.val[i]=*pc;  
                }  
                pc++;  
                p++;  
            }  
        }  
    }  
    cvSubS(src,minV,src);  
  
    int blackNum=0;  
    double total=0;  
    double acc[3]={0};  
    for (int h=0;h<height;h++) {  
        double* p=(double*)(src->imageData+h*src->widthStep);  
        for (int w=0;w<width;w++) {  
            if (p[0]<0.001&&p[1]<0.001&&p[2]<0.001) {  
                blackNum++;  
                p+=3;  
                continue;  
            }  
            double a=p[R];  
            double b=p[R];  
            if (p[B]<a) {  
                a=p[B];  
            }else {  
                b=p[B];  
            }  
            if (p[G]<a) {  
                a=p[G];  
            }else {  
                b=p[G];  
            }  
            total+=a+b;  
  
            for (int i=0;i<3;i++) {  
                acc[i]+=p[i];  
            }  
            p+=3;  
        }  
    }  
    double avgT=total/(2*(width*height-blackNum));  
    CvScalar avg;  
    IplImage* rgb[3];  
    for (int i=0;i<3;i++) {  
        rgb[i]=cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
        avg.val[i]=(double)acc[i]/(width*height);  
        avg.val[i]=avgT/avg.val[i];  
    }  
    cvSplit(src,rgb[0],rgb[1],rgb[2],0);  
    for (int i=0;i<3;i++) {  
        cvScale(rgb[i],rgb[i],avg.val[i]);//bigger than 255  
    }  
    cvMerge(rgb[0],rgb[1],rgb[2],0,src);  
  
    //y component only  
    IplImage* y=cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);  
    for (int h=0;h<height;h++) {  
        double* psrc=(double*)(src->imageData+h*src->widthStep);  
        double* py=(double*)(y->imageData+h*y->widthStep);  
        for (int w=0;w<width;w++) {  
            *py++=psrc[R]*0.299+psrc[G]*0.587+psrc[B]*0.114;  
            psrc+=3;  
        }  
    }  
    double maxY=0;  
    double minY=0;  
    cvMinMaxLoc(y,&minY,&maxY);  
      
    double sumY=0;  
    //scale   
    for (int h=0;h<height;h++) {  
        double* p=(double*)(y->imageData+h*y->widthStep);  
        for (int w=0;w<width;w++) {  
            *p=(*p-minY)/(maxY-minY);  
            sumY+=*p;  
            p++;  
        }  
    }  
    sumY*=255;  
    sumY/=width*height;  
      
    double scale=1.;  
    if (sumY<64) {  
        scale=1.4;  
    }else if (sumY>192) {  
        scale=0.6;  
    }  
  
    if (abs(scale-1.)>0.001) {  
        for (int h=0;h<height;h++) {  
            double* psrc=(double*)(src->imageData+h*src->widthStep);  
            unsigned char* pdst=(unsigned char*)dst->imageData+h*dst->widthStep;  
            double t[3];  
            for (int w=0;w<width;w++) {  
                t[0]=pow(psrc[R],scale);  
                t[1]=pow(psrc[G],scale);  
                t[2]=psrc[B];  
                for (int i=0;i<3;i++) {  
                    if (t[i]>255) {  
                        pdst[i]=255;  
                    }else {  
                        pdst[i]=(unsigned char)t[i];  
                    }  
                }  
                psrc+=3;  
                pdst+=3;  
            }  
        }  
    }else{  
        for (int h=0;h<height;h++) {  
            double* psrc=(double*)(src->imageData+h*src->widthStep);  
            unsigned char* pdst=(unsigned char*)dst->imageData+h*dst->widthStep;  
            for (int w=0;w<width;w++) {  
                for (int i=0;i<3;i++) {  
                    double t=*psrc++;  
                    if (t>255) {  
                        *pdst++=255;  
                    }else {  
                        *pdst++=(unsigned char)t;  
                    }                     
                }  
            }  
        }  
    }  
      
    //free memory  
    cvReleaseImage(&src);  
    cvReleaseImage(&y);  
    for (int i=0;i<3;i++) {  
        cvReleaseImage(&rgb[i]);  
    }  
}  

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

[c-sharp] view plaincopy
%%%%%%%%%%%%%%%RGB normalisation%%%%%%%%%%%%%%%%%%%%%%  
%its cascaded implementain of 1 section of paper "A FAST SKIN REGION DETECTOR" by  
%Phil Chen, Dr.Christos Greecos  
%and  
%section 2.1 of paper "Simple and accurate face detection in color images" by  
%YUI TING PAI et al  
% Coding by Madhava.S.Bhat, Dept of Electronics an communication  
%Dr.Ambedkar Institute of Technology, Bangalore  
%madhava.s@dr-ait.org  
function[]= imag_improve_rgb(IMG)  
figure,imshow(IMG)  
title('original')  
R=double(IMG(:,:,1));  
G=double(IMG(:,:,2));  
B=double(IMG(:,:,3));  
[H,W]=size(R);  
  
% minR=0;  
% minG=0;  
% minB=0;  
% [srow,scol]=find(R==0 & G==0 & B==0);  
% if(isempty(srow) && isempty(scol))  
    minR=min(min(R))  
    minG=min(min(G))  
    minB=min(min(B))  
% end  
R=R-minR;  
G=G-minG;  
B=B-minB;  
  
S=zeros(H,W);  
[srow,scol]=find(R==0 & G==0 & B==0);  
[sm,sn]=size(srow);  
  
for i=1:sm  
     S(srow(i),scol(i))=1;  
end  
mstd=sum(sum(S))  
Nstd=(H*W)-mstd;  
  
Cst=0;  
Cst=double(Cst);  
for i=1:H  
    for j=1:W  
         a=R(i,j);  
         b=R(i,j);  
           
            if(B(i,j)<a)  
               a=B(i,j);  
            else  
               b=B(i,j);  
            end  
           
             if(G(i,j)<a)  
                a=G(i,j);  
             else  
                b=G(i,j);  
             end  
           
         Cst=a+b+Cst;  
           
    end  
end  
%%%%sum of black pixels%%%%%%%%%%%  
blacksumR=0;  
blacksumG=0;  
blacksumB=0;  
for i=1:sm  
    blacksumR=blacksumR+R(srow(i),scol(i));  
    blacksumG=blacksumG+G(srow(i),scol(i));  
    blacksumB=blacksumB+B(srow(i),scol(i));  
end  
Cstd = Cst/(2*Nstd)  
CavgR=sum(sum(R))./(H*W)  
CavgB=sum(sum(B))./(H*W)  
CavgG=sum(sum(G))./(H*W)  
Rsc=Cstd./CavgR  
Gsc=Cstd./CavgG  
Bsc=Cstd/CavgB  
  
R=R.*Rsc;  
G=G.*Gsc;  
B=B.*Bsc;  
  
C(:,:,1)=R;  
C(:,:,2)=G;  
C(:,:,3)=B;  
C=C/255;  
YCbCr=rgb2ycbcr(C);  
Y=YCbCr(:,:,1);  
  
figure,imshow(C)  
title('aft 1st stage of compensation')  
  
%normalize Y  
minY=min(min(Y));  
maxY=max(max(Y));  
Y=255.0*(Y-minY)./(maxY-minY);  
YEye=Y;  
Yavg=sum(sum(Y))/(W*H)  
  
T=1;  
if (Yavg<64)  
    T=1.4  
elseif (Yavg>192)  
    T=0.6  
end  
T  
  
if (T~=1)  
    RI=R.^T;  
    GI=G.^T;  
else  
    RI=R;  
    GI=G;  
end  
  
  
Cfinal(:,:,1)=uint8(RI);  
Cfinal(:,:,2)=uint8(GI);  
Cfinal(:,:,3)=uint8(B);  
figure,imshow(Cfinal)  
title('Light intensity compensated')  
  
% YCbCr=rgb2ycbcr(Cnew);  