本文介绍了各向异性扩散PM模型，并给出了C++代码实现。

一、PM模型原理

其中，                                                     

二、C++代码实现

MATLAB代码可参考：http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/anisodiff.m

http://www.mathworks.com/matlabcentral/fileexchange/14995-anisotropic-diffusion-perona-malik/content/anisodiff_Perona-Malik/anisodiff2D.m

void CImageObj::Perona_Malik(int iter, double dt, double kappa, int option)
{
	int i, j;
	int nx = m_width, ny = m_height;

	double** I_t = NewDoubleMatrix(nx, ny);
	double** I_tmp = NewDoubleMatrix(nx, ny);
	for (i = 0; i < ny; i++)
		for (j = 0; j < nx; j++)
			I_t[i][j] = I_tmp[i][j] = m_imgData[i][j];

	for (int t = 0; t < iter; t++)
	{
		for (i = 0; i < ny; i++)
		{
			for (j = 0; j < nx; j++)
			{
				int iUp = i - 1, iDown = i + 1;
				int jLeft = j - 1, jRight = j + 1;    // 边界处理
				if (0 == i) iUp = i; if (ny - 1 == i) iDown = i;
				if (0 == j) jLeft = j; if (nx - 1 == j) jRight = j;

				double deltaN = I_t[iUp][j] - I_t[i][j];
				double deltaS = I_t[iDown][j] - I_t[i][j];
				double deltaE = I_t[i][jRight] - I_t[i][j];
				double deltaW = I_t[i][jLeft] - I_t[i][j];

				double cN, cS, cE, cW;
				if (1 == option)
				{
					cN = exp(-(deltaN / kappa) * (deltaN / kappa));
					cS = exp(-(deltaS / kappa) * (deltaS / kappa));
					cE = exp(-(deltaE / kappa) * (deltaE / kappa));
					cW = exp(-(deltaW / kappa) * (deltaW / kappa));
				}
				else if (2 == option)
				{
					cN = 1.0 / (1 + (deltaN / kappa) * (deltaN / kappa));
					cS = 1.0 / (1 + (deltaS / kappa) * (deltaS / kappa));
					cE = 1.0 / (1 + (deltaE / kappa) * (deltaE / kappa));
					cW = 1.0 / (1 + (deltaW / kappa) * (deltaW / kappa));
				}

				I_tmp[i][j] += dt * (cN * deltaN + cS * deltaS + cE * deltaE + cW * deltaW);
			}
		}  // 一次迭代

		for (i = 0; i < ny; i++)
			for (j = 0; j < nx; j++)
			{
				I_t[i][j] = I_tmp[i][j];
			}

	} // 迭代结束

	// 给图像赋值
	for (i = 0; i < ny; i++)
		for (j = 0; j < nx; j++)
		{
			double tmp = I_t[i][j];
			tmp = max(0, min(tmp, 255));
			m_imgData[i][j] = (unsigned char)tmp;
		}

	DeleteDoubleMatrix(I_t, nx, ny);
	DeleteDoubleMatrix(I_tmp, nx, ny);
}