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

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

分类： CV相关2013-07-24
20:55 2145人阅读 评论(0) 收藏 举报

申明：本文非笔者原创，原文转载自：/article/7910142.html

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

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

依次为：

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<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 版，算法没问题，但代码写的不太好，我做了点修改

[b][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);