视网膜皮层Retinex算法可以参考下面两篇文章
http://www.cnblogs.com/Imageshop/archive/2013/04/17/3026881.html
https://blog.csdn.net/carson2005/article/details/9502053
下面直接进行讲解基于Opencv 4.10的Retinex算法实现
/*主函数*/
int main()
{
Mat orig;
Mat dst;
unsigned char * sImage, *dImage;
int x, y;
int nWidth, nHeight, step;
orig = imread("test1.jpg", 1); /* 打开图像 */
nWidth = orig.cols;
nHeight = orig.rows;
step = orig.cols * orig.channels();
dst.create(Size(nWidth, nHeight), CV_8UC3);/* 创建目标图像 */
sImage = (unsigned char*)malloc(sizeof(unsigned char)*(nHeight*nWidth * 3)); /* 创建2个图像buffer */
dImage = (unsigned char*)malloc(sizeof(unsigned char)*(nHeight*nWidth * 3));
/* 创建2个显示窗口，一个用于目标图像，一个用于源图像 */
namedWindow("Original Image", WINDOW_AUTOSIZE);
namedWindow("Result Image", WINDOW_AUTOSIZE);
/* 取图像进行处理 */
imshow("Original Image", orig);
for (y = 0; y < orig.rows; y++)
{
uchar* data = orig.ptr<uchar>(y);
for (x = 0; x < orig.cols; x++)
{
sImage[(y*nWidth + x) * 3 + 0] = data[3 * x + 0];
sImage[(y*nWidth + x) * 3 + 1] = data[3 * x + 1];
sImage[(y*nWidth + x) * 3 + 2] = data[3 * x + 2];
}
}
/* 彩色图像增强 */
MSRCR(sImage, nWidth, nHeight, orig.channels());
for (y = 0; y < nHeight; y++)
{
uchar* data = dst.ptr<uchar>(y);
for (x = 0; x < nWidth; x++)
{
data[3 * x + 0] = sImage[(y*nWidth + x) * 3 + 0];
data[3 * x + 1] = sImage[(y*nWidth + x) * 3 + 1];
data[3 * x + 2] = sImage[(y*nWidth + x) * 3 + 2];
}
}
//显示处理图像
imshow("Result Image", dst);
waitKey(0);
free(sImage);
free(dImage);
return 0;
}下面的MSRCR是整个算法的和核心部分
/* 这个函数是算法的核心 */
/* (a) 在多个刻度处过滤，并将结果汇总 */
/* (b) 计算最终结果 */
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; /* SRC图的备份指针 */
float *dst = NULL; /* 算法的浮动缓冲区 */
float *pdst = NULL; /* 浮点缓冲区的备份指针 */
float *in, *out;
int channelsize; /* 一个通道的浮动内存缓存 */
float weight;
gauss3_coefs coef;
float mean, var;
float mini, range, maxi;
float alpha;
float gain;
float offset;
/* 分配算法所需的所有内存 */
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; /* 判断输入*/
}
out = (float *)malloc(channelsize * sizeof(float));
if (out == NULL)
{
free(in);
free(dst);
printf("Failed to allocate memory");
return; /* 判断输出 */
}
/* 根据过滤器数量及其分布计算过滤尺度 */
retinex_scales_distribution(RetinexScales, rvals.nscales, rvals.scales_mode, rvals.scale);
/* 根据不同的尺度数进行过滤 */
/* 根据一个特定的重量总结各种过滤器的结果（这里等同于所有过滤器）*/
weight = 1.0f / (float)rvals.nscales;
/* 根据所选的尺度，递归滤波算法需要不同的系数（~=高斯标准偏差）*/
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]);
/* 过滤（平滑）高斯递归 */
/* 先筛选行 */
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));
/* 过滤（平滑）高斯递归 */
/* 选择列 */
for (col = 0; col < width; col++)
{
pos = col;
gausssmooth(in + pos, out + pos, height, width, &coef);
}
/* 汇总过滤的值 */
//对每个尺度下进行累加计算 Log[R(x,y)] = Log[R(x,y)] + Weight(i)* ( Log[Ii(x,y)]-Log[Li(x,y)]);
//其中Weight(i)表示每个尺度对应的权重，要求各尺度权重之和必须为1，经典的取值为等权重
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);
/* 最终计算原始值和累积滤波器值 */
/* 参数增益、alpha和offset是常量 */
/* 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;
}
/* 根据一阶和二阶的统计数据调整颜色的动态 */
/* 使用方差可以控制颜色的饱和度 */
pdst = dst;
//计算RGB各通道数据的均值Mean和均方差Var
compute_mean_var(pdst, &mean, &var, size, bytes);
mini = mean - rvals.cvar*var;
maxi = mean + rvals.cvar*var;
range = maxi - mini;
if (!range) range = 1.0;
//对Log[R(x,y)]的每一个值，进行线性映射
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);
}计算期望分布的尺度
void retinex_scales_distribution(float* scales, int nscales, int mode, int s)
{
if (nscales == 1)
{
/* 一个参数选择中等尺度 */
scales[0] = (float)s / 2;
}
else if (nscales == 2)
{
/* 两个参数选择中和大尺度*/
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;
}
}
}计算一次平均值和方差
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 */
}高斯模糊的快速计算
关于高斯模糊如果有不清楚的可以参考这篇文章
https://www.cnblogs.com/invisible2/p/9177018.html
void compute_coefs3(gauss3_coefs * c, float sigma)
{
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)
{
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);
}输入的原图
经过retinex算法处理的图，可以看到原图几乎什么都看不清，但是retinex算法可以将图片的颜色较好的还原出来。
代码已经传到了csdn上，下载链接是
https://download.csdn.net/download/weixin_42521239/11165325
没有积分的小伙伴可以在评论中留下你的邮箱，我发给你
扫一扫
82
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
