OpenCV3视觉对比敏感度（CSF）

人类视觉系统特性

对比敏感度函数

基于CSF模型的滤波

人类视觉系统特性

经过生物学家和心理学家长期的探索和研究，通过结合生物学、心理学和简单的人眼视觉模型等理论，发现人类视觉系统（Human Visual System，HVS）的视觉感知与人的心理反应具有紧密联系。由于客观图像质量评价的最终目的是设计出一种评估结果与人类感知高度吻合的数学模型，因而，人类视觉心理学特性的引入具有举足轻重的作用。

多通道分解特性

视觉皮层对不同类型的视觉信息有着各不相同的敏感性。人的视觉感知是一个多通道协同工作的过程，每个通道处理不同的视觉激励，通过相互联系以获得最佳的视觉效果。鉴于各通道的带通性，可以设计一组多尺度多方向的带通滤波器来近似模拟视觉感知的多分辨特性，对输入的视觉信息先在不同的频率和方向上进行分解，然后分配到相应的处理通道。如小波变换、Gabor变换和DCT变换等，已经被广泛应用于图像处理技术中。图像的高频部分表征视觉信息的边缘和结构特征；中频信息最易被区分，从而形成视觉关注点；低频部分提供细节较多的纹理信息。

掩膜效应

HVS在感知时接收到的刺激往往不止一种，而是多种刺激相互作用、彼此联系和互相影响的综合过程。当多种视觉激励信号在空间或时间上同时作用于视觉系统，一种视觉激励信号的出现会削弱或增强其他一种或多种激励的响应，这种现象就是视觉掩膜效应（Visual Masking Effects）。掩膜效应在视觉观察中最典型的现象是，图像中平滑区域的失真由于背景参照物的对比度低而易于发觉，复杂纹理区域的失真因为背景参考物对比度高而不易发现。

亮度特性

相较于激励的绝对亮度，HVS对激励与参照物之间亮度相对差异的敏感性更强。这个特性被称之为对比度，通常有两种基本定义方法：Weber-Fechner 定律和 Michelson 定律。Michelson 定律的表达式为 $C_{M}=\frac{L_{max}-L_{min}}{L_{max}+L_{min}}$ ，其中， $L_{max}$ 和 $L_{min}$ 分别表示图像中的最大和最小亮度。由Weber-Fechner定律表达式可推导出主观亮度与客观亮度的关系为 $S_n=k\textup{\textrm{ln}}B+k_0$ ，表明主观感知亮度同客观亮度对数相关。

对比敏感度函数

对比敏感度函数（CSF）可以衡量HVS对各种视觉刺激频率的敏感性。它反映了在不同的空间频率上，人眼对目标亮度的辨别力不同。

从频域来看，CSF对空间频率的不同成分具有不同的敏感程度。在低频区域和高频区域的视觉感知敏感度明显下降。具体表现为人眼对于图像中表征中频信息的边缘、轮廓和结构部分更加感兴趣，对于低频的平滑区域和纹理区域的敏感性较弱。一种经典的CSF由Campbell F. W.定义为 $H_1(f)=2.6\times (0.192+0.114\times f)\times \textup{exp}(-(0.114\times f)^{1.1})$ ，其中， $f=\sqrt{f_x^2+f_y^2}$ 为空间频率， $f_x$ 和 $f_y$ 为水平方向的空间频率和垂直方向的空间频率。另一种由Mannos和Sakrison提出，并由Daly改进的CSF模型如下：

其中 f 为径向空间频率，单位是(c / deg)；θ∈[-π，π]表示角频率；λ=0.114； $f_\theta = \frac{f}{0.15\textup{cos}(4\theta)+0.85}$ 表示f的基于方向的修改，用来降低沿对角线方向的对比敏感度；H2(f,θ)是CSF模型的频率响应。

基于CSF模型的滤波

网上说CSF原理的很多，但是几乎没有完整实现代码。我搜了一下发现CSF模型的实现代码很少，最后只找到了建立后一种CSF模型（Mannos和Sakrison）的一段Matlab版本例子。这里我在VS2019+OpenCV3把它翻译成C++。

计算公式中的径向频率

