数字图像处理(DIP) C++

【数字图像处理】平滑与锐化smooth/sharpen

DIP平滑与锐化的C++实现

Posted by Raserting_L on July 20, 2019

开篇介绍

这篇文章我想主要说说空间滤波，滤波概念其实是频域概念，即对信号频率进行处理，高于或低于截止频率的将被干掉，或者带通带限，就有了高通滤波器，低通滤波器。频域的相乘对应于时域的卷积(我最早听这句话是在通信原理里，对于为什么我个学计算机的为啥要学通信原理。。。说起来，，，和学校有关，，，毕竟我们学校号称中国通信第一。。。所以，我们也得学，当时觉得我学计算机的学这玩意有啥用，还别说，这里真的用到了。。。)，于是，空域滤波器（空间滤波器也叫卷积核，空间掩膜，核，模板，窗口等）和图像的卷积能达到和频域相同或相近的效果，但值得注意的是空间滤波器只有线性滤波器和频域对应有关联，非线性滤波器在频域无法实现。

空间滤波器组成：

  1. 一个邻域（典型的是以某一点为中心的矩形）
  2. 对该邻域做规定的运算，得到结果赋值到邻域中心位置，特别的，赋值到邻域对应的中心并不是原图上的，而是结果上的，因为本步计算的结果将不会影响下一步计算，所以原图像不能被改变

如果执行线性操作，将可以在频域找到对应操作，如果采用非线性操作，频域则无对应。一般掩膜使用奇数x奇数大小，因为这样中心为整数坐标： 线性计算如下，w为模板，f为图像：

​

图像增强之 平滑

由于之后几乎每一个滤波器的使用都涉及卷积核在图像上移动的操作，这里先把移动操作部分函数贴出来

void RealCOV(double *src, double *dst, double *mask,  int width, int height, int m_width, int m_height){
    //double *Temp = new double[width * height] ;

    int mask_center_h = m_height / 2 ;
    int mask_center_w = m_width / 2 ;

    for(int i = 0 ; i < height ; i++ ){
        for(int j = 0 ; j < width ; j++ ){
            double value = 0.0 ;
            //确保不会越界
            if((i-mask_center_h)>=0 && (j-mask_center_w)>=0 &&(i+mask_center_h)<height && (j+mask_center_w)<width ){
            	for(int n = 0 ; n < m_height ; n++ ){
	                for(int m = 0 ; m < m_width ; m++ ){
						value += src[(i+n-mask_center_h)*width + (j+m-mask_center_w)] * mask[n*m_width+m] ;  
	                }
            	}
			}
            
            dst[i*width+j] = value ;
        }
    }
}

均值滤波

均值滤波模板算是比较简单，原理是，一个规定的邻域内，所有像素的平局值作为最终计算的结果，每个像素的权值相同，为总像素的倒数。

均值滤波的模板越大，图像越模糊。

代码实现如下：

//均值滤波器模板生成
void MeanMask(double *mask,int width,int height){
    double meanvalue = 1.0/(width * height) ;
    for(int i = 0 ; i < width * height ; i++ ){
        mask[i] = meanvalue ;
    }
}

void MeanFilter(double *src, double *dst, int width,int height,int m_width,int m_height){
    double *mask = new double[width * height] ;
    MeanMask(mask, m_width, m_height) ;
    RealCOV(src, dst, mask,  width, height, m_width, m_height) ; 
}

高斯滤波

高斯平滑是均值滤波的升级版本，邻域内每个像素的权值不同，而权值是由高斯函数确定的。 均值平滑和高斯平滑都是线性的，也就是，一旦参数给定，模板就确定下来，不会因为位置和像素分布不同而改变，而线性模板的基本运算是卷积。

一句话简介就是：高斯平滑就是一种加权的均值，权值由高斯函数确定。

高斯函数如下：

​

图像像下面这样：

  • 同等模板大小，标准差越大越模糊
  • 标准差相同，模板越大图像越模糊。

代码实现如下：

//高斯滤波器模板生成
void GaussianMask(double *mask,int width,int height,double deta){
    double deta_2 = deta * deta ;
    double mask_center_h = (double)height / 2 - 0.5 ; 
    double mask_center_w = (double)width / 2 - 0.5 ;
    double param = 1.0 / (2 * M_PI * deta_2) ;

    double sum = 0.0 ;

    for(int i = 0 ; i < height ; i++){
        for(int j = 0 ; j < width ; j++){
            double distance = sqrt((j-mask_center_w)*(j-mask_center_w)+(i-mask_center_h)*(i-mask_center_h)) ;
            mask[i * width + j] = param * exp(-(distance*distance)/(2*deta_2)) ;
            sum += mask[i * width + j] ;
        }
    }

    for(int i = 0 ; i < height * width ; i++){
        mask[i] /= sum ;
    }


}

