NoiseEstimate.rar NoiseEstimateDlg.cpp

www.pudn.com > NoiseEstimate.rar > NoiseEstimateDlg.cpp

// NoiseEstimateDlg.cpp : 实现文件
//

#include "stdafx.h"
#include "NoiseEstimate.h"
#include "NoiseEstimateDlg.h"
#include "HistShow.h"
#include "math.h"
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;

#ifdef _DEBUG
#define new DEBUG_NEW
#endif

#define blocksize 4 //局域标准差法估计噪声时的分块大小

// 用于应用程序“关于”菜单项的 CAboutDlg 对话框

class CAboutDlg : public CDialog
{
public:
CAboutDlg();

// 对话框数据
enum { IDD = IDD_ABOUTBOX };

protected:
virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV 支持

// 实现
protected:
DECLARE_MESSAGE_MAP()
};

CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)
{
}

void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
}

BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)
END_MESSAGE_MAP()

// CNoiseEstimateDlg 对话框

CNoiseEstimateDlg::CNoiseEstimateDlg(CWnd* pParent /*=NULL*/)
: CDialog(CNoiseEstimateDlg::IDD, pParent)
, m_delta(0)
, m_NoiseType(_T("高斯"))
, m_HistType(_T("原图"))
{
m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void CNoiseEstimateDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
DDX_Text(pDX, IDC_EDIT1, m_delta);
DDX_CBString(pDX, IDC_COMBO1, m_NoiseType);
DDX_CBString(pDX, IDC_COMBO2, m_HistType);

}

BEGIN_MESSAGE_MAP(CNoiseEstimateDlg, CDialog)
ON_WM_SYSCOMMAND()
ON_WM_PAINT()
ON_WM_QUERYDRAGICON()
//}}AFX_MSG_MAP
ON_BN_CLICKED(IDC_BUTTONOPEN, &amt;CNoiseEstimateDlg::OnBnClickedButtonopen)
ON_BN_CLICKED(IDC_BUTTONSAVE, &amt;CNoiseEstimateDlg::OnBnClickedButtonsave)
ON_BN_CLICKED(IDC_BUTTONNOISE, &amt;CNoiseEstimateDlg::OnBnClickedButtonnoise)
ON_BN_CLICKED(ID_NOISE_ESTIMATE, &amt;CNoiseEstimateDlg::OnBnClickedNoiseEstimate)
ON_BN_CLICKED(IDC_HISTSHOW, &amt;CNoiseEstimateDlg::OnBnClickedHistshow)
ON_BN_CLICKED(IDC_NOISE_ESTIMATE2, &amt;CNoiseEstimateDlg::OnBnClickedNoiseEstimate2)
END_MESSAGE_MAP()

// CNoiseEstimateDlg 消息处理程序

BOOL CNoiseEstimateDlg::OnInitDialog()
{
CDialog::OnInitDialog();

// 将“关于...”菜单项添加到系统菜单中。

// IDM_ABOUTBOX 必须在系统命令范围内。
ASSERT((IDM_ABOUTBOX &amt; 0xFFF0) == IDM_ABOUTBOX);
ASSERT(IDM_ABOUTBOX < 0xF000);

CMenu* pSysMenu = GetSystemMenu(FALSE);
if (pSysMenu != NULL)
{
CString strAboutMenu;
strAboutMenu.LoadString(IDS_ABOUTBOX);
if (!strAboutMenu.IsEmpty())
{
pSysMenu->AppendMenu(MF_SEPARATOR);
pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
}
}

// 设置此对话框的图标。当应用程序主窗口不是对话框时，框架将自动
// 执行此操作
SetIcon(m_hIcon, TRUE); // 设置大图标
SetIcon(m_hIcon, FALSE); // 设置小图标

// TODO: 在此添加额外的初始化代码

m_picR=m_picB=m_picG=NULL;
m_picO=new CFrameDisplay(this,IDC_STATIC);
m_picA=new CFrameDisplay(this,IDC_STATIC1);

m_sSaveFile="";

return TRUE; // 除非将焦点设置到控件，否则返回 TRUE
}

void CNoiseEstimateDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
if ((nID &amt; 0xFFF0) == IDM_ABOUTBOX)
{
CAboutDlg dlgAbout;
dlgAbout.DoModal();
}
else
{
CDialog::OnSysCommand(nID, lParam);
}
}

// 如果向对话框添加最小化按钮，则需要下面的代码
// 来绘制该图标。对于使用文档/视图模型的 MFC 应用程序，
// 这将由框架自动完成。