const int x = in.rows;
const int y = in.cols;
const double nfreq = 32.0;
Mat gker = (Mat_<double>(7, 7) <<0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,0.0000, 0.0000, 0.0000, 0.0002, 0.0000, 0.0000, 0.0000,0.0000, 0.0000, 0.0113, 0.0837, 0.0113, 0.0000, 0.0000,0.0000, 0.0002, 0.0837, 0.6187, 0.0837, 0.0002, 0.0000,0.0000, 0.0000, 0.0113, 0.0837, 0.0113, 0.0000, 0.0000,0.0000, 0.0000, 0.0000, 0.0002, 0.0000, 0.0000, 0.0000,0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000);
Mat xplane(y, x, CV_64FC1);
Mat yplane(y, x, CV_64FC1);
Mat xplane0(y, x, CV_64FC1);
Mat yplane0(y, x, CV_64FC1);
Mat radfreq(y, x, CV_64FC1);
double xb = -x / 2.0 + 0.5, xe = x / 2 - 0.5, yb = -y / 2.0 + 0.5, ye = y / 2.0 - 0.5;
meshgrid(xplane, yplane, xb, xe, yb, ye);
xplane = xplane / y * 2 * nfreq;
yplane = yplane / y * 2 * nfreq;
xplane0 = xplane.clone();
yplane0 = yplane.clone();
multiply(xplane, xplane, xplane);
multiply(yplane, yplane, yplane);
for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)radfreq.ptr<double>(i)[j] = sqrt(xplane.ptr<double>(i)[j] + yplane.ptr<double>(i)[j]); // radial frequency

计算CSF，结果 $H_{CSF}(u,v)$ 表示基于DFT的 $H_2(f,\theta )$ ，其中u和v是DFT坐标

Mat angleplane(y, x, CV_64FC1);
Mat ICFfrq;
Mat csf(y, x, CV_64FC1);
Mat ICFcx(x, y, CV_64FC2);
double w = 0.7;for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)angleplane.ptr<double>(i)[j] = atan2(yplane0.ptr<double>(i)[j], xplane0.ptr<double>(i)[j]);
angleplane *= 4;
for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)angleplane.ptr<double>(i)[j] = cos(angleplane.ptr<double>(i)[j]);
angleplane = (1 - w) / 2 * angleplane;
for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)angleplane.ptr<double>(i)[j] += (1 + w) / 2;
divide(radfreq, angleplane, radfreq);
xplane = radfreq * 0.114;
yplane = xplane.clone();
for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)xplane.ptr<double>(i)[j] += 0.0192;
xplane *= 2.6;
for (int i = 0; i < y; i++) {for (int j = 0; j < x; j++) {yplane.ptr<double>(i)[j] = pow(yplane.ptr<double>(i)[j], 1.1);yplane.ptr<double>(i)[j] = -yplane.ptr<double>(i)[j];yplane.ptr<double>(i)[j] = exp(yplane.ptr<double>(i)[j]);}
}
multiply(xplane, yplane, angleplane);
std::vector<std::pair<int, int> > p;
for (int i = 0; i < y; i++)for (int j = 0; j < x; j++)if (radfreq.ptr<double>(i)[j] < 7.8909)p.push_back(std::make_pair(i, j));
for (auto iter = p.cbegin(); iter != p.cend(); iter++)angleplane.ptr<double>(iter->first)[iter->second] = 0.9809;
csf = angleplane.t();
for (int i = 0; i < csf.rows; i++)for (int j = 0; j < csf.cols; j++)csf.ptr<double>(i)[j] = floor(csf.ptr<double>(i)[j] * 10000.0 + 0.5) / 10000.0;

通过 $I_{CF}=f^{-1}(H_{CSF}(u,v)f(I))$ 在频域中执行CSF滤波，其中f(·)和其反函数分别表示DFT和IDFT。生成的滤波后图像为 $I_{CF}$ ，滤波后采用上面的7×7的σ=0.5高斯核gker对 $I_{CF}$ 滤波，最后把 $I_{CF}$ 归一化从而能够显示出来。

ICFfrq = fft2(in);
ICFfrq = fftshift(ICFfrq);
std::vector<Mat> csfv;
split(ICFfrq, csfv);
multiply(csfv[0], csf, csfv[0]);
multiply(csfv[1], csf, csfv[1]);
merge(csfv, ICFfrq);
ICFfrq = ifftshift(ICFfrq);
ICFcx = ifft2(ICFfrq);
Mat ICF = ireal(ICFcx); // 取实部
filter2D(ICF, ICF, CV_64FC1, gker); // 可选的滤波
// return ICF;