PM方程源码

cpp 

void Perona_Malik(
	unsigned char* inbuffer,
	int iter,
	double delta_t,
	double kappa,
	int option,
	int width,
	int height,
	unsigned char* outbuffer
	)
{
	double dx = 1.0;
	double dy = 1.0;
	double dd = sqrt(2.0);
	//变换数据类型，准备处理
	double* tmpbuffer1 = new double[width*height];
	double* tmpbuffer2 = new double[width*height];
	for (int kk = 0; kk < width*height; kk++)
	{
		tmpbuffer1[kk] = double(inbuffer[kk]);
		tmpbuffer2[kk] = double(inbuffer[kk]);
	}
	int i, j;
	double cN = 0.0, cS = 0.0, cE = 0.0, cW = 0.0;
	double cNE = 0.0, cSE = 0.0, cSW = 0.0, cNW = 0.0;
	double nablaN = 0.0, nablaS = 0.0, nablaE = 0.0, nablaW = 0.0;
	double nablaNE = 0.0, nablaSE = 0.0, nablaNW = 0.0, nablaSW = 0.0;
	for (int t = 0; t < iter; t++)//迭代次数
	{
		for (i = 0; i < height; i++)
		{
			for (j = 0; j < width; j++)//列
			{
				double pixel = -tmpbuffer1[i*height + j];
				nablaS = nablaE = nablaW = nablaN = pixel;
				nablaNE = nablaSE = nablaNW = nablaSW = pixel;
//八个方向上的求导
				if (i + 1 < height)//不是最后一行
					nablaN += tmpbuffer1[(i + 1)*height + j];
 
				if (i > 0)//不是第一行
					nablaS += tmpbuffer1[(i - 1)*height + j];
 
				if (j + 1 < width)//不是最后一列
					nablaE += tmpbuffer1[i*height + j + 1];
			
				if (j > 0)//不是第一列
					nablaW += tmpbuffer1[i*height + j - 1];
 
				if (i + 1 < height && j > 0)
					nablaNE += tmpbuffer1[(i + 1)*height + j -1 ];
 
				if (i  > 0 && j > 0)
					nablaSE += tmpbuffer1[(i - 1)*height + j - 1];
				
				if (i > 0 && j + 1 < width)
					nablaNW += tmpbuffer1[(i - 1)*height + j + 1];
 
				if (i + 1 < height  && j + 1 < width)
					nablaSW += tmpbuffer1[(i + 1)*height + j + 1];
 
				cN = 0.0, cS = 0.0, cE = 0.0, cW = 0.0;
				if (1 == option)//扩散函数的选择
				{
					cN = exp(-(nablaN / kappa) * (nablaN / kappa));
					cS = exp(-(nablaS / kappa) * (nablaS / kappa));
					cE = exp(-(nablaE / kappa) * (nablaE / kappa));
					cW = exp(-(nablaW / kappa) * (nablaW / kappa));
					cNE = exp(-(nablaNE / kappa) * (nablaNE / kappa));
					cSE = exp(-(nablaSE / kappa) * (nablaSE / kappa));
					cSW = exp(-(nablaSW / kappa) * (nablaSW / kappa));
					cNW = exp(-(nablaNW / kappa) * (nablaNW / kappa));
				}
				else if (2 == option)
				{
					cN = 1.0 / (1 + (nablaN / kappa) * (nablaN / kappa));
					cS = 1.0 / (1 + (nablaS / kappa) * (nablaS / kappa));
					cE = 1.0 / (1 + (nablaE / kappa) * (nablaE / kappa));
					cW = 1.0 / (1 + (nablaW / kappa) * (nablaW / kappa));
					cNE = 1.0 / (1 + (nablaNE / kappa) * (nablaNE / kappa));
					cSE = 1.0 / (1 + (nablaSE / kappa) * (nablaSE / kappa));
					cSW = 1.0 / (1 + (nablaSW / kappa) * (nablaSW / kappa));
					cNW = 1.0 / (1 + (nablaNW / kappa) * (nablaNW / kappa));
				}
 
				tmpbuffer2[i*height + j] += delta_t * ( 1.0 / (dy*dy)*cN * nablaN + 1.0 / (dy*dy)*cS * nablaS +1.0 / (dx*dx)*cE * nablaE + 1.0 / (dx*dx)*cW * nablaW +
					1.0 / (dd *dd)*cNE*nablaNE + 1.0 / (dd *dd)*cSE*nablaSE +1.0 / (dd *dd)*cSW*nablaSW + 1.0 / (dd*dd)*cNW*nablaNW);
			}
		}
		for (int kk = 0; kk < width*height; kk++)
			tmpbuffer1[kk] = tmpbuffer2[kk];
	} // 迭代结束  
 
	// 转换成图像数据 
	for (int kk = 0; kk < height*width; kk++)
		outbuffer[kk] = unsigned char(max(0, min(tmpbuffer2[kk], 255)));
 
	delete[] tmpbuffer2;
	tmpbuffer2 = NULL;
	delete[] tmpbuffer1;
	tmpbuffer1 = NULL;
}

