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

最新推荐文章于 2024-08-17 11:30:52 发布

onezeros

最新推荐文章于 2024-08-17 11:30:52 发布

阅读量4.6w

点赞数 7

分类专栏： computer vision 文章标签：算法 filter image documentation matlab express

本文链接：https://blog.csdn.net/onezeros/article/details/6342661

版权

computer vision 专栏收录该内容

25 篇文章 2 订阅

订阅专栏

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

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

依次为：

1.一种混合算法

2.msr，multi-scale retinex

3.msrcr，multi-scale retinex with color restoration

4.ssr，single scale retinex

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

头文件：

实现：

/* * 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版

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

%%%%%%%%%%%%%%%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);