PM模型用于图像降噪

#include<iostream>
#include<opencv2/opencv.hpp>
#include<math.h>
using namespace std;
using namespace cv;


Mat improved_PM(Mat srcImage);
Mat Demo_PM(Mat srcImage);
Mat Demo_PM (Mat srcImage,float threshold,int Cycle);
Mat Demo_Enhance(Mat srcImage);
int main()
{

//Mat srcImage=imread("F:/33.bmp",0);
Mat dst;
Mat original=imread("F:/高剂量CT图.bmp",0);
Mat srcImage=imread("F:/低剂量CT图.bmp",0);
// Mat srcImage=imread("F:\\病理图片\\示例图片\\分割\\308\\black0.jpg",0);
// Mat srcImage=imread("F:\\病理图片\\示例图片\\分割\\190\\black0.jpg",0);
// Mat srcImage=imread("F:\\病理图片\\示例图片\\分割\\278\\black278.jpg",0);
if(!srcImage.data)
return -1;
imshow("original",original);
imshow("srcImage",srcImage);
imwrite("srcImage.jpg",srcImage);
Mat srcImage1,PM_Image;
Mat gas;
Mat bilter;
/*
bilateralFilter(srcImage,bilter,5,60,60);


imshow("bilterFilter",bilter);
imwrite("bilter.jpg",bilter);
GaussianBlur(srcImage,gas,Size(3,3),1);
double gasps=PSNR(srcImage,gas);
cout<<"gas="<<gasps<<endl;
imshow("gas",gas);
imwrite("gas.jpg",gas);


*/
GaussianBlur(srcImage,gas,Size(3,3),1);
// double gasps=PSNR(srcImage,gas);
double gasps=PSNR(original,gas);
cout<<"gas="<<gasps<<endl;
imshow("gas",gas);
srcImage.convertTo(srcImage1,CV_32FC1);

//原始PM模型
PM_Image=Demo_PM(srcImage1,50,40);//10,13


//改进PM模型
//PM_Image=improved_PM(srcImage1);
PM_Image.convertTo(PM_Image,CV_8UC1);
imshow("dst",PM_Image);
imwrite("PM_image.jpg",PM_Image);

//double pm=PSNR(srcImage,PM_Image);
double pm=PSNR(original,PM_Image);
cout<<"pm="<<pm<<endl;




waitKey(0);
return 0;


}




/*
PM模型
*/
Mat Demo_PM (Mat srcImage,float threshold,int Cycle)
{
//float  threshold =50;
float lamda=1.0/25;
Mat tempImage;
srcImage.copyTo(tempImage);
for(int i=0;i<Cycle;i++)
{
Mat kernelN=(Mat_<float>(3,3)<<0,1,0,0,-1,0,0,0,0);
Mat kernelS=(Mat_<float>(3,3)<<0,0,0,0,-1,0,0,1,0);
Mat kernelE=(Mat_<float>(3,3)<<0,0,0,0,-1,1,0,0,0);
Mat kernelW=(Mat_<float>(3,3)<<0,0,0,1,-1,0,0,0,0);


Mat tempN=Mat::zeros(srcImage.size(),CV_32FC1);
Mat tempS=Mat::zeros(srcImage.size(),CV_32FC1);
Mat tempE=Mat::zeros(srcImage.size(),CV_32FC1);
Mat tempW=Mat::zeros(srcImage.size(),CV_32FC1);


/*
Mat tempN=srcImage*(-1);
Mat tempS=srcImage*(-1);
Mat tempE=srcImage*(-1);
Mat tempW=srcImage*(-1);
*/
filter2D(tempImage,tempN,CV_32FC1,kernelN);
filter2D(tempImage,tempS,CV_32FC1,kernelS);
filter2D(tempImage,tempE,CV_32FC1,kernelE);
filter2D(tempImage,tempW,CV_32FC1,kernelW);


Mat gN=Mat::zeros(srcImage.size(),CV_32FC1);
Mat gS=Mat::zeros(srcImage.size(),CV_32FC1);
Mat gE=Mat::zeros(srcImage.size(),CV_32FC1);
Mat gW=Mat::zeros(srcImage.size(),CV_32FC1);
/*
gN=1.0/(1.0+ ( (1.0/threshold) * tempN ) * ( ( 1.0/threshold ) * tempN ) );
gS=1.0/(1.0+((1.0/threshold)*tempS)*((1.0/threshold)*tempS));
gE=1.0/(1.0+((1.0/threshold)*tempE)*((1.0/threshold)*tempE));
gW=1.0/(1.0+((1.0/threshold)*tempW)*((1.0/threshold)*tempW));
*/
Mat north;
multiply(tempN,tempN,north);
gN=1.0/(1.0+( ( 1.0/(threshold*threshold) ) *north) );


Mat south;
multiply(tempS,tempS,south);
gS=1.0/(1.0+( ( 1.0/(threshold*threshold) ) *south) );


Mat east;
multiply(tempE,tempE,east);
gE=1.0/(1.0+( ( 1.0/(threshold*threshold) ) *east) );


Mat west;
multiply(tempW,tempW,west);
gW=1.0/(1.0+( ( 1.0/(threshold*threshold) ) *west) );


Mat N,S,E,W;
multiply(gN,tempN,N);
multiply(gS,tempS,S);
multiply(gE,tempE,E);
multiply(gW,tempW,W);
//tempImage=tempImage+(pow(lamda,i))*(N+S+W+E);

 

tempImage=tempImage+lamda*(N+S+W+E);

/*tempN.convertTo(tempN,CV_8UC1);
imshow("tempN",tempN);
gN.convertTo(gN,CV_8UC1);
imshow("gN",gN);  
waitKey(0);
*/
//tempImage=tempImage+lamda*(gN*tempN+gS*tempS+gE*tempE+gW*tempW);
}
return tempImage;


}