/*
* 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
- 1
- 2
- 3
- 4
- 5
- 6
前往页