void CNoiseEstimateDlg::OnPaint()
{
CPaintDC dc(this);
m_picO->ShowImage();
m_picA->ShowImage();
}

//当用户拖动最小化窗口时系统调用此函数取得光标显示。
//
HCURSOR CNoiseEstimateDlg::OnQueryDragIcon()
{
return static_cast<HCURSOR>(m_hIcon);
}

void CNoiseEstimateDlg::OnBnClickedButtonopen()
{
// TODO: 在此添加控件通知处理程序代码
CString sFileName;
CFile fBMP;

if (m_picR) free(m_picR);
if (m_picG) free(m_picG);
if (m_picB) free(m_picB);
m_picR=m_picG=m_picB=NULL;

//选择并打开一个原始bmp文件
CFileDialog FileDlg(TRUE,NULL,NULL,OFN_HIDEREADONLY,"Picture File(*.bmp)|*.bmp");
if (FileDlg.DoModal()!=IDOK) return;
sFileName=FileDlg.GetPathName();
if (!fBMP.Open(sFileName,CFile::modeRead|CFile::typeBinary))
{
AfxMessageBox("文件无法打开");
return;
}
//读入文件，返回pos为图像文件头的长度
if (!Read_BMP(&amt;fBMP))
{
AfxMessageBox("文件格式不是无压缩24位色Windows 3格式");
fBMP.Close();
return;
}
fBMP.Close();

m_picO->Copy_RGB24(m_picB,m_picG,m_picR,m_dWidth,m_dHeight);
m_picO->ShowImage();
RedrawWindow();
}

void CNoiseEstimateDlg::OnBnClickedButtonsave()
{
// TODO: 在此添加控件通知处理程序代码
CFileDialog dlg(NULL,"","");
dlg.DoModal();//显示对话框
m_sSaveFile = dlg.GetFileName();
if (m_sSaveFile!="") m_picA->SaveImage24(m_sSaveFile);
m_sSaveFile="";
}

void CNoiseEstimateDlg::OnBnClickedButtonnoise()
{
// TODO: 在此添加控件通知处理程序代码
CWaitCursor waitcursor;

//UpdateData(TRUE);
//CString str;
//str.Format(">f",m_delta);
//MessageBox(str);

//添加噪声
int i,j;
//unsigned char *B,*G,*R;
UpdateData(TRUE);
//ofstream fout;
//fout.open( "C:\\Documents and Settings\\Administrator\\noise_50.txt" );
//if ( !fout.is_open() )
//{
// MessageBox( "Can't open file" );
// return;
//}

B = (unsigned char *)malloc(m_dWidth * m_dHeight);
G = (unsigned char *)malloc(m_dWidth * m_dHeight);
R = (unsigned char *)malloc(m_dWidth * m_dHeight);

if (m_NoiseType == "高斯")
{
//MessageBox( "程序正在运行中，请等待。。。" );
for(i = 0;i < m_dHeight;i++)
{
for(j = 0;j < m_dWidth;j++)
{
/**(B + i * m_dWidth + j) = *(m_picB + i*m_dWidth + j )
+ Discrete_Gauss_Noise(0, m_delta);*/
float result_G = (float)(*(m_picG + i * m_dWidth + j)) + Discrete_Gauss_Noise1(0, m_delta);

if (result_G > 255)
{
*(G + i * m_dWidth + j) = 255;
}
else if (result_G < 0)
{
*(G + i * m_dWidth + j) = 0;
}

else
{
*(G + i * m_dWidth + j) = (unsigned char) result_G;
}
/**(G + i * m_dWidth + j) = *(m_picG + i * m_dWidth + j)
+ Discrete_Gauss_Noise(0, m_delta);*/

/**(R + i * m_dWidth + j) = *(m_picR + i * m_dWidth + j)
+ Discrete_Gauss_Noise(0, m_delta);*/
//fout << setw(5);
//fout << Discrete_Gauss_Noise( 0, m_delta );

//下面三条语句用来生成一幅均匀灰度图像，灰度值由用户给出，大小和先打开的图像相同。
/**(B + i * m_dWidth + j) = (unsigned char)m_delta;
*(G + i * m_dWidth + j) = (unsigned char)m_delta;
*(R + i * m_dWidth + j) = (unsigned char)m_delta;*/
}//end j
//fout << setw(5);
//fout << endl;
}//end i
//fout << flush;
//fout.close();
//MessageBox( "添加噪声已经结束，可以继续其他操作！" );

////以下语句用来生成一幅灰度棋盘图像，棋盘格大小为8x8
//for (i = 0;i != m_dHeight;i += 8)
//{
// for (j = 0;j != m_dWidth;j += 8)
// {
// //每格内像素相同，格之间均匀（100-200）
// int k = 100 + rand() > 100;

// //黑白相间（0或者255）
// //int k = ((i + j) > 16 == 0) ? 0 : 255;

// for (int s = 0;s != 8;s++)
// {
// for (int t = 0;t != 8;t++)
// {
// *(B + (i + s) * m_dWidth + j + t) = k;
// *(R + (i + s) * m_dWidth + j + t) = k;
// *(G + (i + s) * m_dWidth + j + t) = k;
// }//end t
// }//end s
// }//end j
//}//end i
}
else if(m_NoiseType == "椒盐")
{
for(i = 0;i < m_dHeight;i++)
{
for(j = 0;j < m_dWidth;j++)
{
*(B + i * m_dWidth + j) = *(m_picB + i * m_dWidth + j);
*(G + i * m_dWidth + j) = *(m_picG + i * m_dWidth + j);
*(R + i * m_dWidth + j) = *(m_picR + i * m_dWidth + j);

int noise = rand() > 100;

if (noise < 10)
{
*(B + i * m_dWidth + j) = 0;
*(G + i * m_dWidth + j) = 0;
*(R + i * m_dWidth + j) = 0;
}
else if (noise > 40)
{
*(B + i * m_dWidth + j) = 255;
*(G + i * m_dWidth + j) = 255;
*(R + i * m_dWidth + j) = 255;
}
}
}
}
m_picA->Copy_RGB24(B,G,R,m_dWidth,m_dHeight);
m_picA->ShowImage();
RedrawWindow();
}

