将军一去，大树飘零；
壮士不还，寒风萧瑟。
荆璧睨柱，受连城而见欺；
载书横阶，捧珠盘而不定。
[toc]
谱残差显著性检测原理和实现代码
Saliency Detection: A Spectral Residual Approach, 该论文是上交的侯晓迪2017年大三上学期在上交时发表在cvvr上的论文，本科就在cccv着实是大神，更为传奇的是论文核心就是5行代码，当前引用量已近4000。
论文介绍
从信息论的角度来看，图像中的信息$H(Image)$可以分为两部分：
\(H(\text{Image}) = H(\text{Innovation}) + H(\text{Prior Knowledge})\)
其中$H(\text{Innovation})$表示图片中的novelty新奇的部分，$H(\text{Prior Knowledge})$则表示了冗余的部分，在图像领域中，所谓冗余指的是环境中的统计不变性。因此要对图像做显著性检测，就要先去掉冗余信息。
log spectrum表示
如何表示图像中的信息呢？在图像中有一个$1/f$定理，也就是对于大量图片的平均傅立叶谱的幅值$\mathcal{A}(f)$满足以下分布：
\(E\{\mathcal{A}(f)\} \propto 1/f.\)
也就是平均傅立叶谱与频率成反函数的关系，根据下面公式：
\(\begin{aligned}
y &= \frac 1x \\
\log{y} &= \log{\frac 1x} \\
\log{y} &= - \log{x}
\end{aligned}\)
在log-log的尺度上，幅值谱近似的成了一条直线，也就是线性关系。
注意原文中的是：
On the log-log scale, the amplitude spectrum of the ensemble of natural images, after averaging over orientations, lies approximately on a straight line.
也就是这里所谓的平均傅立叶谱不仅是在多个图片上的平均，同时也就是对二维的各个方向频谱的平均。
log-log谱听起来不错，在其他场景上也有广泛的使用，但是由于log-log spectrum存在两个问题：
- 一是单个图像一般不具有尺度不变性；
- 二是在低频处过于稀疏，高频处过于密集。
从论文中给出的图可以看出log-log Spectrum存在上述问题，因此选用log spectrum进行代替（后面会介绍如何计算Log-log谱）。
图像频率值范围大（通常0到10^6，甚至更高），虽然计算机可以处理，图像的显示系统却难以再现如此大的灰度范围，因此许多重要的灰度细节在傅里叶频谱中会丢失（高频量支配显示，低频量丢失），进行对数变换后，频谱值范围被压缩，显示出的图像细节内容更丰富。
剩余谱spectral residual
下图是不同图像中的log谱，可以看到有很大的相似性，注意第二个图中下方的黑边，这是将所有的图像都填充到固定大小$64 \times 64$，这样方面快速计算fft及归一化比较。
因此，图像中的冗余部分$\mathcal{R}(f)$的信息可以用log谱来表示： \(H(\mathcal{R}(f)) = H(\mathcal{L}(f)|\mathcal{A}(f))\) 其中$\mathcal{L}(f)$为log谱，$\mathcal{A}(f)$为log谱中的冗余部分，也就是log谱中一般的形状（在大部分图像中都相似的部分）表示已知的先验信息。$\mathcal{A}(f)$是对log谱做平均滤波得到的： \(\mathcal{A}(f) = h_n(f) * \mathcal{L}(f)\) 其中$h_n(f)$为$n\times n$的均值滤波器： \(h_{n}(f) =\frac{1}{n^{2}}\left(\begin{array}{ccccc} 1 & 1 & 1 & \ldots & 1 \\ 1 & 1 & 1 & \ldots & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & 1 & \ldots & 1 \end{array}\right) \\\) 不同大小滤波器的效果如下：
可以看出滤波器的大小对图像影响不大，因此论文里默认使用$3\times3$的滤波器。
因此对于剩余谱spectral residual $\mathcal{R}(f)$计算如下： \(\mathcal{R}(f) = \mathcal{L}(f) - \mathcal{A}(f)\) 剩余谱中包含了图片中的新奇部分，也就是显著部分。
显著图的计算公式
saliency map就是先在log谱空间内把冗余信息去除得到剩余谱，再将剩余谱从log空间幅值空间转换到频率谱的复数空间中，再利用傅立叶逆变换在原始空间中把图像重构输出。
也就是对应下面的五行公式： \(\begin{aligned} A(f) &=\operatorname{Amplitude}(\mathfrak{F}(\mathbf{x})) \\ P(f) &=\text{Phrase}(\mathfrak{F}(\mathbf{x})) \\ L(f) &=\log (A(f)) \\ R(f) &= L(f)-h_{n}(f) \cdot L(f) \\ S(\mathbf{x}) &= g(x)*\mathfrak{F}^{-1}\left[\exp (R(f)+i P(f))\right]^2 \end{aligned}\)
这里要重点看一下第5个公式，也就是如何从剩余谱恢复到原始图像空间变成显著图。
对于欧拉公式：
\(e^{ix} = \cos x + i \sin x\)
可变形为：
\(e^{a+bi} = e^a(\cos b + i \sin b)\)
对于某一频率的成分：
\(A e^{i \varphi} = A \cos \varphi + A \sin \varphi i\)
其中$A$为该成分对应的对应的振幅，$\varphi$为相位角，$i$为虚数单位
公式4中的剩余项$R(f)$是在log尺度的幅值余项，因此如果要恢复到原始图像空间中就需要先将这个幅值做指数exp恢复到原始的幅值空间。
频域对应的是复数空间，在保持相位不变的情况下，将幅值去除平均项后，要转换到对应的复数空间才能做傅立叶反变换转换到时域中。因此用新的幅值和原来的相位角得到的复数表示为：
\(\underbrace{e^{R(f)}}_{amplitude} \cdot e^{i \overbrace{P(f)}^{angle}} = e^{R(f)+iP(f)}\)
也就是第5个公式的$\exp(R(f)+iP(f))$部分
上式根据欧拉公式可以变换为： \(\begin{aligned} &= e^{R(f)}(\cos(P(f)) + i \sin(P(f))) \\ &= \textcolor{red}{e^{R(f)}\frac{\text{real}}{\text{amp}}} + \textcolor{green}{e^{R(f)}\frac{\text{real}}{\text{amp}}} i \end{aligned}\) 在很多开源实现中实际使用的正是该公式。
下图为在复平面上来看这个过程：
在计算出这一部分后，有一个平方操作$[]^2$，这里是为什么呢？我们知道剩余谱的内容也可以解释为图像中未预料到的部分，因此在逆变换恢复出来的显著性图后，图中每一点的值都被平方以表示估计误差。 同时为了获得更好的视觉效果，我们使用高斯滤波器$g(x)$来平滑显著图。
检测图像中的对象
这里的思路比较直接，根据重建出的显著性图，使用的简单的阈值分隔来根据显著性检测图像中的物体。
给定图像的显著图$S(x)$，目标图$O(x)$可根据如下计算： \(\mathcal{O}(x) = \begin{array}{ccc} 1 &\text{if}\ \mathcal{S}(x) > \text{threshold}, \\ 0 &\text{otherise}.\\ \end{array}\) 论文中根据经验一般可设置$\text{threshold } = E(S(x))\times3$，其中$E(S(x))$为显著图中的平均强度。
下面是一些检测的效果，可以看到在不同的尺度下的检测效果也不同。
实现代码及细节分析
log谱图的绘制
对于论文中二维图像相位谱的绘制搜到的博客里都在好奇，这里给出对应的代码。其主要的思路是：通过径向平均将二维频率对应的功率谱合并到一维的平均频率空间中，再通过histogram来进行统计
import cv2
import matplotlib.pyplot as plt
import numpy as np
WIDTH = 128
img = cv2.imread('./road.jpg', 0)
img = cv2.resize(img, (WIDTH, WIDTH*img.shape[0]//img.shape[1]))
frame = img.copy()
# expand the image to an optimal size for FFT
rows, cols = frame.shape[:2]
# 计算快速DFT的合适的size，DFT在size为二的n次幂时效果最快
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
# 将图像调整到合适的宽高
frame = cv2.copyMakeBorder(frame, 0, nrows-rows, 0, ncols - cols,
cv2.BORDER_CONSTANT, value=0)
img_dft = np.fft.fft2(frame)
# 计算二维Log功率谱，其中平方是功率谱，Log是便于观察高频区间
spectrum = np.log10((np.abs(img_dft))**2)
# 对二维频率做径向平均，对应的为平均频率
dists = np.sqrt(np.fft.fftfreq(frame.shape[0])[:, np.newaxis]**2 +
np.fft.fftfreq(frame.shape[1])**2)
# 展开到一维空间
dists = dists.ravel()
spectrum = spectrum.ravel()
# 计算直方图分布
L = max(frame.shape)
freqs = np.fft.fftfreq(L)[:L//2] # 重点关注的频率范围，根据当前采样率计算的
dcount = np.histogram(dists, bins=freqs)[0] # 当前频率空间中的
histo, bins = np.histogram(dists, bins=freqs,
weights=spectrum)
centers = (bins[:-1] + bins[1:]) / 2 # bin的中心
fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].plot(centers, histo/dcount, marker='.') # 落在每个径向频率空间dists里的log功率大小
ax[0].set(xlabel='frequency', ylabel='log-spectrum')
ax[1].plot(centers, histo/dcount, marker='.') # 落在每个径向频率空间dists里的log功率大小
ax[1].set(xlabel='log-frequency', ylabel='log-log spectrum', xscale='log')
plt.subplots_adjust()
plt.tight_layout()
可以看到，在log频谱中，在log-log频谱中基本上转换成了线性关系，也能看出在log-log谱的主要问题：低频部分很稀疏，高频部分过于密集。
对于上面代码里主要的步骤，我们分别来看一下。
频谱中两个维度分别的频率分布：
import seaborn as sns
freq_x = np.fft.fftfreq(frame.shape[0])
freq_y = np.fft.fftfreq(frame.shape[1])
plt.figure(figsize=(2, 0.1))
sns.heatmap(np.abs(freq_y).reshape(1, -1))
plt.figure(figsize=(0.1, 2))
sns.heatmap(np.abs(freq_x).reshape(-1, 1));
径向平均之后的频率：
dists = np.sqrt(freq_x[:, np.newaxis]**2 +
freq_y**2
)
sns.heatmap(dists);
将其展开用来统计直方图：
plt.figure(figsize=(20, .2))
sns.heatmap(dists.ravel().reshape(1, -1));
log功率谱，注意跟原始图像无关，该图为频域空间。
sns.heatmap(spectrum.reshape(frame.shape));
saliency map
首先看一下均匀滤波器的实现，在python中有很多方法，这里列出三种scipy.ndimage.uniform_filter、cv2.blur和cv2.boxFilter。
from scipy.ndimage import uniform_filter
x = np.random.random(size=(4, 4))
uniform_filter(x, size=3, mode='mirror')
array([[0.5111018 , 0.5432415 , 0.4226726 , 0.4548123 ],
[0.6235194 , 0.57643326, 0.53358509, 0.48649895],
[0.51063538, 0.52194351, 0.4569828 , 0.46829093],
[0.62305298, 0.55513527, 0.56789529, 0.49997758]])
cv2.blur(x, (3, 3))
array([[0.5111018 , 0.5432415 , 0.4226726 , 0.4548123 ],
[0.6235194 , 0.57643326, 0.53358509, 0.48649895],
[0.51063538, 0.52194351, 0.4569828 , 0.46829093],
[0.62305298, 0.55513527, 0.56789529, 0.49997758]])
cv2.boxFilter(x, -1, (3, 3))
array([[0.5111018 , 0.5432415 , 0.4226726 , 0.4548123 ],
[0.6235194 , 0.57643326, 0.53358509, 0.48649895],
[0.51063538, 0.52194351, 0.4569828 , 0.46829093],
[0.62305298, 0.55513527, 0.56789529, 0.49997758]])
注意上面cv2.blur和cv2.boxFilter的padding策略都是镜像填充，因此将uniform_filter也设置为mode='mirror'。
下面给出具体的代码：
c = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
# 平方为求功率谱
mag = np.sqrt(c[:, :, 0]**2 + c[:, :, 1]**2)
# 对剩余谱exp转换回原始空间，均值滤波
spectralResidual = np.exp(
np.log(mag) - cv2.boxFilter(np.log(mag), -1, (3, 3)))
# 转换到复数空间
c[:, :, 0] = c[:, :, 0] * spectralResidual / mag
c[:, :, 1] = c[:, :, 1] * spectralResidual / mag
# 反变换
c = cv2.dft(c, flags=(cv2.DFT_INVERSE | cv2.DFT_SCALE))
# 平方表示估计误差
mag = c[:, :, 0]**2 + c[:, :, 1]**2
# 高斯滤波后处理
cv2.normalize(cv2.GaussianBlur(mag, (9, 9), 3, 3),
mag, 0., 1., cv2.NORM_MINMAX)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(img, cmap='gray')
ax[1].imshow(mag, cmap='gray')
for i, j in enumerate(['raw', 'SaliencyMap']):
ax[i].set(title=j, xticks=[], yticks=[])
整体代码为（本段代码摘抄自opencv-python-blueprints/saliency.py at master · mbeyeler/opencv-python-blueprints (github.com)）：
import cv2
import numpy as np
from matplotlib import pyplot as plt
class Saliency:
"""Generate saliency map from RGB images with the spectral residual method
This class implements an algorithm that is based on the spectral
residual approach (Hou & Zhang, 2007).
"""
def __init__(self, img, use_numpy_fft=True, gauss_kernel=(5, 5)):
"""Constructor
This method initializes the saliency algorithm.
:param img: an RGB input image
:param use_numpy_fft: flag whether to use NumPy's FFT (True) or
OpenCV's FFT (False)
:param gauss_kernel: Kernel size for Gaussian blur
"""
self.use_numpy_fft = use_numpy_fft
self.gauss_kernel = gauss_kernel
self.frame_orig = img
# downsample image for processing
self.small_shape = (64, 64)
self.frame_small = cv2.resize(img, self.small_shape[1::-1])
# whether we need to do the math (True) or it has already
# been done (False)
self.need_saliency_map = True
def get_saliency_map(self):
"""Returns a saliency map
This method generates a saliency map for the image that was
passed to the class constructor.
:returns: grayscale saliency map
"""
if self.need_saliency_map:
# haven't calculated saliency map for this image yet
num_channels = 1
if len(self.frame_orig.shape) == 2:
# single channel
sal = self._get_channel_sal_magn(self.frame_small)
else:
# multiple channels: consider each channel independently
sal = np.zeros_like(self.frame_small).astype(np.float32)
for c in xrange(self.frame_small.shape[2]):
small = self.frame_small[:, :, c]
sal[:, :, c] = self._get_channel_sal_magn(small)
# overall saliency: channel mean
sal = np.mean(sal, 2)
# postprocess: blur, square, and normalize
if self.gauss_kernel is not None:
sal = cv2.GaussianBlur(sal, self.gauss_kernel, sigmaX=8,
sigmaY=0)
sal = sal**2
sal = np.float32(sal)/np.max(sal)
# scale up
sal = cv2.resize(sal, self.frame_orig.shape[1::-1])
# store a copy so we do the work only once per frame
self.saliencyMap = sal
self.need_saliency_map = False
return self.saliencyMap
def _get_channel_sal_magn(self, channel):
"""Returns the log-magnitude of the Fourier spectrum
This method calculates the log-magnitude of the Fourier spectrum
of a single-channel image. This image could be a regular grayscale
image, or a single color channel of an RGB image.
:param channel: single-channel input image
:returns: log-magnitude of Fourier spectrum
"""
# do FFT and get log-spectrum
if self.use_numpy_fft:
img_dft = np.fft.fft2(channel)
magnitude, angle = cv2.cartToPolar(np.real(img_dft),
np.imag(img_dft))
else:
img_dft = cv2.dft(np.float32(channel),
flags=cv2.DFT_COMPLEX_OUTPUT)
magnitude, angle = cv2.cartToPolar(img_dft[:, :, 0],
img_dft[:, :, 1])
# get log amplitude
log_ampl = np.log10(magnitude.clip(min=1e-9))
# blur log amplitude with avg filter
log_ampl_blur = cv2.blur(log_ampl, (3, 3))
# residual
residual = np.exp(log_ampl - log_ampl_blur)
# back to cartesian frequency domain
if self.use_numpy_fft:
real_part, imag_part = cv2.polarToCart(residual, angle)
img_combined = np.fft.ifft2(real_part + 1j*imag_part)
magnitude, _ = cv2.cartToPolar(np.real(img_combined),
np.imag(img_combined))
else:
img_dft[:, :, 0], img_dft[:, :, 1] = cv2.polarToCart(residual,
angle)
img_combined = cv2.idft(img_dft)
magnitude, _ = cv2.cartToPolar(img_combined[:, :, 0],
img_combined[:, :, 1])
return magnitude
def calc_magnitude_spectrum(self):
"""Plots the magnitude spectrum
This method calculates the magnitude spectrum of the image passed
to the class constructor.
:returns: magnitude spectrum
"""
# convert the frame to grayscale if necessary
if len(self.frame_orig.shape) > 2:
frame = cv2.cvtColor(self.frame_orig, cv2.COLOR_BGR2GRAY)
else:
frame = self.frame_orig
# expand the image to an optimal size for FFT
rows, cols = self.frame_orig.shape[:2]
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
frame = cv2.copyMakeBorder(frame, 0, ncols-cols, 0, nrows-rows,
cv2.BORDER_CONSTANT, value=0)
# do FFT and get log-spectrum
img_dft = np.fft.fft2(frame)
spectrum = np.log10(np.abs(np.fft.fftshift(img_dft)))
# return for plotting
return 255*spectrum/np.max(spectrum)
def plot_power_spectrum(self):
"""Plots the power spectrum
This method plots the power spectrum of the image passed to
the class constructor.
:returns: power spectrum
"""
# convert the frame to grayscale if necessary
if len(self.frame_orig.shape) > 2:
frame = cv2.cvtColor(self.frame_orig, cv2.COLOR_BGR2GRAY)
else:
frame = self.frame_orig
# expand the image to an optimal size for FFT
rows, cols = self.frame_orig.shape[:2]
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
frame = cv2.copyMakeBorder(frame, 0, ncols - cols, 0, nrows - rows,
cv2.BORDER_CONSTANT, value=0)
# do FFT and get log-spectrum
if self.use_numpy_fft:
img_dft = np.fft.fft2(frame)
spectrum = np.log10(np.real(np.abs(img_dft))**2)
else:
img_dft = cv2.dft(np.float32(frame), flags=cv2.DFT_COMPLEX_OUTPUT)
spectrum = np.log10(img_dft[:, :, 0]**2+img_dft[:, :, 1]**2)
# radial average
L = max(frame.shape)
freqs = np.fft.fftfreq(L)[:L//2]
dists = np.sqrt(np.fft.fftfreq(frame.shape[0])[:, np.newaxis]**2 +
np.fft.fftfreq(frame.shape[1])**2)
dcount = np.histogram(dists.ravel(), bins=freqs)[0]
histo, bins = np.histogram(dists.ravel(), bins=freqs,
weights=spectrum.ravel())
centers = (bins[:-1] + bins[1:]) / 2
plt.plot(centers, histo/dcount, marker='.')
plt.xlabel('frequency')
plt.ylabel('log-spectrum')
return plt.gca()
def get_proto_objects_map(self, use_otsu=True):
"""Returns the proto-objects map of an RGB image
This method generates a proto-objects map of an RGB image.
Proto-objects are saliency hot spots, generated by thresholding
the saliency map.
:param use_otsu: flag whether to use Otsu thresholding (True) or
a hardcoded threshold value (False)
:returns: proto-objects map
"""
saliency = self.get_saliency_map()
if use_otsu:
_, img_objects = cv2.threshold(np.uint8(saliency*255), 0, 255,
cv2.THRESH_BINARY + cv2.THRESH_OTSU)
else:
thresh = np.mean(saliency)*255*3
_, img_objects = cv2.threshold(np.uint8(saliency*255), thresh, 255,
cv2.THRESH_BINARY)
return img_objects
调用代码的实现：
s = Saliency(img)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(img, cmap='gray')
ax[1].imshow(s.get_saliency_map(), cmap='gray')
for i, j in enumerate(['raw', 'SaliencyMap']):
ax[i].set(title=j, xticks=[], yticks=[])
参考
- (PDF) Saliency Detection: A Spectral Residual Approach (researchgate.net)
- opencv-python-blueprints/saliency.py at master · mbeyeler/opencv-python-blueprints (github.com)
- SpectralResidualSaliency/saliency.py at master · uoip/SpectralResidualSaliency (github.com)
- OpenCV: blur函数的边界到底是怎么处理的?_Felaim的博客-CSDN博客
- Spectral Residual 小记 Grok j
- 基于Spectral Residual的时序异常检测 - 知乎 (zhihu.com)