各向异性扩散滤波(Anisotropic Filter)原理与C++实现
各向异性扩散滤波主要是用来平滑图像,克服了高斯滤波模糊的缺陷,各向异性扩散在平滑图像的同时又能保留图像边缘。
各向异性扩散滤波(Anisotropic Filter,又称 Perona-Malik Filter),主要是用来平滑图像,克服了高斯滤波模糊的缺陷,各向异性扩散在平滑图像的同时又能保留图像边缘;此外,Perona-Malik Filter 在地学影像处理中也得到了相应的应用,如平滑DEM数据、河道提取 等。
将图像看作热量场,每个像元看作热流,根据当前像元和周围像元的关系,来确定是否需要向周围扩散。如果某个邻域像元和当前像元差别较大,那么这个邻域像元可能是边界,当前像元就不向这个方向扩散了,该边界也就得到保留。
主要迭代方程如下:
I t + 1 = I t + λ ( c N x , y ∇ N ( I t ) + c S x , y ∇ S ( I t ) + c E x , y ∇ E ( I t ) + c W x , y ∇ W ( I t ) ) I_{t+1}=I_t+\lambda(cN_{x,y}\nabla_N(I_t)+cS_{x,y}\nabla_S(I_t)+cE_{x,y}\nabla_E(I_t)+cW_{x,y}\nabla_W(I_t)) It+1=It+λ(cNx,y∇N(It)+cSx,y∇S(It)+cEx,y∇E(It)+cWx,y∇W(It))其中, I I I 为图像, t t t 为迭代次数。
四个散度公式是在四个方向上对当前像素求偏导,公式如下:
∇ N ( I x , y ) = I x , y − 1 − I x , y ∇ S ( I x , y ) = I x , y + 1 − I x , y ∇ E ( I x , y ) = I x − 1 , y − I x , y ∇ W ( I x , y ) = I x + 1 , y − I x , y \begin{array}{lcl} \nabla_N(I_{x,y})=I_{x,y-1}-I_{x,y} \\ \nabla_S(I_{x,y})=I_{x,y+1}-I_{x,y} \\ \nabla_E(I_{x,y})=I_{x-1,y}-I_{x,y} \\ \nabla_W(I_{x,y})=I_{x+1,y}-I_{x,y} \end{array} ∇N(Ix,y)=Ix,y−1−Ix,y∇S(Ix,y)=Ix,y+1−Ix,y∇E(Ix,y)=Ix−1,y−Ix,y∇W(Ix,y)=Ix+1,y−Ix,y
而 c N cN cN、 c S cS cS、 c E cE cE 和 c W cW cW 代表四个方向上的导热系数,公式如下:
c N x , y = e − ∥ ∇ N ( I ) ∥ 2 k 2 c S x , y = e − ∥ ∇ S ( I ) ∥ 2 k 2 c E x , y = e − ∥ ∇ E ( I ) ∥ 2 k 2 c W x , y = e − ∥ ∇ W ( I ) ∥ 2 k 2 \begin{array}{lcl} cN_{x,y}=e^{-\frac{\rVert\nabla_N(I)\rVert^2}{k^2}} \\ cS_{x,y}=e^{-\frac{\rVert\nabla_S(I)\rVert^2}{k^2}} \\ cE_{x,y}=e^{-\frac{\rVert\nabla_E(I)\rVert^2}{k^2}} \\ cW_{x,y}=e^{-\frac{\rVert\nabla_W(I)\rVert^2}{k^2}} \end{array} cNx,y=e−k2∥∇N(I)∥2cSx,y=e−k2∥∇S(I)∥2cEx,y=e−k2∥∇E(I)∥2cWx,y=e−k2∥∇W(I)∥2
整个公式需要设置的参数主要有三个,迭代次数 t t t,根据情况设置;导热系数相关的 k k k,取值越大,结果越平滑,越不易保留边缘; λ \lambda λ 取值越大,结果越平滑。
C++实现(未使用opencv):
void AnisotropicFilter(float *ppafScan, float *outppafScan, int nImgWidth, int nImgHeight, int iterationSum, float k, float lambda)
{
//ppafScan 输入数据
//outppafScan 输出数据
//nImgWidth, nImgHeight 数据的宽和高
//iterationSum 迭代总次数
//k 导热系数,控制平滑
//lambda 控制平滑
for(int iter = 0; iter < iterationSum; ++iter)
{
for(int i = 0; i < nImgHeight; ++i)
{
for (int j = 0; j < nImgWidth; ++j)
{
if(i == 0 || i == nImgHeight - 1 || j == 0 || j == nImgWidth - 1)
{
outppafScan[i * nImgWidth + j] = ppafScan[i * nImgWidth + j];
continue;
}
float NI = ppafScan[i * nImgWidth + (j - 1)] - ppafScan[i * nImgWidth + j];
float EI = ppafScan[(i - 1) * nImgWidth + j] - ppafScan[i * nImgWidth + j];
float WI = ppafScan[(i + 1) * nImgWidth + j] - ppafScan[i * nImgWidth + j];
float SI = ppafScan[i * nImgWidth + (j + 1)] - ppafScan[i * nImgWidth + j];
float cN = exp(-NI * NI / (k * k));
float cE = exp(-EI * EI / (k * k));
float cW = exp(-WI * WI / (k * k));
float cS = exp(-SI * SI / (k * k));
outppafScan[i * nImgWidth + j] = ppafScan[i * nImgWidth + j] + lambda * (cN * NI + cS * SI + cE * EI + cW * WI);
}
}
}
return;
}
当 t = 20 t=20 t=20, k = 30 k=30 k=30, λ = 0.20 \lambda=0.20 λ=0.20 时,结果如下图所示:
左:原始图像,右:滤波后的影像
高斯滤波(线性)和Perona-Malik滤波(非线性)的比较:(a)原始影像;(b) Gaussian filtering (the kernels = 7 m);© Perona-Malik filter (t = 50);(d) Gaussian filtering (the kernels = 14 m);(e) Perona-Malik filter (t = 200) [原图链接]
欢迎大家批评指正。
更多推荐
所有评论(0)