void CNoiseEstimateDlg::OnBnClickedNoiseEstimate()
{
//M2M4方法估计噪声

double kx,ky;
long double M2;
double N,sigma;

M2 = Two_Stage_Squares(G,m_dWidth * m_dHeight);

kx = Kurtosis(m_picG,m_dWidth * m_dHeight);
ky = Kurtosis(G,m_dWidth * m_dHeight);

if (fabs((kx - 3)) < 1)
{
MessageBox("原图为纯噪声图像，该算法不适用！");
return;
}
else if ((((ky - 3) / (kx - 3)) < 0))
{
MessageBox("((ky - 3) / (kx - 3)) < 0，该算法不适用！");
return;
}
else if((((ky - 3) / (kx - 3)) > 1))
{
MessageBox("((ky - 3) / (kx - 3)) > 1),该算法不适用！");
return;
}
else
{
N = M2 * (1 - sqrt((ky - 3) / (kx - 3)));
sigma = sqrt(N);

CString str;
str.Format(">f",sigma );
//str.Format(">f",ky);
MessageBox(str);
}
}

void CNoiseEstimateDlg::OnBnClickedNoiseEstimate2()
{
// 局域标准差法估计噪声
//首先计算出各图像块的局部标准差，然后找出频数最大的那个作为噪声估计

int lsdmax,lsdhistmax;
double lsdesitimate;
int M = int(m_dHeight * m_dWidth / blocksize / blocksize);
int N = blocksize * blocksize;
double* variance = new double [m_dHeight * m_dWidth / blocksize / blocksize];//存放计算出的各个局域方差
int* lsd = new int [m_dHeight * m_dWidth / blocksize / blocksize];//存放计算出的各个局域标准差的近似值（取整）
unsigned char* buf = new unsigned char [blocksize * blocksize];//临时数组，存放每一个图像块的数据

for (int i = 0; i <= m_dHeight;i += blocksize)
{
for (int j = 0;j <= m_dWidth;j += blocksize)
{
for (int s = 0; s != blocksize;s++)
{
for (int t = 0;t != blocksize;t++)
{
*(buf + s * blocksize + t) = *(G + (i + s) * blocksize + j + t);
}//end t
}//end s

*(variance + (i / blocksize) * (m_dWidth / blocksize) + (j / blocksize) ) = var(buf);//计算出了局域方差

}//end j
}//end i

for (int i = 0; i != M;i++)
{
*(lsd + i) = int(sqrt(*(variance + i)) * 100);//计算出了局域标准差
}

//求lsd的最大值
lsdmax = 0;lsdesitimate = 0;
for (int i = 0;i != M;i++)
{
if (lsdmax < *(lsd + i))
{
lsdmax = *(lsd + i);
}
}

//lsd的直方图数组lsdhist初始化
for (int i = 0;i != lsdmax;i++)
{
lsdhist[i] = 0;
}

//求直方图lsdhist(其实不是真正的直方图，因为不是概率，而是出现的次数)
for (int i = 0;i != M;i++)
{
lsdhist[*(lsd + i)]++;
}

//求直方图lsdhist的最大值和lsd的估计值
lsdhistmax = 0;
for (int i = 0;i != lsdmax;i++)
{
if (lsdhistmax < lsdhist[i])
{
lsdhistmax = lsdhist[i];
lsdesitimate = i;
}
}
CString str;
str.Format(">f",lsdesitimate / 100 * (double (sqrt(double((N - 1)/(N - 3))))));
AfxMessageBox(str);

}
double CNoiseEstimateDlg::var(unsigned char *buf)
{
//计算数组的方差值，数组大小为图像块的大小：blocksize * blocksize

double ave,var;
ave = 0;var = 0;

for (int i = 0; i != blocksize * blocksize;i++)
{
ave += *(buf + i);
}
ave /= (blocksize * blocksize);

for (int i = 0;i != blocksize * blocksize;i++)
{
var += pow((*(buf + i) - ave),2);
}
var /= (blocksize * blocksize - 1);

return var;
}