python

import cv2
import numpy as np
import math
class anisodiff2D(object):
 
    def __init__(self, num_iter=5, delta_t=1/7, kappa=30, option=2):
 
        super(anisodiff2D, self).__init__()
 
        self.num_iter = num_iter
        self.delta_t = delta_t
        self.kappa = kappa
        self.option = option
 
        self.hN = np.array([[0, 1, 0], [0, -1, 0], [0, 0, 0]])
        self.hS = np.array([[0, 0, 0], [0, -1, 0], [0, 1, 0]])
        self.hE = np.array([[0, 0, 0], [0, -1, 1], [0, 0, 0]])
        self.hW = np.array([[0, 0, 0], [1, -1, 0], [0, 0, 0]])
        self.hNE = np.array([[0, 0, 1], [0, -1, 0], [0, 0, 0]])
        self.hSE = np.array([[0, 0, 0], [0, -1, 0], [0, 0, 1]])
        self.hSW = np.array([[0, 0, 0], [0, -1, 0], [1, 0, 0]])
        self.hNW = np.array([[1, 0, 0], [0, -1, 0], [0, 0, 0]])
 
    def fit(self, img):
 
        diff_im = img.copy()
 
        dx = 1; dy = 1; dd = math.sqrt(2)
 
        for i in range(self.num_iter):
            nablaN = cv2.filter2D(diff_im, -1, self.hN)
            nablaS = cv2.filter2D(diff_im, -1, self.hS)
            nablaW = cv2.filter2D(diff_im, -1, self.hW)
            nablaE = cv2.filter2D(diff_im, -1, self.hE)
            nablaNE = cv2.filter2D(diff_im, -1, self.hNE)
            nablaSE = cv2.filter2D(diff_im, -1, self.hSE)
            nablaSW = cv2.filter2D(diff_im, -1, self.hSW)
            nablaNW = cv2.filter2D(diff_im, -1, self.hNW)
 
            cN = 0; cS = 0; cW = 0; cE = 0; cNE = 0; cSE = 0; cSW = 0; cNW = 0
 
            if self.option == 1:
                cN = np.exp(-(nablaN/self.kappa)**2)
                cS = np.exp(-(nablaS/self.kappa)**2)
                cW = np.exp(-(nablaW/self.kappa)**2)
                cE = np.exp(-(nablaE/self.kappa)**2)
                cNE = np.exp(-(nablaNE/self.kappa)**2)
                cSE = np.exp(-(nablaSE/self.kappa)**2)
                cSW = np.exp(-(nablaSW/self.kappa)**2)
                cNW = np.exp(-(nablaNW/self.kappa)**2)
            elif self.option == 2:
                cN = 1/(1+(nablaN/self.kappa)**2)
                cS = 1/(1+(nablaS/self.kappa)**2)
                cW = 1/(1+(nablaW/self.kappa)**2)
                cE = 1/(1+(nablaE/self.kappa)**2)
                cNE = 1/(1+(nablaNE/self.kappa)**2)
                cSE = 1/(1+(nablaSE/self.kappa)**2)
                cSW = 1/(1+(nablaSW/self.kappa)**2)
                cNW = 1/(1+(nablaNW/self.kappa)**2)
 
            diff_im = diff_im + self.delta_t * (
 
                (1/dy**2)*cN*nablaN +
                (1/dy**2)*cS*nablaS +
                (1/dx**2)*cW*nablaW +
                (1/dx**2)*cE*nablaE +
 
                (1/dd**2)*cNE*nablaNE +
                (1/dd**2)*cSE*nablaSE +
                (1/dd**2)*cSW*nablaSW +
                (1/dd**2)*cNW*nablaNW
            )
 
        return diff_im