c++\OpenCV\多尺度Retinex图像增强算法

·retinex算法原理参考：
1、 https://www.cnblogs.com/sleepwalker/p/3676600.html
2、 https://blog.csdn.net/carson2005/article/details/9502053
retinex算法代码参考：
https://github.com/BBuf/Image-processing-algorithm/blob/master/Retinex%20MSRCR.cpp

	//opencv接口函数
	void process(const char * sFilename/*, const char * rFilename*/)
	{
IplImage * orig = NULL;
IplImage * dst = NULL;
unsigned char * sImage, *dImage;
int x, y, i;
int nWidth, nHeight, step;

orig = cvLoadImage(sFilename, 1); //打开图像   
if (orig == NULL)
{
	printf("Could not get image. Program exits!\n");
	exit(0);
}
nWidth = orig->width;
nHeight = orig->height;
step = orig->widthStep / sizeof(unsigned char);
dst = cvCreateImage(cvSize(nWidth, nHeight), IPL_DEPTH_8U, 3); // 创建目标图像   

sImage = (unsigned char*)malloc(sizeof(unsigned char)*(nHeight*nWidth * 3)); // 创建2个图像buffer   
dImage = (unsigned char*)malloc(sizeof(unsigned char)*(nHeight*nWidth * 3));

// 创建2个显示窗口，一个用于目标图像，一个用于源图像   
//cvNamedWindow( "Original Video", CV_WINDOW_AUTOSIZE );
cvNamedWindow( "Result Video", CV_WINDOW_AUTOSIZE );   
// 取图像进行处理   
//cvShowImage( "Original Video", orig );
//waitKey(0);
if (orig->nChannels == 3)
{
	for (y = 0; y < nHeight; y++)
		for (x = 0; x < nWidth; x++)
		{
			sImage[(y*nWidth + x)*orig->nChannels] = orig->imageData[y*step + x * orig->nChannels];
			sImage[(y*nWidth + x)*orig->nChannels + 1] = orig->imageData[y*step + x * orig->nChannels + 1];
			sImage[(y*nWidth + x)*orig->nChannels + 2] = orig->imageData[y*step + x * orig->nChannels + 2];
		}
}
memcpy(dImage, sImage, nWidth*nHeight*orig->nChannels);
/* 彩色图像增强 */
MSRCR(dImage, nWidth, nHeight, orig->nChannels);
/**/printf( "MSRCR parameters:\n" );
printf( "number of scales: \t\t%d\n", rvals.nscales );
printf( "each scale is respectively: \t\t" );
for( i = 0; i < rvals.nscales; i++)
{
	printf("[scale_%d] = %7.4f ",i+1,RetinexScales[i]);
}
printf( "\n" );

for (y = 0; y < nHeight; y++)
	for (x = 0; x < nWidth; x++)
	{
		dst->imageData[y*step + x * 3] = dImage[(y*nWidth + x) * 3];
		dst->imageData[y*step + x * 3 + 1] = dImage[(y*nWidth + x) * 3 + 1];
		dst->imageData[y*step + x * 3 + 2] = dImage[(y*nWidth + x) * 3 + 2];
	}
/************** 显示处理图像 ****************/
/**/cvShowImage( "Result Video", dst );
waitKey(0);

 //if (rFilename != NULL)
//cvSaveImage("C:\\Users\\h\\Desktop\\retinex\\result.jpg", dst, 0);  //并保存结果文件   
Mat matimg;
matimg = cvarrToMat(dst);
imwrite("C:\\Users\\h\\Desktop\\retinex\\04033_r.jpg", matimg);
cvWaitKey( 0 ); /**/ 
// 清除   
cvReleaseImage(&orig);
cvReleaseImage(&dst);
free(sImage);
free(dImage);
}

主函数

··
void main(int argc, char * argv[])
{
   
process("E:\\HRZBS_PRO\\Object-Detection-Using-Pretrained-YOLO-v2-main\\FLIR_03942.jpg"/*, "C:\\Users\\h\\Desktop\\retinex\\result.jpg"*/);
}

MSRCR