void CNoiseEstimateDlg::OnBnClickedHistshow()
{
//显示直方图
UpdateData();

DWORD height,width,max;
DWORD i;
unsigned char *p;

//求长度
height = m_dHeight;

//求宽度
width = m_dWidth;

//指向图像数据
if (m_HistType == "原图")
{
p = m_picG;
}
else if (m_HistType == "噪声图")
{
p = G;
}

if (p == NULL)
{
MessageBox("请先打开一幅图像。");
return;
}

//直方图数组初始化
for (i = 0;i < 256;i++)
{
hist[i] = 0;
}

//求直方图(其实不是真正的直方图，因为不是概率，而是出现的次数)
for (i = 0;i < height * width;i++)
{
hist[*(p + i)]++;
}

//求直方图最大值
max = 0;
for (i = 0;i < 256;i++)
{
if (max < hist[i])
{
max = hist[i];
}
}

//将直方图对最大值归一化
for (i = 0;i < 256;i++)
{
hist[i] = hist[i] * 280 / max;
}

CHistShow dlg(NULL,hist);
dlg.DoModal();

}

double CNoiseEstimateDlg::Two_Stage_Squares(unsigned char *Y,long n)
{
//求一组数的二阶矩（统计量）

double M2 = 0;
/*n = m_dWidth * m_dHeight;*/

//double EY = 0;
//for (int i = 0; i < n; i++)
//EY += *(Y + i);
//EY = EY / n;//均值

//M2 = 0;
for (int i = 0;i < n;i++)
//M2 += (*(Y + i) - EY) * (*(Y + i) - EY);
M2 += (*(Y + i)) * (*(Y + i));
//由于图像的像素值均值肯定大于0，不满足论文中均值为0的条件，所以这里计算的二阶矩是原点二阶矩，不是中心二阶矩。
M2 = M2 / n;

return (M2);
}

double CNoiseEstimateDlg::Kurtosis(unsigned char *Y,long n)
{
//求某一块(大小为n)的峰度系数

double k;
double M2,M4;
int i;

//double EY = 0;
//for (i = 0; i < n; i++)
//EY += *(Y + i);
//EY = EY / n;//均值

M2 = 0;M4 = 0;
for (i = 0;i < n;i++)
{
/*M2 += (*(Y + i) - EY) * (*(Y + i) - EY);
M4 += (*(Y + i) - EY) * (*(Y + i) - EY) * (*(Y + i) - EY) * (*(Y + i) - EY);*/
M2 += (*(Y + i)) * (*(Y + i));//同样，这里都是原点矩，不是中心矩。
M4 += (*(Y + i)) * (*(Y + i)) * (*(Y + i)) * (*(Y + i));
}
M2 = M2 / n; M4 = M4 / n;
k = M4 / M2; k = k / M2;

return (k);
}

int CNoiseEstimateDlg::Discrete_Gauss_Noise(double mu, double sigma)
{
//产生符合所给参数的高斯分布的离散值

int k;
double MAX=256;
double r,y,t;
t = 0;
unsigned int number;

for (int i = 0;i < MAX * 12;i++)
{
r = 1.0 * rand_s(&amt;number) / RAND_MAX;
t += r - 0.5;
}
//根据大数定理，上述操作产生了服从标准正态分布的值，范围是：[-MAX*12/2,MAX*12/2]

y = mu + sigma * t / sqrt(MAX);
if (y >= 0)
t = 1;
else
t = -1;

//取绝对值再加上0.5，使-0.5和0.5区间没有值存在（左右边分别移动0.5），避免了取整后0值加倍的情况
k = (int)((fabs(y) + 0.5) * t);
return(k);
}

