#include "stdafx.h" #include "NLmeans.h" #include //全局变量 int g_nSigma = 50; // 高斯曲线参数sigma int g_nSWBlocX = 9; // 搜索窗口，行方向 int g_nSWBlocY = 9; // 搜索窗口，列方向 int g_nPWBloc = 1; // 匹配窗口大小 // 高斯滤波不同权值拟合直线k、b，10组数据分别表示滤波参数sigma=10:10:100的拟合结果 int g_nGaussK[10] = {-159, -80, -53, -40, -32, -27, -23, -20, -18, -16}; int g_nGaussB[10] = {4100, 4114, 4118, 4120, 4121, 4122, 4123, 4123, 4123, 4124}; /*************************************************************************** * NlMeansDenoiseCenter() * 功能描述：NL-MEANS 图像去噪算法函数 一次只更新中心像素 * 输入参数： * unsigned char bitWitdth 输入数据位宽8位/16位 * unsigned short pSrc16 16位输入数据 * unsigned char pSrc8 8位输入数据 * int nHeight 图像高度 * int nWidth 图像宽度 * 输出参数：unsigned short pSrc16 16位输出数据 输入是16位数据，输出也为16位数据 * unsigned char pSrc8 8位输出数据 输入是8位数据，输出也为8位数据 * 返 回 值： * 调用关系 * 调用函数：图像边缘扩展ImageExpandINT() 图像边缘切割ImageExpandInverseINT() * 被调函数：ArithProcessFunc() ***************************************************************************/ void NlMeansDenoiseCenter(unsigned char bitWitdth, unsigned short *pSrc16, unsigned char* pSrc8, int nWidth, int nHeight) { //搜索窗口大小 int nSWBlocX = g_nSWBlocX; int nSWBlocY = g_nSWBlocY; //int nSWBlocX = 8; //int nSWBlocY = 8; //高斯滤波标准差 int nSigma = g_nSigma; //匹配窗口参数设置 int nPWBloc = g_nPWBloc; //nPWBloc = 1 //高斯线性拟合K、B int k = g_nGaussK[nSigma/10 - 1]; int b = g_nGaussB[nSigma/10 - 1]; //分配空间及初始化-->change into 16bit datablock int nLen = nWidth * nHeight; unsigned short *pImgTemp= new unsigned short[nLen]; memset(pImgTemp,0,sizeof(unsigned short) * nLen); if(bitWitdth == 8) { for(int i = 0;i < nLen;i++) { pImgTemp[i] = pSrc8[i]; } } else { memcpy(pImgTemp,pSrc16,sizeof(unsigned short) * nLen); } // 距离、权重相关的变量 int nWeight = 0; int nMaxWeight = 0; int nSumWeight = 0; int nDif = 0; int nSumDif = 0; double nDstValue = 0; // 图像扩展 将原始图像进行扩展 int nWidthNew = nWidth + 2 * (nSWBlocX + nPWBloc); int nHeightNew = nHeight + 2 * (nSWBlocY + nPWBloc); int nLenNew = nWidthNew * nHeightNew; //分配扩展后图像空间并初始化 unsigned short *pImgExpand = new unsigned short[nLenNew]; memset(pImgExpand,0,sizeof(unsigned short) * nLenNew); //存放降噪结果的数组 unsigned short *pImgOut = new unsigned short[nLenNew]; memset(pImgOut,0,sizeof(unsigned short) * nLenNew); //图像边缘扩展 ImageExpandINT(pImgTemp,nWidth,nHeight,pImgExpand,nWidthNew,nHeightNew); ////////////////////////////////////// PROCESS STARTS////////////////////////////////////////////////// // 针对每个像素进行去噪处理 int nSpos = 0; int nPpos = 0; for ( int y = nSWBlocY + nPWBloc; y < nHeightNew - nSWBlocY - nPWBloc; y++) //第一个去噪待处理点y坐标 原始图像的第一个点 pImgTemp(0,0) { for ( int x = nSWBlocX + nPWBloc ; x < nWidthNew - nSWBlocX - nPWBloc; x++)//第一个去噪待处理点x坐标 { nWeight = 0; nMaxWeight = 0; nSumWeight = 0; nDstValue = 0; //单元区域内处理，结果存放到WinsDenoised for (int j = y - nSWBlocY;j <= y + nSWBlocY;j++) { for ( int i = x - nSWBlocX ; i <= x + nSWBlocX; i++) { nDif = 0; nSumDif = 0; for (int s = -nPWBloc; s <= nPWBloc; s++) { nSpos = (j + s) * nWidthNew + i; //搜索窗口内的其他窗口 nPpos = (y + s ) * nWidthNew + x ;//待处理像素窗口 for (int r = -nPWBloc; r <= nPWBloc; r++) { nDif = pImgExpand[nSpos + r] - pImgExpand[nPpos + r]; nSumDif += abs(nDif); } } //通过拟合直线计算高斯权值 , k，b在拟合时已放大4096倍 if (nSumDif >= 3*nSigma) { nWeight = 0; } else { nWeight = (b - k * nSumDif); } if (nWeight < 0) { nWeight = 0; } else if (nWeight > 4096) { nWeight = 4096; } //记录最大权重fmaxWeight和总权重fTotalWeight nMaxWeight = __max(nMaxWeight , nWeight); nSumWeight += nWeight; nDstValue += nWeight*pImgExpand[j*nWidthNew+i]; } } //如果总权重为正数，则计算该单元区域中心点的估计值 if (nSumWeight > 0) { pImgOut[y * nWidthNew + x] = (unsigned short)(nDstValue / nSumWeight); } else { pImgOut[y * nWidthNew + x] = pImgExpand[y * nWidthNew + x]; } }//y }//x //截取去噪后图像 ImageExpandInverseINT(pImgOut,nWidth,nHeight,pImgOut,nWidthNew,nHeightNew); if(bitWitdth == 8) { for(int i = 0;i < nLen;i++) { pSrc8[i] =(unsigned char)pImgOut[i]; } } else { memcpy(pSrc16,pImgOut,sizeof(unsigned short) * nLen); } //释放内存 delete [] pImgTemp; delete [] pImgExpand; delete [] pImgOut; } /*************************************************************************** * ImageExpandINT() * 功能描述：图像边缘扩展函数，将原图像向周围扩展边缘 * 输入参数： * unsigned short* pSrc 输入数据 * int nHeight 输入图像高度 * int nWidth 输入图像宽度 * int NewWidth 边缘扩展后图像宽度 * int NewHeight 边缘扩展后图像高度 * * 输出参数：unsigned short *pDst 输出边缘扩展后数据 * 返 回 值：无 * 调用关系 * 调用函数：无 * 被调函数：NLM去噪算法NlMeansDenoiseCenter() ***************************************************************************/ void ImageExpandINT(unsigned short* pSrc,int nWidth,int nHeight,unsigned short *pDst,int NewWidth,int NewHeight) { int nstepW= (int)((NewWidth - nWidth) / 2); int nstepH = (int)((NewHeight - nHeight) / 2); int i = 0; int j = 0; for (i = nstepH;i < nstepH + nHeight;i++) { for(j = 0;j < nstepW;j++) { pDst[i * NewWidth+j] = pSrc[(i - nstepH) * nWidth + 0]; } for(j = nstepW;j < nstepW + nWidth;j++) { pDst[i * NewWidth + j] = pSrc[(i - nstepH) * nWidth + (j - nstepW)]; } for(j = nstepW + nWidth;j < NewWidth;j++) { pDst[i * NewWidth + j] = pSrc[(i - nstepH) * nWidth + nWidth - 1]; } } for(i = 0;i < nstepH;i++) { for( j = 0;j < NewWidth;j++) { pDst[i * NewWidth + j] = pDst[nstepH * NewWidth + j]; } } for(i = nstepH + nHeight;i < NewHeight;i++) { for(j = 0;j < NewWidth;j++) { pDst[i * NewWidth + j] = pDst[(nstepH + nHeight - 1) * NewWidth + j]; } } } /*************************************************************************** * ImageExpandInverseINT() * 功能描述：图像边缘切割, 将图像扩展的边缘去掉，获取与原图像相同大小和位置的数据 * 输入参数： * unsigned short* pSrc 输入数据 * int nHeight 输入图像高度 * int nWidth 输入图像宽度 * int NewWidth 边缘扩展后图像宽度 * int NewHeight 边缘扩展后图像高度 * * 输出参数：unsigned short *pDst 输出边缘扩展后数据 * 返 回 值：无 * 调用关系 * 调用函数：无 * 被调函数：NLM去噪算法NlMeansDenoiseCenter() ***************************************************************************/ void ImageExpandInverseINT(unsigned short* pSrc,int nWidth,int nHeight,unsigned short *pDst,int NewWidth,int NewHeight) { int nstepW = (int)((NewWidth - nWidth) / 2); int nstepH = (int)((NewHeight - nHeight) / 2); unsigned short * pos = NULL; pos = pDst + nstepH * NewWidth + nstepW; for (int i = 0;i < nHeight;i++) { memcpy(pSrc + i * nWidth, pos,sizeof(unsigned short) * nWidth); pos = pos + NewWidth; } }