void GaussianFilter(double *src, double *dst, int width, int height,int m_width,int m_height, double deta){
    double *mask = new double[width * height] ;
    GaussianMask(mask, m_width, m_height, deta) ;
    RealCOV(src, dst, mask,  width, height, m_width, m_height) ; 

}

中值滤波

中值滤波时典型的非线性方法，与前面介绍的方法不同，中值滤波更接近于灰度图像的腐蚀和膨胀，是在一定区域内比较大小，找出中值，也就是排序后中间那个数，也就是中学的中位数。

另外中值滤波对椒盐噪声和斑点噪声效果显著，而且中值滤波具有较好的边缘保持特性(反正书上给的使用中值滤波前后的效果图，很惊艳，像黑魔法)。

另外对于中值滤波算法的实现还有一个快速算法，这里先不做介绍，如果有兴趣可以参考博客

下面是代码实现(这里排序部分我简单使用了个冒泡，emmm，是Low了点，如果想改改的话也很方便，毕竟真正用起来一定是用其他很牛的快速的排序方法)。

//中值滤波器
void MedianFilter(double *src, double *dst, int width, int height, int m_width, int m_height){
    int mask_center_h = m_height / 2 ;
    int mask_center_w = m_width / 2 ;

    for(int i = 0 ; i < height ; i++ ){
        for(int j = 0 ; j < width ; j++ ){
            double value = 0.0 ;
            //确保不会越界
            if((i-mask_center_h)>=0 && (j-mask_center_w)>=0 &&(i+mask_center_h)<height && (j+mask_center_w)<width ){
                //对从src[(i-mask_center_h)*width+(j-mask_center_w)]开始的m_height*m_width个元素进行排序取中值
                int index = 0 ;
                double *sort = new double[m_width * m_height] ;
            	for(int n = 0 ; n < m_height ; n++ ){
	                for(int m = 0 ; m < m_width ; m++ ){
						sort[index] = src[(i+n-mask_center_h)*width + (j+m-mask_center_w)] ;
                        index ++ ;
	                }
            	}
                for(int n = 0 ; n < index-1 ; n++){
                    for(int m = 0 ; m < index-1-n ; m++){
                        if(sort[m] > sort[m+1]){
                            double temp_sort = sort[m] ;
                            sort[m] = sort[m+1] ;
                            sort[m+1] = temp_sort ;
                        }
                    }
                }
                value = sort[index/2] ;
			}
            
            dst[i*width+j] = value ;
        }
    }
}

图像锐化

简介

图像锐化是图像增强的一部分，增强的目的是使观察者看起来更容易识别某些模式，要观察的模式从频率域来分就有低频模式和高频模式，低频模式，也就是之前一直在讲的相对变化缓慢的部分，或者根本没有灰度变化的大片区域，高频部分，也就是灰度发生急剧变化的部分，其对应的就是图像边界，轮廓，还有就是。。。噪声。而锐化就是把这些边界，轮廓等灰度变化剧烈的地方强调出来。图像锐化步骤可以抽象为下面几步：

​ 原图——>提取图像细节——>突出原图细节——>锐化后的图像

第一步找细节，数字图像处理冈萨雷斯版第三版中给出了三种方法：一阶微分，二阶微分，非锐化掩蔽。

第二步突出灰度急剧变化的部分：将提取的部分乘预定的系数与原图像相加。

最后得到锐化后的图像。

数学基础

一阶微分

​

当然上面是对$x$的偏导，对于$y$的同理

二阶微分

​

同样，上面是对$x$的偏导，对于$y$同理

拉普拉斯算子

最后我们需要的就是将对$x$和$y$的二阶偏导都加起来得到拉普拉斯算子：

结合操作

最后的锐化公式是输出$g(x,y)$

其中的c为系数，使用的滤波器中心为负时，$c=-1$

使用的滤波器中心为正时$c=1$

代码

/*
 实数矩阵的180度旋转，对于Sobel算子卷积核分解为两个1*n矩阵时使用
 */
void RotateRealMatrix(double *matrix,double *dst,int width,int height){
    double *temp = new double[width*height] ;

    for(int i=0;i<width*height;i++){
        dst[width*height-1-i]=matrix[i];
    }
}