int CNoiseEstimateDlg::Discrete_Gauss_Noise1(double mu, double sigma)
{
//产生符合所给参数的高斯分布的离散值
int noise;
//noise = floor(gen_normal(mu,sigma));
noise = Poisson(sigma);
//noise = floor(Normal(mu,sigma));
//noise = floor(gen_uniform(100,200));
return noise;
}

BOOL CNoiseEstimateDlg::Read_BMP(CFile *fBMP)
{
// 读入BMP图像，将各信息赋给响应的变量或指针

unsigned char temp[4];
int i,j,m;
long t;

fBMP->Read(temp,2);//首先读入两个字节
if ((temp[0] != 'B')&amt;&amt;(temp[1] != 'M')) return(FALSE);
//位图格式中文件头结构体的第一个变量是bfType,两个字节，代表文件类型，必须是0x424,
//即字符串"MB"

fBMP->Read(temp,4); fBMP->Read(temp,4); fBMP->Read(temp,4); fBMP->Read(temp,4);
if (Bytes_to_Long(temp,4) != 40) return(FALSE);
//读到了位图信息头结构体，第一个变量是biSize,WORD型，表示结构体的大小，必须是40个字节

fBMP->Read(temp,4); m_dWidth = Bytes_to_Long(temp,4);
//图像宽度

fBMP->Read(temp,4); m_dHeight = Bytes_to_Long(temp,4);
//图像高度

fBMP->Read(temp,2);
if (Bytes_to_Long(temp,2) != 1) return(FALSE);
//biPlanes必须是1

fBMP->Read(temp,2);
if (Bytes_to_Long(temp,2) != 24) return(FALSE);
//biBitCount是图像数据位数，这里要求是24，是真彩色图像

fBMP->Read(temp,4);
if (Bytes_to_Long(temp,4) != 0) return(FALSE);
//biCompression指示位图是否被压缩，这里要求非压缩

fBMP->Read(temp,4); fBMP->Read(temp,4); fBMP->Read(temp,4);
fBMP->Read(temp,4);
if (Bytes_to_Long(temp,4) != 0) return(FALSE);
//biClrUsed指示图像用到的颜色数，为0表示用到的颜色数是2的biBitCount次方

fBMP->Read(temp,4);
if (Bytes_to_Long(temp,4) != 0) return(FALSE);
//biClrImpression指示图像中重要的颜色数，为0则认为所有的颜色都重要

//分配内存
m_picB = (unsigned char *)malloc(m_dWidth * m_dHeight);
m_picG = (unsigned char *)malloc(m_dWidth * m_dHeight);
m_picR = (unsigned char *)malloc(m_dWidth * m_dHeight);

//读入图像数据

m = (4 - (m_dWidth * 3) > 4) > 4;
//避免图像宽度不是3的倍数，造成读入混乱，但是这样会忽略边缘的这些零头

for (i = m_dHeight;i > 0;i--)
{
//BMP文件从下到上，从左到右排列，最先读到的数据是最下面一行的左边的第一个像素
//这里i是递减的，则在数组中改变了这种排列方式，即变为从上到下，从左到右了
t=(i - 1) * m_dWidth;
for (j = 0;j < m_dWidth;j++)
{
fBMP->Read(temp,3);
//注意真彩色图像的3字节图像数据依次代表B,G,R分量
*(m_picB + t + j) = temp[0];
*(m_picG + t + j) = temp[1];
*(m_picR + t + j) = temp[2];
}
if (m != 0) fBMP->Read(temp,m);
}
return(TRUE);
}

long CNoiseEstimateDlg::Bytes_to_Long(unsigned char *temp, int n)
{
//把字节型数据-->长整型数据
long k=0;
for (;n>0;n--)
k=k*256+temp[n-1];
//一个字节一个字节地读，故乘于256
return(k);
}

void CNoiseEstimateDlg::Histogram_Equalizor(unsigned char *pointer)
{
//直方图均衡

int i,j,p[256];
long n,q[256];

n = m_dWidth * m_dHeight;

for (i = 0;i < 256;i++) p[i] = 0;

for (i = 0;i < n;i++)
p[*(pointer + i)]++;

for (i = 0;i < 256;i++)
{
q[i]=0;
for (j = 0;j <= i;j++)
q[i] += p[j];
}

for (i = 0;i < n;i++)
*(pointer + i) = (unsigned char)(q[*(pointer + i)] * 255 / n);
}