色彩恢复的多尺度Retinex算法（MSRCR）-VS2010版

/* * MSRCR：带彩色恢复的多尺度Retinex图像增强 * (Multi-Scale Retinex with Color Restoration) * 改写自：2003 Fabien Pelisson <Fabien.Pelisson@inrialpes.fr>的 * GRetinex GIMP plug-in * * Copyright (C) 2009 MAO Y.B * 2009. 3. 3 * Visual Information Processing (VIP) Group, NJUST * * 算法细节请参考下面论文： * D. J. Jobson, Z. Rahman, and G. A. Woodell. A multi-scale * Retinex for bridging the gap between color images and the * human observation of scenes. IEEE Transactions on Image Processing, * 1997, 6(7): 965-976 * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * */ # include <stdlib.h> # include <stdio.h> # include <math.h> # include <string.h> # include <opencv\cxcore.h> # include <opencv\cv.h> # include <opencv\highgui.h> #ifdef _DEBUG #pragma comment(lib,"opencv_imgproc246d.lib") #pragma comment(lib,"opencv_highgui246d.lib") #pragma comment(lib,"opencv_core246d.lib") #else #pragma comment(lib,"opencv_imgproc246.lib") #pragma comment(lib,"opencv_highgui246.lib") #pragma comment(lib,"opencv_core246.lib") #endif # define MAX_RETINEX_SCALES 8 /* Retinex最多可采用的尺度的数目 */ # define MIN_GAUSSIAN_SCALE 16 /* 最小Gaussian尺度 */ # define MAX_GAUSSIAN_SCALE 250 /* 最大Gaussian尺度 */ typedef struct { int scale; /* 最大Retinex尺度 */ int nscales; /* 尺度个数 */ int scales_mode; /* Retinex尺度计算模式，有3种：UNIFORM, LOW, HIGH */ float cvar; /* 用于调整色彩动态范围的方差的倍乘系数 */ } RetinexParams; /* 3种Retinex尺度计算模式，均匀，低和高，它们决定RetinexScales中的尺度数据 */ # define RETINEX_UNIFORM 0 # define RETINEX_LOW 1 # define RETINEX_HIGH 2 /* 多尺度Retinex中需要的各个Retinex尺度保存在下面数组中 */ static float RetinexScales[MAX_RETINEX_SCALES]; typedef struct { int N; float sigma; double B; double b[4]; } gauss3_coefs; /* * Private variables. */ static RetinexParams rvals = { 240, /* Scale */ 3, /* Scales */ RETINEX_UNIFORM, /* Retinex processing mode */ 1.2f /* A variant */ }; # define clip( val, minv, maxv ) (( val = (val < minv ? minv : val ) ) > maxv ? maxv : val ) /* * calculate scale values for desired distribution. */ void retinex_scales_distribution( float* scales, int nscales, int mode, int s) { if ( nscales == 1 ) { /* For one filter we choose the median scale */ scales[0] = (float)s / 2; } else if (nscales == 2) { /* For two filters we choose the median and maximum scale */ scales[0] = (float) s / 2; scales[1] = (float) s; } else { float size_step = (float) s / (float) nscales; int i; switch( mode ) { case RETINEX_UNIFORM: for ( i = 0; i < nscales; ++i ) scales[i] = 2.0f + (float)i * size_step; break; case RETINEX_LOW: size_step = (float)log(s - 2.0f) / (float) nscales; for ( i = 0; i < nscales; ++i ) scales[i] = 2.0f + (float)pow (10, (i * size_step) / log (10)); break; case RETINEX_HIGH: size_step = (float) log(s - 2.0) / (float) nscales; for ( i = 0; i < nscales; ++i ) scales[i] = s - (float)pow (10, (i * size_step) / log (10)); break; default: break; } } } /* * Calculate the average and variance in one go. */ void compute_mean_var( float *src, float *mean, float *var, int size, int bytes ) { float vsquared; int i, j; float *psrc; vsquared = 0.0f; *mean = 0.0f; for ( i = 0; i < size; i+= bytes ) { psrc = src+i; for ( j = 0 ; j < 3 ; j++ ) { *mean += psrc[j]; vsquared += psrc[j] * psrc[j]; } } *mean /= (float) size; /* mean */ vsquared /= (float) size; /* mean (x^2) */ *var = ( vsquared - (*mean * *mean) ); *var = (float)sqrt(*var); /* var */ } void compute_mean_var_test( unsigned char*src, float *mean, float *var, int size, int bytes ) { float vsquared; int i, j; unsigned char *psrc; vsquared = 0.0f; *mean = 0.0f; for ( i = 0; i < size; i+= bytes ) { psrc = src+i; for ( j = 0 ; j < 3 ; j++ ) { *mean += psrc[j]; vsquared += psrc[j] * psrc[j]; } } *mean /= (float) size; /* mean */ vsquared /= (float) size; /* mean (x^2) */ *var = ( vsquared - (*mean * *mean) ); *var = (float)sqrt(*var); /* var */ } /* * Calculate the coefficients for the recursive filter algorithm * Fast Computation of gaussian blurring. */ void compute_coefs3( gauss3_coefs * c, float sigma ) { /* * Papers: "Recursive Implementation of the gaussian filter.", * Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995. * formula: 11b computation of q * 8c computation of b0..b1 * 10 alpha is normalization constant B */ float q, q2, q3; q = 0; if ( sigma >= 2.5f ) { q = 0.98711f * sigma - 0.96330f; } else if ( (sigma >= 0.5f) && (sigma < 2.5f) ) { q = 3.97156f - 4.14554f * (float) sqrt ((double) 1 - 0.26891 * sigma); } else { q = 0.1147705018520355224609375f; } q2 = q * q; q3 = q * q2; c->b[0] = (1.57825f+(2.44413f*q)+(1.4281f *q2)+(0.422205f*q3)); c->b[1] = ( (2.44413f*q)+(2.85619f*q2)+(1.26661f *q3)); c->b[2] = ( -((1.4281f*q2)+(1.26661f *q3))); c->b[3] = ( (0.422205f*q3)); c->B = 1.0f-((c->b[1]+c->b[2]+c->b[3])/c->b[0]); c->sigma = sigma; c->N = 3; } void gausssmooth( float *in, float *out, int size, int rowstride, gauss3_coefs *c ) { /* * Papers: "Recursive Implementation of the gaussian filter.", * Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995. * formula: 9a forward filter * 9b backward filter * fig7 algorithm */ int i,n, bufsize; float *w1,*w2; /* forward pass */ bufsize = size+3; size -= 1; w1 = (float *)malloc (bufsize * sizeof (float)); w2 = (float *)malloc (bufsize * sizeof (float)); w1[0] = in[0]; w1[1] = in[0]; w1[2] = in[0]; for ( i = 0 , n=3; i <= size ; i++, n++) { w1[n] = (float)(c->B*in[i*rowstride] + ((c->b[1]*w1[n-1] + c->b[2]*w1[n-2] + c->b[3]*w1[n-3] ) / c->b[0])); } /* backward pass */ w2[size+1]= w1[size+3]; w2[size+2]= w1[size+3]; w2[size+3]= w1[size+3]; for ( i = size, n = i; i >= 0; i--, n-- ) { w2[n]= out[i * rowstride] = (float)(c->B*w1[n] + ((c->b[1]*w2[n+1] + c->b[2]*w2[n+2] + c->b[3]*w2[n+3] ) / c->b[0])); } free (w1); free (w2); } /* * This function is the heart of the algo. * (a) Filterings at several scale