void MSRCR(unsigned char * src, int width, int height, int bytes)
	{
int           scale, row, col;
int           i, j;
int           size;
int           pos;
int           channel;
unsigned char *psrc = NULL;            /* backup pointer for src buffer */
float         *dst = NULL;            /* float buffer for algorithm */
float         *pdst = NULL;            /* backup pointer for float buffer */
float         *in, *out;
int           channelsize;            /* Float memory cache for one channel */
float         weight;
gauss3_coefs  coef;
float         mean, var;
float         mini, range, maxi;
	float         alpha;
	float         gain;
	float         offset;

#ifdef _DEBUG
IplImage * imL = NULL;
#endif

/* Allocate all the memory needed for algorithm*/
size = width * height * bytes;
dst = (float *)malloc(size * sizeof(float));
if (dst == NULL)
{
	printf("Failed to allocate memory");
	return;
}
memset(dst, 0, size * sizeof(float));

channelsize = (width * height);
in = (float *)malloc(channelsize * sizeof(float));
if (in == NULL)
{
	free(dst);
	printf("Failed to allocate memory");
	return; /* do some clever stuff */
}

out = (float *)malloc(channelsize * sizeof(float));
if (out == NULL)
{
	free(in);
	free(dst);
	printf("Failed to allocate memory");
	return; /* do some clever stuff */
}

/*
Calculate the scales of filtering according to the
number of filter and their distribution.
*/
retinex_scales_distribution(RetinexScales,
	rvals.nscales, rvals.scales_mode, rvals.scale);

/*
Filtering according to the various scales.
Summerize the results of the various filters according to a
specific weight(here equivalent for all).
*/
weight = 1.0f / (float)rvals.nscales;

/*
The recursive filtering algorithm needs different coefficients according
to the selected scale (~ = standard deviation of Gaussian).
*/
pos = 0;
for (channel = 0; channel < 3; channel++)
{
	for (i = 0, pos = channel; i < channelsize; i++, pos += bytes)
	{
		/* 0-255 => 1-256 */
		in[i] = (float)(src[pos] + 1.0);
	}
	for (scale = 0; scale < rvals.nscales; scale++)
	{
		compute_coefs3(&coef, RetinexScales[scale]);
		/*
		*  Filtering (smoothing) Gaussian recursive.
		*
		*  Filter rows first
		*/
		for (row = 0; row < height; row++)
		{
			pos = row * width;
			gausssmooth(in + pos, out + pos, width, 1, &coef);
		}

		memcpy(in, out, channelsize * sizeof(float));
		memset(out, 0, channelsize * sizeof(float));

		/*
		*  Filtering (smoothing) Gaussian recursive.
		*
		*  Second columns
		*/
		for (col = 0; col < width; col++)
		{
			pos = col;
			gausssmooth(in + pos, out + pos, height, width, &coef);
		}


		/*
		Summarize the filtered values.
		In fact one calculates a ratio between the original values and the filtered values.
		*/
		for (i = 0, pos = channel; i < channelsize; i++, pos += bytes)
		{
			dst[pos] += weight * (float)(log(src[pos] + 1.0f) - log(out[i]));
		}
	}
}
free(in);
free(out);

/*
Final calculation with original value and cumulated filter values.
The parameters gain, alpha and offset are constants.
*/
/* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */

alpha = 128.0f;
gain = 1.0f;
offset = 0.0f;

for (i = 0; i < size; i += bytes)
{
	float logl;

	psrc = src + i;
	pdst = dst + i;

	logl = (float)log((float)psrc[0] + (float)psrc[1] + (float)psrc[2] + 3.0f);

	pdst[0] = gain * ((float)(log(alpha * (psrc[0] + 1.0f)) - logl) * pdst[0]) + offset;
	pdst[1] = gain * ((float)(log(alpha * (psrc[1] + 1.0f)) - logl) * pdst[1]) + offset;
	pdst[2] = gain * ((float)(log(alpha * (psrc[2] + 1.0f)) - logl) * pdst[2]) + offset;
}

/*
Adapt the dynamics of the colors according to the statistics of the first and second order.
The use of the variance makes it possible to control the degree of saturation of the colors.
*/
pdst = dst;


//compute_mean_var_test(psrc, &mean, &var, size, bytes );  
compute_mean_var(pdst, &mean, &var, size, bytes);
mini = mean - rvals.cvar*var;
maxi = mean + rvals.cvar*var;
range = maxi - mini;

/*
printf( "variance: \t\t%7.4f\n", var * rvals.cvar );
printf( "mean: \t\t%7.4f\n", mean );
printf( "min: \t\t%7.4f\n", mini );
printf( "max: \t\t%7.4f\n", maxi );
printf( "range: \t\t%7.4f\n", range );
*/

if (!range) range = 1.0;

for (i = 0; i < size; i += bytes)
{
	psrc = src + i;
	pdst = dst + i;

	for (j = 0; j < 3; j++)
	{
		float c = 255 * (pdst[j] - mini) / range;

		psrc[j] = (unsigned char)clip(c, 0, 255);
	}
}

free(dst);

}

retinex图像增强算法结果展示