retinex 的水下图像增强算法_Retinex图像增强算法
前一段时间深入学习了一下图像增强算法后发现Retinex理论在彩色图像增强图像去雾以及彩色图像恢复等方面表现卓越下面将介绍我对该算法的见解
Retinex理论
Retinex理论源于Land及McCann在20世纪60年代的一系列重要贡献。其核心理念指出，在感知某一特定点的颜色与亮度时，并非仅由该点进入人眼的绝对光线决定，而与其周边区域的颜色亮度水平密切相关。‘Retinex’这一术语由‘视网膜’（Retina）与‘大脑皮层’（Cortex）两个关键词组合而成。为了表达他对视觉系统特性的认识仍有待进一步明确。
Land的Retinex模型是建立在以下的基础之上的：
一、真实世界的本体没有色彩，人类感官所感知的颜色源于光与物质相互作用的过程。纯净水呈无色状态，在此情况下观察到泡沫状的水膜呈现斑斓色彩，则是因为薄膜表面发生光波干涉现象的结果。
二、每一颜色区域由给定波长的红、绿、蓝三原色构成的；
三、三原色决定了每个单位区域的颜色。
根据Retinex理论的核心观点，在确定颜色时，并非依赖于反射光强度的具体数值（绝对值），而是基于物体对不同色光（红、绿、蓝）的不同反射能力。该理论指出，在色彩感知方面与光源分布不均匀无关。其呈现一致性的特征。如图所示，请注意以下内容：观察到的对象图像S是由于表面对入射光L进行反射而产生的。其中R代表表面相对于光源L的比例关系（即反映比率），这一参数完全由表面特性决定而不受光源变化影响
图1.Retinex理论中图像的构成
Retinex理论的核心理论基础认为原始输入图像S是由光源反射图L和表面反射特性图R的乘积构成,其数学模型可表示为以下公式：
旨在通过推导原始图像S中的光照参数L来分离出反射成分R，并通过补偿光照不均匀性来优化图像的视觉效果。这一过程与人类视觉系统的方式相似，在实际操作中，默认情况下会将输入图像转换为对数域进行处理。
，从而将乘积关系转换为和的关系：
Retinex的核心原理在于计算图像中的亮度值 L 。该方法通过从图像 S 中确定亮度分量 L ，并去除这些亮度影响来获得原始反射比值 R 。具体而言，在这一过程中我们首先从图像 S 中确定亮度分量 L ；接着通过去除这些亮度影响来获得原始反射比值 R 。
函数
实现对照度L的估计(可以去这么理解，实际很多都是直接估计r分量)。
Retinex理论的理解
若读者关注学术论文，在后续内容中必然涉及介绍两个经典的 Retinex 算法：即基于路径的 Retinex 以及基于中心/环绕 Retinex（Central/Cascaded）。然而，在阐述这两个经典算法之前, 我可以分享一些个人见解以帮助初次接触该理论的朋友更快理解其核心概念。如有误, 请各位指正。
我对Retinex理论的理解是类似于降噪的技术，在计算机视觉领域有着重要的应用价值。其核心在于对图像构成做出合理的假设：观察到的图像可被视为由入射光与乘性噪声相乘得到的结果，在这种模型下，默认情况下认为入射光是一种具有乘性特性的、相对均匀且变化缓慢的空间域噪声。基于这一前提构建起来的方法论框架，则是 Retinex 算法所追求的目标——准确估计并消除各像素位置上的残留噪声。
在极端条件下，我们可以合理地假设整幅图像中的各个分量呈现均匀分布状态。基于此设想，在将图像转换为对数域后计算其均值是一种直接的方法。基于此设想，我提出了一种验证方案来验证自己的猜想。流程如下：
(1) 将图像变换到对数域
；
(2) 归一化去除加性分量
；
(3) 对步骤3得到的结果求指数，反变换到实数域
。
这里为了简化描述，省略了对图像本身格式的变换，算法用Matlab实现：
% ImOriginal:原始图像
%type:'add' 指明了成分的叠加特性，在图像场景中得到体现；
'mult' 则指明了成分间的乘积特性，在亮度计算中得到应用。
[m,n,z] = size(ImOriginal);ImOut = uint8(zeros(m,n,z));for i = 1:z
if strcmp(type,'add')
ImChannel = double(ImOriginal(:,:,i))+eps; elseif strcmp(type,'mult')
ImChannel = log(double(ImOriginal(:,:,i))+eps);else
error('type must be ''add'' or ''mult''');end
ImOut(:,:,i) = EnhanceOneChannel(ImChannel);end
ImOut = max(min(ImOut,255), 0);end
function ImOut = EnhanceOneChannel(ImChannel)
% 计算计算单个通道的反射分量
%1.对全图进行照射分量估计
%2.减去照射分量
%3.灰度拉伸
计算 ImChannel 等于其自身元素与该矩阵所有元素的最大值的商。
经过处理后得到 ImRetinex 作为指数函数的结果。
将结果转换为 uint8 类型。
end
旨在检验该算法的有效性。为了比较该经典Retinex算法与其他方法之间的性能差异，在这里我们采用了对比试验的方法。结果显示该对比试验表明
图2.测试原图
图3.经典Retinex算法结果
图4.上述方法结果
通过对比可以看出，在去除非均匀光照方面该算法仍有一定的可行性，并且能够有效避免边缘模糊现象。其缺点在于，在去除非均匀光照分量L的过程中会保留反射成分R，在实际应用中我们采用归一化处理后再进行反变换运算。这一步骤的作用相当于消除了一个全局均值亮度偏移。由于整个过程采用了全局操作，默认假设所有位置上的光照强度分布是相同的，在灰度拉伸过程中未能充分考虑局部细节的变化而导致整体亮度偏低。值得注意的是这种方法虽然能在一定程度上提高图像质量但全局化的光照估计策略仍然限制了其在细节增强方面的表现因此在色彩恢复和亮度提升方面仍存在一定的局限性
我认为Retinex算法的核心在于准确地分析噪声特性。很多人都注意到基于Retinex技术的水下增强和去雾效果，并且我也很好奇便尝试了一下。大雾照片确实不少，在最近一次大雾天气里随手拍了几张照片后终于派上了用场，请看对比图
图5.有雾原图
图6.经典Retinex去雾效果
图7.上述方法去雾效果
还是老方法，
此时通过对比试验来验证效果最为直接，
因此特意选择了幅极具干扰性的图像作为测试样本。
对于显示器质量稍逊的用户而言，
在原始图像中难以辨识雾气后的细节内容。
就去雾效果而言，
与传统经典算法相比，
本实验的结果在主观感受上依然令人满意。
在上述案例中核心在于准确解析有雾图像的结构特征与Retinex理论一开始的核心思想存在明显不同之处即在线性叠加模型下针对这种残留性干扰成分传统Retinex算法仅能实现对其残留性干扰成分的估算而未能有效解决这一问题值得注意的是无论从基本理论框架还是具体实现细节上都应摒弃传统的对照度和反射率主导型方法进而最为关键的是对这一残留性干扰成分进行精确估算
对于有雾图像而言，在这种情况下我们可以将其视为通过磨砂玻璃观察到清晰图像的过程。这样大家就能很好地理解为什么认为在这种情况下将雾气干扰视为加性影响的原因。至于后面提到的经典算法，在本质上都是通过原图像中的像素点来推断原始照度值。从上述程序运行结果来看，在全局估计方法中对局部增强效果较差的情况较为常见：当存在照度不均匀（如雾气浓度分布不均）或背景颜色亮度差异较大时会导致处理效果显著下降。
然而，在经典方法中也存在明显的局限性。通过查看图3可知，在实际应用中经典的图像增强算法往往会产生明显的光晕效应（例如，在处理教科书中深色区域时会出现蓝色背景带）。由此可见，在解决光照估计问题以及去除光晕现象方面又衍生出了众多基于Retinex理论的变体算法。这些 Retinex 理论的变体算法主要针对的是光照估计问题以及去除光晕现象等具体应用场景。接下来将重点介绍其中两个具有代表性的经典算法，请关注后续内容更新！如需进一步了解这些算法的具体实现细节，请参考附录中的 Matlab 代码部分。
该算法由McCann与Frankle共同发展而来。该方法遵循多重迭代策略。每个像素点的颜色值由特定环路的影响决定。通过反复迭代逐步趋近于理想数值。通过对螺旋形路径上各个像素点亮度的变化进行分析计算，并结合多次迭代的结果，能够有效估算并消除图像中的光照偏差。
图8.McCann Retinex算法路径选择示意图
这张图取自某篇参考文献。不过我打算按照论文中的方式来讲解它，因为它内容较为复杂且不易于理解。通过观察该图可以看出，在算法设计中我们采用螺旋式路径选择用于估计的关键点，并且靠近预测中心区域的像素点被选中的数量更多。这在实际应用中也是合理的。
该算法估测图像经过以下几个步骤：
1. 将原图像变换到对数域
，彩色图像对各通道都进行对数变换；
2. 初始化常数图像矩阵
，该矩阵作为进行迭代运算的初始值；
3. 计算路径，如上图8所示，这里令
为路径上的点，从远到近排列；
4. 对路径上的像素点按照如下公式运算：
公式所表示的大致意思为：由远及近的过程中,中心点像素值减去路径上的像素值得到的差值的一半量与前一时刻的估计值之间形成了总和.最终而言,在图像处理中,中心像素点的像素通常表现为这样的形式.
其中
表示中心位置最终的反射率估计，
为常数值为转换后的图像中的最大值，在步骤2中被确定。从这里将
基于Retinex理论对图像进行分解后可得，通过查看最终公式可以看出所有非照度分量已被去除。其中中心位置处的反射率则由路径上各点反射率差异计算得出。此外，在观察轨迹时会发现距离越近的位置在估算过程中占据的比例越大。
从我的角度来看，这类估计算法其本质即为中心像素与其周边像素亮度值之差来推断该像素处反射率的变化情况。在实际应用中发现，在以下两种情况下表现尤为突出：首先，在图像区域整体明亮的情况下（即当该区域的平均亮度较高时），算法能够较为准确地预测出同样呈现高亮特征的结果；其次，在区域平均亮度偏低的情况下，则该区域与邻近区域间的对比度必然为负值（即出现反向变化），此时预测结果同样呈现高亮特征。
其亮度相对较低,这就要求路径上的各个点必须确保能够充分反映整幅图像的独特特征.假设输入图像具有高度规律性,则可能导致中心像素周围的计算结果不足以准确反映该像素自身的灰度特性.可见于图9,算法其实质是通过计算中心区域与周边区域之间的差异来实现.黑色区域中的图像由于其与周边区域差异较小而呈现出较高的亮度,而随着距离逐渐增大差异也随之减小,最终导致无法准确体现原本黑色像素点处应有的亮度特性.
图9.算法的测试图(来自Retinex in Matlab)
从原始算法中可以看出，在其中还有一个关键步骤。这个过程即为迭代过程。其主要功能是尽可能地保留中央区域及其周边区域的能量分布特征，在每次运算过程中只会将中央区域及其周边区域的能量分布值减半（如步骤4中的除以2操作）。经过反复迭代次数后能够保持下来的主要分量大约是原始能量分布特征值的\frac{1}{2^{n}}倍（其中n表示迭代次数）。
通过这种方式会使局部区域的信息更加丰富，并非单纯地降低动态范围的压缩程度。类似地，在进行这一操作时能够有效提升图像的表现力的同时也会带来一定的计算负担即运算量会有所上升这可以通过下图的具体数据对比来直观体现不同迭代次数对收敛速度的影响
图10.迭代次数的影响
为了方便各位自己研究，下面我给出该算法源码供大家参考：
function Retinex = retinex_frankle_mccann(L, nIterations)
The %RETINEX_FRANKLE_McCANN% variable calculates the raw Retinex output from an intensity image according to a specific algorithm or method.
% original model describedin:
Franklin A. Frankle and James P. McCarthy ("FRANKLE AND McCARTHY," "METHOD AND APPARATUS FOR LIGHTNESS IMAGING")
%INPUT:L - logarithmic single-channel intensity image to be processed
% nIterations - number of Retinex iterations
%
%OUTPUT:Retinex - raw Retinex output
%
%NOTES: - The input image is assumed to be logarithmic and in the range [0..1]
% - To obtain the retinex"sensation"prediction, a look-up-table needs to
% be applied to the raw retinex output
% - For colour images, apply the algorithm individually for each channel
%
%AUTHORS: Florian Ciurea, Brian Funt andJohn McCann.
% Code developed at Simon Fraser University.
%
%关于codesee的信息由Florian Ciurea、Brian Funt和John McCann提供
%"Retinex in Matlab,"by Proceedings of the IS&T/SID Eighth Color Imaging
%Conference: Color Science, Systems and Applications, 2000, pp 112-121.
%
% paper available online athttp://www.cs.sfu.ca/~colour/publications/IST-2000/
%
/* 版权信息：*/
此软件仅供研究性和教育用途复制和使用授权；不得出售该软件；本软件仅限于内部测试版或研究用途。
% redistributed so long as the original sourceandauthors are cited.
global RR IP OP NP Maximum
RR = L;Maximum = max(L(:));% maximum color value in the image
[nrows, ncols] = size(L);shift =2^(fix(log2(min(nrows, ncols)))-1);% initial shift
OP = Maximum*ones(nrows, ncols);% initialize Old Product
while (abs(shift) >=1)
for i =1:nIterations
CompareWith(0, shift);% horizontal step
CompareWith(shift, 0);% vertical step
end
shift = -shift/2;% update the shift
end
Retinex = NP;function CompareWith(s_row, s_col)
global RR IP OP NP Maximum
IP = OP;if (s_row + s_col > 0)
IP((s_row+1):end, (s_col+1):end) = OP(1:(end-s_row), 1:(end-s_col)) + ...
RR((s_row+1):end, (s_col+1):end) - RR(1:(end-s_row), 1:(end-s_col));else
IP(1:(end+s_row), 1:(end+s_col)) = OP((1-s_row):end, (1-s_col):end) + ...
RR(1:(end+s_row),1:(end+s_col)) - RR((1-s_row):end, (1-s_col):end);end
IP(IP > Maximum) = Maximum;% The Reset operation
NP = (IP + OP)/2;% average with the previous Old Product
OP = NP;% get ready for the next comparison
McCann99 Retinex算法
McCann 99Retinex 算法本质上与 McCann Retinex 算法是相同的，在 McCann 99 算法中不再以螺旋路径选择估计像素点而是采用了图像金字塔分层选取像素的方式 算法同样通过依次取样 比较评估并计算平均值完成迭代计算过程
图像金字塔的概念较为基础，即通过降采样的方法对图像进行处理，以多层次的表现形式呈现其细节信息。金字塔顶端层具有最低分辨率，而底部层通常保留原始图像是最高分辨率的部分。
图像金字塔的概念较为基础, 即通过降采样的方法对原始图片进行处理, 以多层次的表现形式呈现其细节信息. 金字塔顶端层具有最低分辨率, 而底部层通常保留原始图像是最高分辨率的部分.
图11.图像金字塔示意图(来自网络)
基于图1的展示，McCann99算法通过逐层分析每个像素及其周边8个像素的关系来评估其反射率分量。
在前一层计算结束后通过插值运算对估计的反射率进行分类处理,从而实现上一层结果的优化
在金字塔层中的图像是与其下方一层具有相同尺寸，并且再次执行同样的运算过程；最后能得到结果的就是增强后的图像是正确的。其中，在上一节步骤4中描述了相减后求平均这一操作。
因此，McCann99算法，此处简化描述为以下几步：
1. 将原图像变换到对数域
，彩色图像对各通道都进行对数变换；
2. 初始化(计算图像金字塔层数；初始化常数图像矩阵
，该矩阵作为进行迭代运算的初始值)；
从顶层结构开始,一直到最下层执行8邻域对比计算.该计算方式与经典的MccCann Retinex算法一致,请参阅上一节所述的步骤4.
4. 第n层运算结束后对第n层的运算结果
进行插值，变成原来的两倍，与n+1层大小相同(此处默认n越大越靠近底层)；
5. 当最底层计算完毕得到的
即最终增强后的图像。
为了方便各位自己研究，下面我给出该算法源码供大家参考：
function Retinex = retinex_mccann99(L, nIterations)
global OPE RRE Maximum
[nrows ncols] = size(L);% get size of the input image
通过调用ComputeLayers函数计算给定行数和列数得到金字塔层数nLayers。此代码段用于计算金字塔结构中的层数参数
nrows = nrows/(2^nLayers);% size of image to process for layer 0
cols = cols / pow(2, n_layers); if (row_count * col_count > MAX_IMAGE_SIZE) % Exceeding image size for initial layer throw an error indicating invalid dimensions for the first layer
end
Maximum = max(L(:));% maximum color value in the image
OP = Maximum*ones([nrows ncols]);% initialize Old Product
for layer = 0:nLayers
RR为图像分辨率调整后的结果，并通过ImageDownResolution函数从输入L中计算得到；指数项2^(nLayers-layer)用于确定调整的比例或层次。注释说明：该操作旨在将输入数据缩减至所需处理的层规模。
OPE = [zeros(nrows,1) OP zeros(nrows,1)];% pad OP with additional columns
OPE = [zeros(1,ncols+2);OPE; zeros(1,ncols+2)]; % and rows
RRE = [RR(:,1) RR RR(:,end)];% pad RR with additional columns
RRE = [RRE(1,:);RRE; RRE(end,:)]; % and rows
for iter =1:nIterations
CompareWithNeighbor(-1, 0);% North
CompareWithNeighbor(-1, 1);% North-East
CompareWithNeighbor(0, 1);% East
CompareWithNeighbor(1, 1);% South-East
CompareWithNeighbor(1, 0);% South
CompareWithNeighbor(1, -1);% South-West
CompareWithNeighbor(0, -1);% West
CompareWithNeighbor(-1, -1);% North-West
end
NP = OPE(2:(end-1), 2:(end-1)); OP = NP(:, [fix(1:0.5:ncols) ncols]);%%% these two lines are equivalent with
OP = OP([fix(1:0.5:nrows) nrows], :);%%% OP = imresize(NP, 2) if using Image
nrows = 2nrows;ncols = 2ncols; % Processing Toolbox in MATLAB
end
Retinex = NP;function CompareWithNeighbor(dif_row, dif_col)
global OPE RRE Maximum
% Ratio-Product operation
IP = OPE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col)) + ...
RRE(from 2 to end - 1, from 2 to end - 1) subtracts RRE(from 2 + dif\_row to end - 1 + dif\_row, from 2 + dif\_col to end - 1 + dif\_col); the Irradiance Penquet is set to Maximum when its value exceeds Maximum;// This line ensures the reset operation is performed.
% skip computing outputs in row or column with undefined neighbors
if (dif_col == -1) IP(:,1) = OPE(2:(end-1),2);end
if (dif_col == +1) IP(:,end) = OPE(2:(end-1),end-1);end
if (dif_row == -1) IP(1,:) = OPE(2, 2:(end-1));end
if (dif_row == +1) IP(end,:) = OPE(end-1, 2:(end-1));end
NP = (OPE(2:(end-1),2:(end-1)) + IP)/2;% The Averaging operation
OPE(2:(end-1), 2:(end-1)) = NP;function Layers = ComputeLayers(nrows, ncols)
power =2^fix(log2(gcd(nrows, ncols)));% start from the Greatest Common Divisor
while(power > 1 & ((rem(nrows, power) ~= 0) | (rem(ncols, power) ~= 0)))
power = power/2;% and find the greatest common divisor
end % that indicates it's a power of two Layers = log2(power)
[rows, cols] = size(A);% the input matrix A is viewed as
result_rows = rows/blocksize;% a series of square blocks
result_cols = cols/blocksize;% of size = blocksize
Result = zeros([result_rows result_cols]); for crt_row从1到result_rows %然后每个像素则按以下方式计算
for crt_col =1:result_cols % the average of each such block
Result(crt_row, crt_col) 是由参数 crt_row 和 crt_col 所决定的计算结果，并赋值给 mean2 函数作用于数组 A 中从 1+(crt_row-1)blocksize 到 crt_rowblocksize 的行以及从 1+(crt_col-1)blocksize 到 crt_colblocksize 的列的位置。
en
综上所述，在上述两个经典的估计算法的基础上，在经过反对数变换后就可以恢复为原始格式的图像。通过实验结果表明，在这些算法中仍存在一定的缺陷，在对光照度估计方面后续还存在许多改进方法如：单尺度Retinex (SSR)、多尺度Retinex (MSR)以及可变框架的Retinex等；此外还有针对增强后的光晕现象首先进行光照度分割然后再进行相关计算的方法等；这些改进方法本质上大同小异。鉴于篇幅限制以及详细讨论超出了本文范围，在此不做深入探讨。如有疑问或发现文中错误，请随时指正。