基于双边滤波（BF）的红外弱小目标检测之背景抑制Opencv
1. 前言
红外图像中的弱小目标，目标属性包涵“弱"和“小’’两个方面：“弱’’是指目标在红外波长上所表现出来的强度，具体反映到所拍摄的红外图像上，就是指目标的灰度值；而“小’’是指目标的尺寸大小，也就是前面所述的成像面积很小，反映到红外图像上就是指目标所占的像素数目很少。
SPIE国际光学工程学会(Societyof Photo一0ptical Instrumentation Engineers，简记为SPIE)从1989年开始，几乎每年都会举办有关弱小目标检测技术的国际会议，研讨弱小目标检测技术的最新研究成果。根据SPIE的定义，成像尺寸小于整个成像区域0.12％的目标均可称为弱小目标(即当成像尺寸为256×256时，弱小目标应不超过81个像素，其目标尺寸在9*9以内)。
2. 双边滤波算法
双边滤波的思想是抑制与中心像素差别太大的像素。高斯滤波是以距离为权重，设计滤波模板作为滤波系数，只考虑了像素间的空间位置上的关系，因此滤波的结果会丢失边缘的信息。高斯滤波的缺陷：平坦区域正常滤波，图像细节没有变化，而在突变的边缘上，因为只使用了距离来确定滤波权重，导致边缘被模糊。
在高斯基础上，进一步优化，叠加了像素值的考虑，因此也就引出了双边滤波，一种非线性滤波，滤波效果对保留边缘更有效。为了理解双边滤波的距离和像素差两个影响因素，先说明下面两个概念帮助理解。
（1）空间距离：当前点距离滤波模板中心点的欧式距离。
（2）灰度距离：当前点距离滤波模板中心点的灰度的差值的绝对值。
因此，双边滤波的核函数是空间域核与像素范围域核的综合结果：
（1）在图像的平坦区域，像素值变化很小，那么像素差值接近于0，对应的像素范围域权重接近于1，此时空间域权重起主要作用，相当于进行高斯模糊；
（2）在图像的边缘区域，像素值变化很大，那么像素差值大，对应的像素范围域权重变大，即使距离远空间域权重小，加上像素域权重总的系数也较大，从而保护了边缘的信息。
双边滤波在突变的边缘上，使用了像素差权重，所以很好的保留了边缘。
3. Opencv实例
需要创建VS工程，后续会补充介绍VS2015配置OPENCV环境，创建工程后，添加main.cpp文件即可生成，调试运行。
mian.cpp
/*
sigma_color参数是式中的σd, sigma_space参数是式中的σr
*/
#pragma once
#include "core/core.hpp" //3个opencv库的头文件
#include "highgui/highgui.hpp"
#include "imgproc/imgproc.hpp"
#include <iostream>
using namespace std;
using namespace cv;
static void bilateralFilter_8u(const Mat& src, Mat& dst, int d, double sigma_color, double sigma_space, int borderType)
{
int cn = src.channels();//获得图片的通道数
int i, j, k, maxk, radius;
Size size = src.size();
CV_Assert((src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
src.type() == dst.type() && src.size() == dst.size() &&
src.data != dst.data);
if (sigma_color <= 0)
sigma_color = 1;
if (sigma_space <= 0)
sigma_space = 1;
double gauss_color_coeff = -0.5 / (sigma_color*sigma_color);//分母值
double gauss_space_coeff = -0.5 / (sigma_space*sigma_space);//分母值
if (d <= 0)
radius = cvRound(sigma_space*1.5);//进行四舍五入
else
radius = d / 2;
radius = MAX(radius, 1);
d = radius * 2 + 1;
Mat temp;
copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
/*复制图像并且制作边界。（处理边界卷积）
目的是为了放大图像好做图像的边界处理*/
vector<float> _color_weight(cn * 256);//用来存放值域差值对应的权重
vector<float> _space_weight(d*d);//用来存放空间距离对应的权重
vector<int> _space_ofs(d*d);//用来存放模板各点与锚点（中心点）的偏移量
float* color_weight = &_color_weight[0];
float* space_weight = &_space_weight[0];
int* space_ofs = &_space_ofs[0];
// initialize color-related bilateral filter coefficients 函数1
//由于sigma_color已经给定，所以可以先算出差值为0-255时，对应的高斯相似度权重。
for (i = 0; i < 256 * cn; i++)
color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
// initialize space-related bilateral filter coefficients 函数2
//由于sigma_space已经给定，所以一旦选定好模板就可计算出高斯距离权重了
for (i = -radius, maxk = 0; i <= radius; i++)
for (j = -radius; j <= radius; j++)
{
double r = std::sqrt((double)i*i + (double)j*j);
/*if (r > radius)
continue;*/
space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
space_ofs[maxk++] = (int)(i*temp.step + j*cn);
}
for (i = 0; i < size.height; i++)
{
const uchar* sptr = temp.data + (i + radius)*temp.step + radius*cn;
uchar* dptr = dst.data + i*dst.step;
///灰度图
if (cn == 1)
{
for (j = 0; j < size.width; j++)
{
float sum = 0, wsum = 0;
int val0 = sptr[j];
for (k = 0; k < maxk; k++)
{
int val = sptr[j + space_ofs[k]];
float w = space_weight[k] * color_weight[std::abs(val - val0)];
sum += val*w;
wsum += w;
}
// overflow is not possible here => there is no need to use CV_CAST_8U
dptr[j] = (uchar)cvRound(sum / wsum);
}
}
else
{
//彩色图
assert(cn == 3);
for (j = 0; j < size.width * 3; j += 3)
{
float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
int b0 = sptr[j], g0 = sptr[j + 1], r0 = sptr[j + 2];
for (k = 0; k < maxk; k++)
{
const uchar* sptr_k = sptr + j + space_ofs[k];
int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
float w = space_weight[k] * color_weight[std::abs(b - b0) +
std::abs(g - g0) + std::abs(r - r0)];
sum_b += b*w; sum_g += g*w; sum_r += r*w;
wsum += w;
}
wsum = 1.f / wsum;
b0 = cvRound(sum_b*wsum);
g0 = cvRound(sum_g*wsum);
r0 = cvRound(sum_r*wsum);
dptr[j] = (uchar)b0; dptr[j + 1] = (uchar)g0; dptr[j + 2] = (uchar)r0;
}
}
}
}
//获得自适应值
double adaptive_gausscolor(float sum, float sumsqr, int n, double maxSigma_color)
{
double sigma_color = ((sumsqr*n) - sum*sum) / ((double)(n*n));
if (sigma_color < 0.01)sigma_color = 0.01;
else
if (sigma_color >(double)maxSigma_color)
sigma_color = (double)maxSigma_color;
double gauss_color_coeff = -0.5 / (double)(sigma_color*sigma_color);
return gauss_color_coeff;
}
自适应双边滤波算法
void adaptivebilateralFilter(Mat& src, Mat& dst, int d, double sigma_space, double sigmacolor_max, int borderType)
{
int cn = src.channels();//获得图片的通道数
int i, j, k, maxk, radius;
Size size = src.size();
CV_Assert((src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
src.type() == dst.type() && src.size() == dst.size() &&
src.data != dst.data);
if (sigma_space <= 0)
sigma_space = 1;
CV_Assert(d & 1);//确保为奇数
double gauss_space_coeff = -0.5 / (sigma_space*sigma_space);
if (d <= 0)
radius = cvRound(sigma_space*1.5);//进行四舍五入
else
radius = d / 2;
radius = MAX(radius, 1);
d = radius * 2 + 1;
Mat temp;
copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
/*复制图像并且制作边界。（处理边界卷积）
目的是为了放大图像好做图像的边界处理*/
vector<float> _space_weight(d*d);
vector<int> _space_ofs(d*d);
float* space_weight = &_space_weight[0];
int* space_ofs = &_space_ofs[0];
// initialize color-related bilateral filter coefficients 函数1
//由于sigma_color没有给定，所以无法算出差值为0-255时，对应的高斯相似度权重。
/*for (i = 0; i < 256 * cn; i++)
color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);*/
// initialize space-related bilateral filter coefficients
//由于sigma_space已经给定，所以一旦选定好模板就可计算出高斯距离权重了
for (i = -radius, maxk = 0; i <= radius; i++)
for (j = -radius; j <= radius; j++)
{
double r = std::sqrt((double)i*i + (double)j*j);
space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
space_ofs[maxk++] = (int)(i*temp.step + j*cn);
}
for (i = 0; i < size.height; i++)
{
const uchar* sptr = temp.data + (i + radius)*temp.step + radius*cn;
uchar* dptr = dst.data + i*dst.step;
if (cn == 1)
{
for (j = 0; j < size.width; j++)
{
float sum = 0, wsum = 0, sumValsqr = 0;
int val0 = sptr[j];
//获得自适应的高斯相似度的sigma值并计算相似度
for (k = 0; k < maxk; k++)
{
int val = sptr[j + space_ofs[k]];
sum += val;
sumValsqr += (val*val);
}
double gauss_color_coeff = adaptive_gausscolor(sum, sumValsqr, maxk, sigmacolor_max);
sum = 0;
for (k = 0; k < maxk; k++)
{
int val = sptr[j + space_ofs[k]];
int temp = std::abs(val - val0);
float color_Weight = (float)std::exp((float)temp*temp*gauss_color_coeff);
float w = space_weight[k] * color_Weight;
sum += val*w;
wsum += w;
}
// overflow is not possible here => there is no need to use CV_CAST_8U
dptr[j] = (uchar)cvRound(sum / wsum);
}
}
else
{
assert(cn == 3);
for (j = 0; j < size.width * 3; j += 3)
{
float sum_b = 0, sum_g = 0, sum_r = 0, wbsum = 0, wgsum = 0, wrsum = 0, sum_bsqr = 0, sum_gsqr = 0, sum_rsqr = 0;
int b0 = sptr[j], g0 = sptr[j + 1], r0 = sptr[j + 2];
//获得自适应的高斯相似度的sigma值并计算相似度
for (k = 0; k < maxk; k++)
{
const uchar* sptr_k = sptr + j + space_ofs[k];
int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
sum_b += b; sum_g += g; sum_r += r;
sum_bsqr += b*b; sum_gsqr += g*g; sum_rsqr += r*r;
}
double gauss_colorb_coeff = adaptive_gausscolor(sum_b, sum_bsqr, maxk, sigmacolor_max);
double gauss_colorg_coeff = adaptive_gausscolor(sum_g, sum_gsqr, maxk, sigmacolor_max);
double gauss_colorr_coeff = adaptive_gausscolor(sum_r, sum_rsqr, maxk, sigmacolor_max);
sum_b = 0; sum_g = 0; sum_r = 0;
for (k = 0; k < maxk; k++)
{
const uchar* sptr_k = sptr + j + space_ofs[k];
int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
double colorb_Weight = (double)std::exp((std::abs(b - b0))*(std::abs(b - b0))*gauss_colorb_coeff);
double colorg_Weight = (double)std::exp((std::abs(g - g0))*(std::abs(g - g0))*gauss_colorg_coeff);
double colorr_Weight = (double)std::exp((std::abs(r - r0))*(std::abs(r - r0))*gauss_colorr_coeff);
float wb = space_weight[k] * colorb_Weight;
float wg = space_weight[k] * colorg_Weight;
float wr = space_weight[k] * colorr_Weight;
sum_b += b*wb; sum_g += g*wg; sum_r += r*wr;
wbsum += wb;
wgsum += wg;
wrsum += wr;
}
wbsum = 1.f / wbsum;
wgsum = 1.f / wgsum;
wrsum = 1.f / wrsum;
b0 = cvRound(sum_b*wbsum);
g0 = cvRound(sum_g*wgsum);
r0 = cvRound(sum_r*wrsum);
dptr[j] = (uchar)b0; dptr[j + 1] = (uchar)g0; dptr[j + 2] = (uchar)r0;
}
}
}
}
int main()
{
Mat srcImage, dstImageBF, dstImageABF, ImageBF, ImageABF; //源图像，输出图像
//---------【1】读取源图像并检查图像是否读取成功---------
srcImage = imread("./src\\69.bmp", 0);
if (!srcImage.data)
{
cout << "读取图片错误，请重新输入正确路径！\n";
system("pause");
return-1;
}
int type = srcImage.type();
dstImageBF.create(srcImage.rows, srcImage.cols, type);
imshow("【源图像】", srcImage);
imwrite("./dst\\srcImage.png", srcImage);
//---------【2】进行双边滤波操作---------
bilateralFilter_8u(srcImage, dstImageBF,9,150,150, BORDER_REFLECT);
//---------【3】显示效果图---------
imshow("【双边滤波】", dstImageBF);
//---------【4】TopHAT变换，并显示效果图---------
ImageBF = srcImage - dstImageBF;
imshow("【双边滤波抑制背景】", ImageBF);
//---------【5】进行自适应双边滤波操作---------
dstImageABF.create(srcImage.rows, srcImage.cols, type);
adaptivebilateralFilter(srcImage, dstImageABF,9,150, 200, BORDER_REFLECT);
//---------【6】进行自适应双边滤波操作---------
imshow("【自适应双边滤波】", dstImageABF);
//---------【7】显示自适应双边滤波抑制效果图---------
ImageABF = srcImage - dstImageABF;
imshow("【自适应双边滤波抑制】", ImageABF);
imwrite("./dst\\dstImageBF.png", dstImageBF);
imwrite("./dst\\ImageBF.png", ImageBF);
imwrite("./dst\\dstImageABF.png", dstImageABF);
imwrite("./dst\\ImageABF.png", ImageABF);
waitKey(0);
return 0;
}
4. 小结
双边滤波是一种非线性滤波器，它可以达到保持边缘、降噪平滑的效果。在视网膜、血管检测中应用较多，本次应用在复杂背景下弱小目标的背景抑制中，后续试经过优化，选出综合效果最好的进行详细介绍。
本系列文章列表如下：
弱小目标检测跟踪算法研究（1）红外弱小目标数据集准备（红外弱小目标的数据集来了）
弱小目标检测跟踪算法研究（2）红外弱小目标数据集准备（红外弱小目标的数据集又来了）
弱小目标检测跟踪算法研究（3） 基于局部对比度（LCM）的红外弱小目标检测之背景抑制
弱小目标检测跟踪算法研究（4） 基于双边滤波（BF）的红外弱小目标检测之背景抑制
弱小目标检测跟踪算法研究（5） 基于顶帽变换（Top_hat）算法的红外弱小目标检测之背景抑制
弱小目标检测跟踪算法研究（6） 基于小波变换的红外弱小目标检测之背景抑制
弱小目标检测跟踪算法研究（7） 基于简单平滑滤波算法的红外弱小目标检测之背景抑制
