基于各向异性扩散的数字图像处理 作业

基于各向异性扩散的数字图像处理
1 Perona-Malik 算法概述
各向异性扩散作为现行的一种非常流行的偏微分方程数字图像处理技术，是由传统的 Gaussian 滤波发展而来的，有着强大的理论基础，并有着传统的数字图像方法无法企及的良好特性，其特点是可以在平滑的同时保持边缘特征。

由于这种优良的特性，使其在图像的平滑、去噪、恢复、增强和分割等方面得到了广泛的应用。

Koenderink 和 Witkin 两位学者把尺度空间的严格理论引入到了数字图像处理之中，而尺度空间理论正是现代偏微分方程数字图像处理的理论基础，他们的工作主要是将多尺度图像表示为 Gaussian 滤波器处理的结果，相当于将原图像输入热传导方程进而得到的序列时间图像。

多尺度滤波得到的图像序列，可以看作各向同性热扩散方程的解，这就是用扩散方程的观点来看待高斯滤波：
()
yy xx u u c t u
-=∂∂ (1.1)
不仅是热扩散方程可以产生尺度空间，其他的抛物线方程也可以产生尺度空间，进而满足极大值原理的演化方程也能定义一类尺度空间。

各向异性扩散方程的提出是这个领域中一个里程碑式的标志，它开辟了数字图像处理中偏微分方程理论和应用的一个新领域，提出了各向异性扩散（Anisotropic Diffusion ）方程，一般称为 Perona-Malik （简称 P-M ）扩散模型。

)
)((u u g div t u
∇∇=∂∂
(1.2) 其中div 是散度算子，u ∇是图像的梯度，)(u g ∇是扩散系数。

各向异性扩散是一个能量散发的过程，这个过程与能量曲面的形状有关，从数学上来看，其等价于一个能量最小化问题的求解。

从数学意义上讲Perona-Malik 模型是改进的热传导偏微分方程。

它是通过函数(,,)c x y t 自适应地控制扩散速度。

理论上希望边缘内部的区域，c 远离0，图像可以平滑，而在边缘附近，0c ≈，图像不再平滑。

铃铛形的径向函数g(w)可取
2
1()exp[(
)]u K
g u ∇∇=- (1.3)
或者
211(),01(
)
u K
g u αα∇+∇=
>+ (1.4)
当(,,)()c x y t g u =∇时形成了各向异性的扩散模型：
()·)()()?u
g u u g u u g u u t
∂=∇∇∇=∇∆+∇∇∇∂ (1.5) 可以利用Taylor 级数展开进行线性近似
(,,)(,,0)(,,0)((,,))u
u x y t u x y t u x y t c x y t u t
∂=+∆=+∆∇⋅∇∂ (1.6)
由于
()··
c u c u c u ∇∇=∇∇+∆ (1.7) 所以可以得出P-K 模型简化推导：
0(,,)(,,0)()(,,0)((,,0))t u
u x y t u x y t u x y t c x y u t δδδ=∂≈+=+∇⋅∇∂ (1.8)
这样把较长的时间t 分割为 t n t δ=⋅，始终从n t 步计算到1n t +步即可
112(,,)(,,)((,,)(,,))
(,,)[(,,)(,,)(,,)(,,)](,,)[]
n n n n n n n n n n n
n u x y t u x y t t c x y t u x y t u x y t t c x y t u x y t c x y t u x y t u x y t t I I δδδ+≈+∇⋅∇=+∇∇+∆=++ (1.9)
把(,,)n u x y t 改写为,n i j u ，(,,)n c x y t 改写为,n i j c ，其中
11,.1,,,1.,1,, 1.,1,,.1,,11
[()()()()2()()()()]
n n n n n n n n n i j i j i j i j i j i j i j i j n n n n n n n n i j i j i j i j i j i j i j i j I c c u u c c u u c c u u c c u u ++++----=--+--+
--+-- (1.10) 2,1,1,,1,1,1[4]2
n
n n n n n n i j i j i j i j i j i j I c u u u u u +-+-=+++-
(1.11)
则
121,,,1,1,,,1,1[]2
n n
n
n n n n n n n i j S i j i j E i j i j N i j i j W i j I I c u c u c u c u ++--+=
∇+∇+∇+∇ (1.12) 其中
,1,,,,1,,1,,,,1,n n n S i j i j i j
n n n E i j i j i j n n n N i j
i j
i j
n n n W i j i j i j
u u u u u u u u
u
u u u ++--∇=-∇=-∇=-∇=- (1.13)
则得到Perona-Malik 算法的下述迭代形式
1
,,1,,,1,1,,,1,[]n n
n
n
n
n
n
n
n
n
i j i j i j S i j i i E i j i j N i j i j W i j u u c u c u c u c u λ+++--=+∇+∇+∇+∇ (1.14)
其中，2
t
δλ=
以Perona-Malik 算法为代表的各向异性扩散滤波模型既能有效地去除噪声，又能较好地保持图像的边缘和重要的细节信息。

但是考虑到PM 模型对小尺度噪声敏感以及部分边缘和细节失真问题的存在，对模型的扩散系数的改进仍然需要进一步的改进。

2算法实现
本算法的实验在VC++6.0环境下实现。

先建一个类CDib ，用来实现位图的载入以及一些参数的获取。

新建一个类CPeronaMalik ，类中函数PM( CDib* pDib )用来实现PM 算法。

void CPeronaMalik:: PM( CDib* pDib ) ｛
int i=0,j=0; int nWidth = pDib->GetWidth(); int nHeight = pDib->GetHeight(); int nWidthBytes = pDib->GetWidth();
unsigned char* pDibBits = pDib->GetpDibBits(); //每行字节数
int lineByte=(bmpWidth*biBitCount/8+3)/4*4;
//循环变量，针对彩色图像，遍历每像素的三个分量 int k; int n;
double temp1=0; double temp2=0; int nn;//迭代次数 double deltat;
double delta;//步长，deltat 小于0.25 //c1、c2表示扩散函数 int p;//范数
double sigma;//平均灰度
double K;
if(biBitCount==8)
{//对于灰度图像
for(i=0;i<nHeight;i++)
{
for(j=0;j<nWidth;j++)
{
temp1=double(pDibBits[i*nWidthBytes+j])+temp1;
}
}
//添加图像处理的代码
////////////定参数////////////////////////
nn=200;//迭代次数
deltat=0.02;delta=deltat/2;//步长，deltat小于0.25
//c1、c2表示扩散函数
p=2;//范数
sigma=0.4*temp1/(double(bmpHeight)*double(bmpWidth));//平均灰度int a=1;//参数,a>0
//////////////////////////////////////////
for(i=0;i<nHeight;i++)
{
for(j=0;j<bmpWidth;j++)
{
I1[i][j]=double(*(pBmpBuf+i*lineByte+j));
I2[i][j]=double(*(pBmpBuf+i*lineByte+j));//交替存储}
}
for(n=0;n<nn;n++)
{
temp2=0;
K=1;
for(i=0;i<bmpHeight;i++)
{
for(j=0;j<bmpWidth;j++)
{
temp2=pow(I1[i][j],p)+temp2;
}
}
K=sigma*pow(temp2,1/double(p))/(bmpHeight*bmpWidth);
////////////////
for(i=0;i<nHeight;i++)
{
for(j=0;j<nWidth;j++)
{
if(i==0||i==bmpHeight-1||j==0||j==bmpWidth-1)
RF[i][j]=0;
else
RF[i][j]=1/(1+pow((pow(pow((I1[i+1][j]-I1[i-1][j])/2,2)+pow((I1[i][j+1]-I1[i][j-1])/2,2),0.
5)/K),1+a));//扩散函数
}
}
////////////////
for(i=1;i<(nHeight-1);i++)
{
for(j=1;j<(nWidth-1);j++)
I2[i][j]=I1[i][j]+delta*(RF[i+1][j]*(I1[i+1][j]-I1[i][j])+RF[i-1][j]*(I1[i-1][j]-I1[i][j])+RF[i][j+ 1]*(I1[i][j+1]-I1[i][j])+RF[i][j-1]*(I1[i][j-1]-I1[i][j]));
//迭代
}
for(i=0;i<bmpHeight;i++)
{
for(j=0;j<bmpWidth;j++)
{
I1[i][j]=I2[i][j];
}
}
}
for(i=0;i<nHeight;i++)
{
for(j=0;j<nWidth;j++)
{
if(I1[i][j]>255)
*(pBmpBuf+i*lineByte+j)=255;
else
if(I1[i][j]<0)
*(pBmpBuf+i*lineByte+j)=0;
else
*(pBmpBuf+i*lineByte+j)=unsigned char(int(I1[i][j]));
//cout<<int(*(pBmpBuf+i*lineByte+j))<<endl;
}
}
}
else if(biBitCount==24)
{//彩色图像
for(i=0;i<nHeight/2;i++){
for(j=0;j<nWidth/2;j++){
for(k=0;k<3;k++)//每像素RGB三个分量分别置0才变成黑色
*(pBmpBuf+i*lineByte+j*3+k)=0;
}
｝
3 实验结果
源图采用lena256*256 灰度图像，经过处理后如右图
可以看出，平滑效果很好，边缘保持的效果也是不错的。

对图像加椒盐噪声及处理后，
效果不是很好，噪声没能去除。