一种简单的图像显著性计算模型

《Saliency Detection: A Spectral Residual Approach》是上交高材生侯晓迪在07年的CVPR上发表的一篇论文。这篇文章提出了一个图像视觉显著性的简单计算模型。有关算法原理可以参见此篇论文以及作者主页：http://www.klab.caltech.edu/~xhou/projects/spectralResidual/spectralresidual.html

这篇文章的核心思想如下：

从信息理论角度：信息可分为冗余部分和变化部分。人们的视觉对变化部分更敏感。视觉系统的一个基本原则就是抑制对频繁出现的特征的响应，同时对非常规的特征保持敏感。那么就将图像分为如下两部分：

 而作者对图像的log频谱发下了如下规律（log是自然对数）:

 大量图像的log频谱的平均值是和频率呈现正比关系的。

如图可知大量图像的log频谱和频率的曲线形状，在log-log scale上，几乎是一条直线。文中的log频谱就是对图像傅里叶变换后的振幅谱取自然对数。然后作者又提出了既然大量图像的log振幅谱都差不多趋近一条直线，那么一幅图像的log振幅谱减去平均log振幅谱不就是显著性部分了吗？这就是作者提出的：Spectral Residual理论。

 作者如下的定义Spectral Residual：

 其中L(f)就是图像f的log振幅谱，A(f)是平均log振幅谱，作者就是将L(f)进行3*3的均值滤波得到的A(f)。R(f)就是图像f的Spectral Residual。计算过程如下表示：

其中输入图像为I(x)，对其傅里叶变换，并且求出振幅谱为A(f)。P(f)是其相位谱（复数x+i*y的相位是arctan(y/x)）。L(f)是log振幅谱。h是一个n*n均值滤波的卷积核，作者设n=3。R(f)就是Spectral Residual谱。再将R(f)+i*P(f)求出自然指数exp。注意：由欧拉公式可知，exp(r+i*Θ) = exp(r)*(cos(Θ) + i*sin(Θ)) = exp(r)*cos(Θ) + i*exp(r)*sin(Θ)，Θ是相位谱。其实Θ的正余弦也可由傅里叶谱的实部和虚部求出：sin(Θ) = ImageIm/Amplitude; cos(Θ) = ImageRe/Amplitude。然后对其，傅里叶反变换，在进行一个高斯模糊滤波就得到了所谓的显著性区域。效果如下：

作者给出的matlab代码只有5行：

clear
clc

%% Read image from file
inImg = im2double(rgb2gray(imread('256.png')));
%%inImg = imresize(inImg, 64/size(inImg, 2));

%% Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySpectralResidual = myLogAmplitude - imfilter(myLogAmplitude, fspecial('average', 3), 'replicate');
saliencyMap = abs(ifft2(exp(mySpectralResidual + i*myPhase))).^2;

%% After Effect
saliencyMap = mat2gray(imfilter(saliencyMap, fspecial('gaussian', [10, 10], 2.5)));
imshow(saliencyMap);

我实现的基于OpenCV的C++代码：

/***********************************************************************
 * OpenCV 2.1 测试例程
 ***********************************************************************/
#include "highgui.h"
#include "cv.h"
#include <string>
 
//所有的以新风格命名的函数都在 cv 命名空间中
//如果希望不要每次都输入 cv:: ，则可使用下面语句
using namespace cv;
using namespace std;

//傅里叶正变换
void fft2(IplImage *src, IplImage *dst);