//卷积操作
void RealCOV(double *src, double *dst, double *mask,  int width, int height, int m_width, int m_height){
    //double *Temp = new double[width * height] ;

    int mask_center_h = m_height / 2 ;
    int mask_center_w = m_width / 2 ;

    for(int i = 0 ; i < height ; i++ ){
        for(int j = 0 ; j < width ; j++ ){
            double value = 0.0 ;
            //确保不会越界
            if((i-mask_center_h)>=0 && (j-mask_center_w)>=0 &&(i+mask_center_h)<height && (j+mask_center_w)<width ){
            	for(int n = 0 ; n < m_height ; n++ ){
	                for(int m = 0 ; m < m_width ; m++ ){
						value += src[(i+n-mask_center_h)*width + (j+m-mask_center_w)] * mask[n*m_width+m] ;  
	                }
            	}
			}
            
            dst[i*width+j] = value ;
        }
    }
}

void RealConvolution(double *src,double *dst,double *mask, int width,int height,int m_width,int m_height){
    double *temp = new double[width * height] ;
    RotateRealMatrix(mask,temp,m_width,m_height);
    RealCOV(src, dst, temp, width, height, m_width, m_height);
}

//最简单的Sobel算子操作
double Sobel(double *src, double *dst, int width, int height, int sobel_size){
    //最后一个参数可以指定卷积核大小，这里我们实现一个为3的并且为别的留下接口便于后面扩展
    //这里可以直接定义简单的3*3的卷积核，也可以将其分成两个1*3的。
    //另外注意对于不同的大小的卷积核，其分解得到的1*n矩阵也不同

    double *dst_x = new double[width * height] ;
    double *dst_y = new double[width * height] ;

    if(sobel_size == 3){
        //double SobelMask_x[3]={-1,-2,-1,0,0,0,1,2,1};
        double SobelMask1[3] = {0.25, 0.5, 0.25} ;
        double SobelMask2[3] = {1, 0, -1} ;
        RealConvolution(src, dst_y, SobelMask2, width, height, 1, 3) ;
        RealConvolution(dst_y, dst_y, SobelMask1, width, height, 3, 1) ;

        RealConvolution(src, dst_x, SobelMask1, width, height, 1, 3);
        RealConvolution(dst_x, dst_x, SobelMask2, width, height, 3, 1);
    }
    // else if(sobel_size==5){
    //     double SobelMask1[5]={1/16.,4/16.,6/16.,4/16.,1/16.};
    //     double SobelMask2[5]={1/3.,2/3.,0,-2/3.,-1/3.};
    //     RealConvolution(src, dst_x, SobelMask1, width, height, 1, 5);
    //     RealConvolution(dst_x, dst_x, SobelMask2, width, height, 5, 1);
        
    //     RealConvolution(src, dst_y, SobelMask2, width, height, 1, 5);
    //     RealConvolution(dst_y, dst_y, SobelMask1, width, height, 5, 1);
    
    // }else if(sobel_size==7){
    //     double SobelMask1[7]={1/64.,6/64.,15/64.,20/64.,15/64.,6/64.,1/64.};
    //     double SobelMask2[7]={0.1,0.4,0.5,0,-0.5,-0.4,-0.1};
    //     RealConvolution(src, dst_x, SobelMask1, width, height, 1, 7);
    //     RealConvolution(dst_x, dst_x, SobelMask2, width, height, 7, 1);
        
    //     RealConvolution(src, dst_y, SobelMask2, width, height, 1, 7);
    //     RealConvolution(dst_y, dst_y, SobelMask1, width, height, 7, 1);
        
    // }
    for(int i = 0 ; i < height ; i++ ){
        for(int j = 0 ; j < width ; j++ ){
            dst[i*width+j] = abs(dst_x[i*width+j]) + abs(-dst_y[i*width+j]) ; 
        }
    }
}

本博文参考：https://face2ai.com/dip-5-0-%E7%81%B0%E5%BA%A6%E5%9B%BE%E5%83%8F-%E7%A9%BA%E5%9F%9F%E6%BB%A4%E6%B3%A2%E5%9F%BA%E7%A1%80-%E5%8D%B7%E7%A7%AF%E5%92%8C%E7%9B%B8%E5%85%B3/

CATALOG
    FEATURED TAGS
    机器学习深度学习 前期博客迁移 推荐算法 统计学习方法 opencv+dlib 微信小程序 Knowledge_Graph OpenCV C++ SLAM 数字图像处理(DIP) C++ Spring Boot Java
    FRIENDS