项目地址：
照度反射模型
入射光线 (L(x,y)) 经过物体 (R(x,y)) 反射后，被相机捕获，形成图像 (S(x,y)) 。即
单尺度 Retinex
原理
单尺度 Retinex 即 Single Scale Retinex，它是所有基于 Retinex 方法的最简版本，原理如下：
上式两边取对数变换：
于是得到物体真实图像 (R)
但是入射光线 (L) 真实情况下是无法获取的，于是通过高斯函数进行卷积，获取高斯模糊图像，来近似地估算光照分量：
于是，对于照度反射模型，可整理出一个通用的表达式：
即给定一个输入图像 (S) 和一个高斯函数 (G) 就可以根据上述公式计算出物体真实图像 (R) ，从而一定程度上对将拍摄得到的图像中光照带来的影响减弱。这里的高斯函数表达式如下：
且要满足归一化条件
在实际操作中，卷积的实现有多种方法，其中比较快速的做法是转换到频域做乘积，因为
时域中（图像中通常称为空域）的卷积等于频域中的乘积
而且通常得到的最后的 (log{R}) 在对数域中，尽管做指数变换可以还原得到 (R)，但图像已经减去了反应光照情况的分量，所以得到的图像视觉效果并不好，为此需要做增益补偿，通常可直接对 (log{R}) 做一个线性拉伸处理。
实现
function [R,R_temp,L] = ssr4gray(S,sigma)
% S = 0 ~ 255 , 单通道图像
% R = 0 ~ 1 , 单通道图像
% R_temp = any , 用于进一步计算，避免做归一化影响结果
% S = R×L
% log(R) = log(S) - log(L)
% = log(S) - log(S*L)
G = fspecial('gaussian',size(S),sigma); % 高斯函数
% 时域卷积等于频域乘积，估计光照图像 L
%
L = real(ifft2(fft2(S).*fftshift(fft2(G,size(S,1),size(S,2)))));
L = L - min(L(:));
%}
% L = imfilter(S,G,'replicate','conv'); % 使用滤波卷积效果不好
R_temp = log(S+1) - log(L+1); % 去除光照影响
R = mat2gray(exp(R_temp)); % 线性拉伸
end
上述过程只能处理单通道的图像，对于三通道 RGB 图像，每一个通道做一遍即可获得 SSR 增强图像。
多尺度 Retinex
原理
可以看出单尺度 Retinex 最大的缺陷在于图像的增强效果非常依赖于高斯函数的 $sigma $ ，即所谓的 Retienx 尺度，$sigma $ 取值越小，高斯模糊效果越弱，光照估计的越不准，图像细节保留的越多，于是处理后的图像细节较好，但色彩严重失真；取值越大，则图像细节丢失，光照估计的越准确，处理后的图像颜色越自然，但细节丢失较多，且容易出现光晕现象。
所以为了弥补单一尺度下的这种缺陷，采取多个尺度的 SSR 进行处理，然后再对处理后的多个图像做加权求和，一定程度上实现细节和颜色方面的效果兼容。
这里的 (R) 同样也是仅针对的是一个通道的处理，多通道的重复使用即可。
实现
function [R,R_temp] = ssr(S,sigma)
% S = 0 ~ 255 , 三通道彩色图像
% Si = Ri×Li
% log(Ri) = log(Si) - log(Li)
% = log(Si) - log(S*Li)
% 结果初始化
R = zeros(size(S));
R_temp = zeros(size(S));
for i = 1:size(S,3)
Si = S(:,:,i); % 提取单通道图像
[~,Ri,~]= ssr4gray(Si,sigma); % 单通道 ssr
R_temp(:,:,i) = Ri; % 单通道结果存储
R(:,:,i) = mat2gray(Ri); % 单通道结果存储
end
end
其中，使用到了前面已经写好的单通道单尺度 Retinex 函数 ssr4gray
当然，我们也可以将图像转换到 HSV 颜色空间，使用多尺度 Retinex 进行增强
function R = msr4hsv(S,sigma)
% S = 0 ~ 255 , 三通道彩色图像
% R = 0 ~ 1 , 三通道彩色图像
% RGB --> HSV
[h,s,v] = rgb2hsv(S);
V = zeros(size(v)); % 增强后亮度分量初始化
for i = 1:length(sigma)
[~, Vi]= ssr(v,sigma(i)); % 多尺度 ssr
V = V + Vi/length(sigma); % 加权求和
end
V = mat2gray(V); % 归一化
% HSV --> RGGB
R = hsv2rgb(cat(3,h,s,V));
end
带颜色恢复的多尺度 Retinex
原理
根据前面介绍，Retienx 算法一直致力于单通道图像的处理，对于其他通道的处理都是同样的操作，这使得通道间的分量比例明显不协调，从而失去了真彩色，图像颜色不能令人满意，于是再在多尺度 Retinex 的基础上加入了各通道之间的色彩比例调节系数，以此来恢复色彩比例。
对于第 (i) 个通道图像 (S_i(x,y)) 其通道色彩恢复调节系数 (C_i(x,y)) 被表示为
于是，将其与之前的多尺度 Retinex 增强图像进行相乘，调整像素取值，实现通道间色彩比例调节。
这里的 (R_j) 表示第 (j) 个色彩通道。
实现
function [R,R_temp] = msrcr(S,sigma,beta,alpha,G,b)
% S = 0 ~ 255 , 三通道彩色图像
% R = 0 ~ 1 , 三通道彩色图像
[~,R_msr] = msr(S,sigma); % 多尺度 msr
% 结果初始化
R = zeros(size(S));
R_temp = zeros(size(S));
for i = 1:size(S,3)
Si = S(:,:,i); % 提取单通道图像
Sa = sum(S,3); % 按通道求和
Ci = beta*(log(alpha*Si+1) - log(Sa+1)); % 色彩恢复系数
R_temp(:,:,i) = G*(R_msr(:,:,i).*Ci+b); % 单通道恢复结果存储
R(:,:,i) = mat2gray(R_temp(:,:,i)); % 单通道结果存储
end
end
其中，使用到了前面已经写好的单通道多尺度 Retinex 函数 msr
基于二维伽马函数的光照不均匀图像自适应校正算法
二维伽马函数
我们知道，传统的伽马矫正，就是在原图像的基础上，(gamma) 次幂运算，(gamma) 若大于 1 则图像整体会被变暗，等于 1 则保持不变，小于 1 则图像整体变亮。
这就导致，图像是整体变化的，缺乏局部像素的考虑，因此可以针对每一个像素设置一个合适的 (gamma) 值，自适应地调整图像亮度。在该文章中，提出了一种新的二维伽马函数（文中公式写反了）。
即，根据根据多尺度 Retinex 方法，计算出光照分量 (L)
然后依据光照分量 (L) 计算每一个像素值的 (gamma) ，再将其应用在要处理的图像 (S) 上，进行逐像素伽马矫正。
这里，论文选择将图像转换到 HSV 颜色空间，然后针对亮度分量 (V) 使用上述方法增强后，还原图像，得到RGB增强图像。
实现
function img_out = gamma2d(img_in)
% rgb --> hsv
img_in = double(img_in);
[h,s,v] = rgb2hsv(img_in);
% 多尺度光照估计
[~,~,l1] = ssr4gray(v,15);
[~,~,l2] = ssr4gray(v,80);
[~,~,l3] = ssr4gray(v,250);
l = (l1+l2+l3)/3;
% 自适应二维伽马矫正
% m = mean2(l);
m=128;
r = (1/2).^((m-l)./m);
v1 = (v/255).^r;
% v1 = mat2gray(v1);
% hsv --> rgb
img_out = hsv2rgb(cat(3,h,s,v1));
end
论文整体代码
论文首先绘制了亮度图像，然后给出了二维伽马函数的曲线对比图，以及三个不同场景下的图像经直方图均衡化（HE），带颜色恢复的多尺度 Retinex ，以及二维伽马函数矫正的方法。下面是论文插图的整体实现。
clear;close all;clc
addpath('./model');
addpath('./src');
%% 图2
figure(2)
img_gray = imread('Fig0310(a)(Moon Phobos).tif');
subplot(1,2,1);imshow(img_gray);xlabel('(a) 灰度图像')
img_gray = double(img_gray);
[~,~,l1] = ssr4gray(img_gray,15);
[~,~,l2] = ssr4gray(img_gray,80);
[~,~,l3] = ssr4gray(img_gray,250);
img_gray_l = (l1+l2+l3)/3;
img_gray_l = wcodemat(img_gray_l,255);
subplot(1,2,2);imshow(uint8(img_gray_l));xlabel('(b) 光照分量')
%% 图2
figure(3)
img_rgb = imread('彩色图像.png');
subplot(1,2,1);imshow(img_rgb);xlabel('(a) 彩色图像')
[h,s,v] = rgb2hsv(img_rgb); % 提取亮度分量
[~,~,l1] = ssr4gray(v,15);
[~,~,l2] = ssr4gray(v,80);
[~,~,l3] = ssr4gray(v,250);
img_rgb_l = (l1+l2+l3)/3;
img_rgb_l = wcodemat(img_rgb_l,255);
subplot(1,2,2);imshow(uint8(img_rgb_l));xlabel('(b) 光照分量')
%% 图4
figure(4)
F = 0:10:255;
I = [0, 64, 128, 192, 255];
m = 128;
func = @(F,I) 255*(F/255).^((1/2).^((m-I)/m));
lineset = {
'r--s','$L(x,y)$=0';
'g-o','$L(x,y)$=64';
'b-*','$L(x,y)$=128';
'c-x','$L(x,y)$=192';
'k-p','$L(x,y)$=255';
};
for n = 1:length(I)
plot(F,func(F,I(n)),lineset{n,1},'linewidth',1,'MarkerSize',5,'MarkerEdgeColor','k','MarkerFaceColor',lineset{n,1}(1));
hold on
end
legend(lineset(:,2),'location','NorthWest','box','off','interpreter','latex')
xlabel('输入图像亮度');ylabel('校正后输出图像亮度');
% axis square
axis tight
%% 模型对比
models = {
'(a)原图像',@(x) x;
'(b)HE算法处理结果', @(x) he(x);
'(c)MSRCR算法处理结果',@(x) msrcr(double(x),[15,80,250],46,125,192,-30);
'(d)本文算法处理结果',@(x) gamma2d(x);
};
%% 图7
figure(7)
img = imread('场景1.png');
for i = 1:length(models)
model = models{i,2};
img_en = model(img);
subplot(3,size(models,1),i);
imshow(img_en);
xlabel(models{i,1});
end
%% 图8
% figure(8)
img = imread('场景2.png');
for i = 1:length(models)
model = models{i,2};
img_en = model(img);
subplot(3,size(models,1),i+4);
imshow(img_en);
xlabel(models{i,1});
end
%% 图9
% figure(9)
img = imread('场景3.png');
for i = 1:length(models)
model = models{i,2};
img_en = model(img);
subplot(3,size(models,1),i+8);
imshow(img_en);
xlabel(models{i,1});
end
这里，文中并没有取均值作为 (m) ，而是直接使用了 128 作为固定的 (m) 。并且在求光照分量时，也达不到图像那样的高斯模糊效果（只有当高斯模糊的卷积用 imfilter 实现时，才会和论文中那样，但复现的结果表明作者使用了频域乘积去计算时域卷积，所以到底怎么回事？）
这里的所有图像，均来自论文截图，并未找到对应的标准测试图像，可能是作者自己拍的吧。
参考文献
- 刘志成,王殿伟,刘颖,刘学杰.基于二维伽马函数的光照不均匀图像自适应校正算法[J].北京理工大学学报,2016,36(02):191-196+214.
- 何艳. 基于Retinex图像增强的研究与应用[D].合肥工业大学,2014.
- 闫保中,韩旭东,何伟.基于Retinex理论改进的低照度图像增强算法[J].应用科技,2020,47(05):74-78.
- 刘远仲,张强,杨嘉,谭鹤毅,张海波,唐天国.一种改进的低照度图像增强算法研究[J].四川轻化工大学学报(自然科学版),2020,33(06):46-52.
- 姜雪松,姚鸿勋.夜晚图像增强方法综述[J].智能计算机与应用,2020,10(03):394-403.
- 刘军. 基于Retinex理论的彩色图像增强技术研究[D].中国科学院研究生院(西安光学精密机械研究所),2015.
- 林宝栋. 基于人眼视觉特性的低照度图像增强算法研究[D].南京邮电大学,2017.
- 严亮. 基于Retinex理论的图像增强技术研究[D].华中科技大学,2016.
- 王浩,张叶,沈宏海,张景忠.图像增强算法综述[J].中国光学,2017,10(04):438-448.
- 陈雾. 基于Retinex理论的图像增强算法研究[D].南京理工大学,2006.
- 肖晓. 基于Retinex的图像增强算法研究与实现[D].四川师范大学,2016.
- 宋瑞霞,李达,王小春.基于HSI色彩空间的低照度图像增强算法[J].图学学报,2017,38(02):217-223.
- 王殿伟,王晶,许志杰,刘颖.一种光照不均匀图像的自适应校正算法[J].系统工程与电子技术,2017,39(06):1383-1390.
- 秦钟,杨建国,王海默,杨佳睿,崔春晖.基于Retinex理论的低照度下输电线路图像增强方法及应用[J].电力系统保护与控制,2021,49(03):150-157.
心得
输入输出数据类型要当心
在做图像处理时，要十分注意数据类型的转换，图像一般取值 0 - 255 或者 0 - 1 ，但很多情况下，论文或者代码对输入的图像数据取值范围是有条件的，必须弄清楚，否则效果或出现很大的异常，这也是为什么我要在函数输出时，写一个 R_temp，接住这个中间计算结果。因为在其他函数调用该函数时，我们希望它不是经过转型或者归一化处理的，这样才能保证数据计算一致性。
函数功能单一化、可扩展很重要
从原理上来讲，Retinex 系列算法之间是存在递进关系的，因此在编写函数时，尤其要关注算法的输入输出，既保证函数编写完成后，能够实现自身功能的同时，也能为其他调用它的函数提供模块化的处理，这样一层套一层，代码变得很结构化、看起来会比较有层次。这些都需要前期对函数功能进行仔细地分析。
句柄和元胞是个好东西
本文在复现效果对比时，使用了句柄和元胞来存放一些同类型的函数或变量，这样既方便了修改，也减少了代码的冗余度，个人觉得这个做法值得推广。