int main()
{
    string imagePath = "41.jpg";

    //////////////////////////////////////////////////////////////////////////
    //显著性计算
    //参考论文：Saliency Detection: A Spectral Residual Approach
    
    //amplitude和phase分别是图像的振幅谱和相位谱
    IplImage *src, *ImageRe, *ImageIm, *Fourier, *Inverse, *LogAmplitude, *Sine, *Cosine;
    IplImage *Saliency, *Residual;
    IplImage *tmp1, *tmp2, *tmp3;
    double minNum = 0, maxNum = 0, scale, shift;
    int i, j,     nRow, nCol;
    
    //加载源图像，第二个参数为0表示将输入的图片转为单通道，大于0为三通道
    src = cvLoadImage(imagePath.c_str(),0);
    //注意Fourier是一个两通道的图像，一个通道为实部，一个为虚部
    Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
    Inverse = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
    //频谱的实部
    ImageRe = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    //频谱的虚部
    ImageIm = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    //log振幅谱
    LogAmplitude = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    //正弦谱
    Sine = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    //余弦谱
    Cosine = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);

    //频谱冗余（spectral residual）
    Residual = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    //特征map(Saliency map)
    Saliency = cvCreateImage(cvGetSize(src), src->depth, src->nChannels);

    //临时的空间
    tmp1 = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    tmp2 = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
    tmp3 = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);

    nRow = src->height;
    nCol = src->width;

    //归一化一下
    scale = 1.0/255.0;
    cvConvertScale(src, tmp1, scale, 0);

    //傅里叶变换，得到的Fourier有两个通道，一个是实部，一个是虚部
    fft2(tmp1, Fourier);

    //将傅里叶谱的实部和虚部存放到对应的图像中去。
    cvSplit(Fourier, ImageRe, ImageIm, 0, 0);

    //计算傅里叶振幅谱，实部和虚部平方和再开方，得到振幅谱存到tmp3中
    cvPow( ImageRe, tmp1, 2.0);
    cvPow( ImageIm, tmp2, 2.0);
    cvAdd( tmp1, tmp2, tmp3);
    cvPow( tmp3, tmp3, 0.5 );

    //计算正弦谱和余弦谱和自然对数谱
    cvLog( tmp3, LogAmplitude );
    cvDiv(ImageIm, tmp3, Sine);
    cvDiv(ImageRe, tmp3, Cosine);

    //对LogAmplitude做3*3均值滤波
    cvSmooth(LogAmplitude, tmp3, CV_BLUR, 3, 3);

    //计算出剩余普
    cvSub(LogAmplitude, tmp3, Residual);

    /************************************************************************/
    /*
    对exp(Residual+i*Phase)作傅里叶反变换
    由欧拉公式：
    exp(r+i*Θ) = exp(r)*(cos(Θ) + i*sin(Θ)) = exp(r)*cos(Θ) + i*exp(r)*sin(Θ)
    其中Θ是相位。

    而sin(Θ) = ImageIm/Amplitude; cos(Θ) = ImageRe/Amplitude;
    */
    /************************************************************************/
    cvExp(Residual, Residual);
    cvMul(Residual, Cosine, tmp1);
    cvMul(Residual, Sine, tmp2);

    //将剩余普Residual作为实部，相位谱Phase作为虚部
    cvMerge(tmp1, tmp2, 0, 0, Fourier);
    
    //实现傅里叶逆变换
    cvDFT(Fourier, Inverse, CV_DXT_INV_SCALE);

    cvSplit(Inverse ,tmp1, tmp2,0,0);

    //求出对应的实部虚部平方和
    cvPow(tmp1, tmp1, 2);
    cvPow(tmp2, tmp2, 2);
    cvAdd(tmp1, tmp2, tmp3, NULL);

    //高斯滤波
    cvSmooth(tmp3, tmp3, CV_GAUSSIAN, 7, 7);

    cvMinMaxLoc(tmp3, &minNum, &maxNum, NULL, NULL);
    scale = 255/(maxNum - minNum);
    shift = -minNum * scale;

    //将shift加在ImageRe各元素按比例缩放的结果上，存储为ImageDst
    cvConvertScale(tmp3, Saliency, scale, shift);

    cvNamedWindow("Saliency", 1);
    cvShowImage("Saliency",Saliency);

    cvWaitKey(0);

    //释放图像
    cvReleaseImage(&src);
    cvReleaseImage(&ImageIm);
    cvReleaseImage(&ImageRe);
    cvReleaseImage(&Fourier);
    cvReleaseImage(&Inverse);
    cvReleaseImage(&LogAmplitude);
    cvReleaseImage(&Sine);
    cvReleaseImage(&Cosine);
    cvReleaseImage(&Saliency);
    cvReleaseImage(&Residual);
    cvReleaseImage(&tmp1);
    cvReleaseImage(&tmp2);
    cvReleaseImage(&tmp3);

    cvDestroyAllWindows();

    return 0;
}

/**************************************************************************
//src IPL_DEPTH_8U
//dst IPL_DEPTH_64F
**************************************************************************/
//傅里叶正变换
void fft2(IplImage *src, IplImage *dst)
{   //实部、虚部
    IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;
    //   int i, j;
    image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //实部
    //Imaginary part
    image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虚部
    //2 channels (image_Re, image_Im)
    Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);

    /************************************************************************/
    /*
    void cvConvertScale( const CvArr* src, CvArr* dst, double scale=1, double shift=0 );
    src
    原数组. 
    dst
    输出数组 
    scale
    比例因子. 
    shift
    原数组元素按比例缩放后添加的值。 

    函数 cvConvertScale 有多个不同的目的因此就有多个意义，
    函数按比例从一个数组中拷贝元素到另一个元素这种操作是最先执行的，
    或者任意的类型转换，正如下面的操作：

    dst(I)=src(I)*scale + (shift,shift,...) 

    多通道的数组对各个地区通道是独立处理的。
    */
    /************************************************************************/

    //实部的值初始设为源图像，虚部的值初始设为0
    // Real part conversion from u8 to 64f (double)
    cvConvertScale(src, image_Re, 1, 0);
    // Imaginary part (zeros)
    cvZero(image_Im);
    // Join real and imaginary parts and stock them in Fourier image
    cvMerge(image_Re, image_Im, 0, 0, Fourier);
    // Application of the forward Fourier transform
    cvDFT(Fourier, dst, CV_DXT_FORWARD);
    cvReleaseImage(&image_Re);
    cvReleaseImage(&image_Im);
    cvReleaseImage(&Fourier);
}

效果如下，分别是输入图像和显著性区域标识：