Python语言不调用OpenCV函数实现Laplace算子锐化图像
 
编译环境
前言
离散Laplace滤波模板
不带对角项的滤波模板
带有对角项的扩展模板
Lpalace算子的数学定义
  
实例分析
代码展示
运行结果
结果对比
  
结论

 
编译环境
 
编程语言：Python
 IDE：PyCharm2017
 
前言
 
锐化处理的主要目的是突出灰度的过渡部分。图像锐化的用途多种多样，应用范围从电子印刷和医学成像到工业检测和军事系统的制导等。由之前的博客了解到，图像模糊可通过在空间域用像素领域平均法实现。因为均值处理和积分类似，在逻辑上，我们可以得出锐化处理可由空间微分来实现这一结论。基本上，微分算子的响应强度与图像在用算子操作的这一点的突变成都成正比，这样图像微风会增强边缘和其他突变（如噪声），削弱灰度变换缓慢的区域。
 
离散Laplace滤波模板
 
不带对角项的滤波模板
 
w=
 0 1 0
 1 -4 1
 0 1 0
 
带有对角项的扩展模板
 
w=
 1 1 1
 1 -8 1
 1 1 1
 
Lpalace算子的数学定义
 

 
 如果使用的模板中心系数为正时，c=1，如果中心系数为负时，c= -1
 
实例分析
 
下图是月球北极的一副略显模糊的图像，要对其进行空间锐化处理。
 
 
代码展示
 
import cv2
from pylab import *
import numpy as np

def Image_cal(Image,k):         # Image是输入图像，k是标定图像灰度最大值
    im=array(Image)
    img=[]
    fm=im-np.min(im)
    maxfm=np.max(fm)
    for i in range(im.shape[0]):
        for j in range(im.shape[1]):
            img[i,j]=k * (fm[i,j] / maxfm)
    return img
def Laplace(Image,a):                   #输入图像Image，模板a
    im=array(Image)
    img=array(Image)
    dim=math.sqrt(len(a))                   # 模板的维度dim
    w=im.shape[0]                           # 计算输入图像的宽高
    h=im.shape[1]
    b=[]                                    # 待处理的与模板等大小的图像块,分BGR通道
    g=[]
    r=[]
    if sum(a)==0:                           # 判断模板a的和是否为0
        A=1
    else:
        A=sum(a)
    if size(im.shape)==3:                   # 判断图像为灰度图像还是RGB图像
        for i in range(int(dim/2),w-int(dim/2)):
            for j in range(int(dim/2),h-int(dim/2)):
                for m in range(-int(dim/2),-int(dim/2)+int(dim)):
                    for n in range(-int(dim / 2), -int(dim / 2) + int(dim)):
                        b.append(im[i+m,j+n,0])
                        g.append(im[i+m,j+n,1])
                        r.append(im[i+m,j+n,2])
                img[i,j,0]=sum(np.multiply(np.array(a),np.array(b)))/A
                img[i, j, 1] =sum(np.multiply(np.array(a),np.array(g)))/A
                img[i, j, 2] =sum(np.multiply(np.array(a),np.array(r)))/A
                b=[];g=[];r=[]
    else:
        for i in range(int(dim/2),w-int(dim/2)):
            for j in range(int(dim/2),h-int(dim/2)):
                for m in range(-int(dim/2),-int(dim/2)+int(dim)):
                    for n in range(-int(dim / 2), -int(dim / 2) + int(dim)):
                        b.append(im[i+m,j+n])
                img[i,j]=sum(np.multiply(np.array(a),np.array(b)))/A
                b=[]
    return img
img=cv2.imread('book\\0338.tif')
a=[0,-1,0,-1,5,-1,0,-1,0]
im=Laplace(img,a)
cv2.imshow('Origin',img)
cv2.imshow('Laplace',im)
cv2.waitKey()
cv2.destroyAllWindows()
 
运行结果
 
 
结果对比
 

 可以看出处理之后的图像轮廓更清晰
 
结论
 
拉普拉斯算子是锐化数字图像的一种重要工具，能够有效的应用于各个领域。

       代码包括椒盐噪声，高斯噪声，均值滤波，中值滤波，高斯滤波，Sobel滤波，Laplace滤波和对应的系统函数以及三种常用的参数分析，MSE，PSNR，SSIM。
 
         太久没写博客了，懒得把知识点都写下来，有缘人看到需要报告的移步传送门去下载报告，报告上写的很详细。
 

 
图5 原图与取灰度值后的图片
 
 
 
 
图6 原图加信噪比0.2的椒盐噪声
 
 
图7原图加sigma=3的高斯噪声
 
 
图8 椒盐噪声，3*3均值滤波，系统3*3均值滤波
 
 
图9椒盐噪声，3*3均值滤波，7*7均值滤波
 
 
图10 高斯噪声，3*3均值滤波，5*5均值滤波
 

高斯滤波（加权均值滤波）：
 

   图11 高斯噪声，5*5高斯滤波，5*5系统高斯滤波
 
 
图12 高斯噪声，3*3高斯滤波，5*5高斯滤波
 
 
图13 椒盐噪声，3*3高斯滤波，5*5高斯滤波
 
 
 
中值滤波：
 

 图15 椒盐噪声，3*3中值滤波，3*3系统中值滤波
 
 
16 椒盐噪声，3*3中值滤波，5*5中值滤波 10*10中值滤波
 
 
 
Sobel算子滤波：
 
 
图17 原图，Sobel滤波，系统Sobel滤波
 
 
 
Laplace滤波：
 
 
图18 原图，Laplace滤波，系统Laplace滤波
 
#########所有函数###################
Mse 两个图片均方差
Psnr 两个图片峰值信噪比
Ssim 两个图片相似度
matrix_convolve 两矩阵卷积
add_salt_noise 加椒盐噪声	
add_gauss_noise	加高斯噪声
get_mid	取矩阵均值
get_ave 取矩阵中位数
mid_filter 中值滤波
mid_function 系统中值滤波函数
ave_filter 均值滤波
ave_function 系统均值滤波函数
gauss_filter_self 高斯滤波
gauss_function 系统高斯滤波函数
sobel_filter_self  Sobel算子锐化
sobel_function   系统Sobel算子锐化
laplacian_filter_self 拉普拉斯算子锐化
laplacian_function 系统拉普拉斯算子锐化
#######################################
import cv2 as cv
import numpy as np
import random
import math
import copy

###########均方误差(MSE)###########
def mse(pc1,pc2):
    mse=np.mean((pc1-pc2)**2)
    return float(mse)
###########峰值信噪比(PSNR)##########
# PSNR高于40dB说明图像质量极好（即非常接近原始图像）
# 在30—40dB通常表示图像质量是好的（即失真可以察觉但可以接受）
# 在20—30dB说明图像质量差
# 最后 PSNR低于20dB图像不可接受
def psnr(img1, img2):
    mse = np.mean( (img1 - img2) ** 2 )
    if mse == 0:
        return 100
    plxel_max = 255.0
    return 20 * math.log10(plxel_max / math.sqrt(mse))
###########结构相似度(SSIM):range~[-1:1]##########
#-1表示完全不相似,1表示完全相似
def ssim(y_true, y_pred):
    u_true = np.mean(y_true)
    u_pred = np.mean(y_pred)
    var_true = np.var(y_true)
    var_pred = np.var(y_pred)
    std_true = np.sqrt(var_true)
    std_pred = np.sqrt(var_pred)
    c1 = np.square(0.01 * 7)
    c2 = np.square(0.03 * 7)
    ssim = (2 * u_true * u_pred + c1) * (2 * std_pred * std_true + c2)
    denom = (u_true ** 2 + u_pred ** 2 + c1) * (var_pred + var_true + c2)
    return ssim / denom
###############矩阵卷积##############
def matrix_convolve(pc,mode):
    n,m=pc.shape
    c=np.zeros((n,m),dtype=np.float)
    step=mode.shape[0]
    mode=mode/mode.sum()  # 除于加权平均
    for i in range(n):
        for j in range(m):
            if i-int(step/2)<0 or i+int(step/2)>=n:
                c[i][j]=pc[i][j]
            elif j-int(step/2)<0 or j+int(step/2)>= m:
                c[i][j]=pc[i][j]
            else:
                x=int(step/2)
                c[i][j]=np.sum(pc[i-x:i-x+step,j-x:j-x+step]*mode)
    c=c.clip(0,255)
    c=cv.convertScaleAbs(c) #将结果转化为8位int
    return c
###############加椒盐噪声################
def add_salt_noise(pc,maybe):#图片,噪声比
    n,m=pc.shape
    for i in range(n):
        for j in range(m):
            if np.random.random(1)>maybe:
                continue
            else:
                pc[i,j]=0
###############加高斯噪声################
def add_gauss_noise(pc,mu,sigma,k,maybe=1):#means 均值 sigma 方差
    n,m=pc.shape
    for i in range(n):
        for j in range(m):
            if np.random.random(1)<=maybe:
                pc[i][j]+=k*random.gauss(mu,sigma)
                pc[i][j]=min(pc[i][j],255)
                pc[i][j]=max(pc[i][j],0)
##############求中值##############
def get_mid(pc,x,y,cnt):
    ans=[]
    for i in range(x-int(cnt/2),x+int(cnt/2)+1):
        for j in range(y-int(cnt/2),y+int(cnt/2)+1):
            ans.append(pc[i][j])
    ans.sort()
    return ans[int(len(ans)/2)+1]
##############求均值##############
def get_ave(pc,x,y,cnt):
    ans=0
    for i in range(x-int(cnt/2),x+int(cnt/2)+1):
        for j in range(y-int(cnt/2),y+int(cnt/2)+1):
            ans+=pc[i][j]
    return int(ans/cnt/cnt)
############中值滤波#############
def mid_filter(pc,step): #图片,几位滤波
    n,m=pc.shape
    c=np.zeros((n,m),dtype="uint8")
    for i in range(0,n):
        for j in range(0,m):
            if i-int(step/2)<0 or i+int(step/2)>=n:
                c[i][j]=pc[i][j]
            elif j-int(step/2)<0 or j+int(step/2)>=m:
                c[i][j]=pc[i][j]
            else:
                c[i][j]=get_mid(pc,i,j,step)
    return c
def mid_function(pc,step):
    c=cv.medianBlur(pc,step)
    return c
############均值滤波#############
def ave_filter(pc,step):
    n,m=pc.shape
    c = np.zeros((n, m), dtype="uint8")
    for i in range(n):
        for j in range(m):
            if i - int(step / 2) < 0 or i + int(step / 2) >= n:
                c[i][j] = pc[i][j]
            elif j - int(step / 2) < 0 or j + int(step / 2) >= m:
                c[i][j] = pc[i][j]
            else:
                c[i][j] = get_ave(pc,i,j,step)
    return c
def ave_function(pc,step):
    c=cv.blur(pc,(step,step))
    return c
###########高斯平滑#############
def gauss_filter_self(pc,step):
    if step==3:
        mode=np.array([[1,2,1],[2,4,2],[1,2,1]])
    if step==5:
        mode=np.array([[1,4,7,4,1],[4,16,26,16,4],[7,26,41,26,7],[4,16,26,16,4],[1,4,7,4,1]])
    return matrix_convolve(pc,mode)
def gauss_function(pc,step):
    c=cv.GaussianBlur(pc,(step,step),0)
    return c
###########Sobel###############
#sobel算子
# Gx =-1 0 1   Gy =1 2 1
#     -2 0 2       0 0 0
#     -1 0 1      -1-2-1
def sobel_filter_self(pc):
    c=copy.deepcopy(pc)
    n,m=pc.shape
    for i in range(1,n-1):
        for j in range(1,m-1):
            x=int(pc[i+1,j+1])-int(pc[i-1,j+1])+int(pc[i+1,j-1])-int(pc[i-1,j-1])+int(2*pc[i+1,j])-int(2*pc[i-1,j])
            y=int(pc[i+1,j+1])-int(pc[i+1,j-1])+int(pc[i-1,j+1])-int(pc[i-1,j-1])+int(2*pc[i,j+1])-int(2*pc[i,j-1])
            c[i,j]=min(255,int(math.sqrt(x*x+y*y)))
        c=cv.convertScaleAbs(c)
    return c
def sobel_function(pc):
    edges=cv.Sobel(pc,cv.CV_16S,1,1)
    edgesh=cv.convertScaleAbs(edges)
    return edgesh
###########Laplacian###############
#Laplacian算子
# 0  1  0
# 1 -4  1
# 0  1  0
def laplacian_filter_self(pc):
    c=copy.deepcopy(pc)
    n,m=pc.shape
    for i in range(1,n-1):
        for j in range(1,m-1):
            c[i,j]=abs(int(pc[i+1,j])+int(pc[i-1,j])+int(pc[i,j-1])+int(pc[i,j+1])-int(4*pc[i,j]))
            c[i,j]=min(255,c[i,j])
    c=cv.convertScaleAbs(c)
    return c
def laplacian_function(pc):
    edges=cv.Laplacian(pc,-1)
    return edges
################################
begin=cv.imread('E:/PC/4.jpg')
im=cv.imread('E:/PC/4.jpg',0) #读取图片并取灰度值
cv.imshow("begin",begin)
a=copy.deepcopy(im)
b=copy.deepcopy(im)
cv.imshow("initial",im) #输入图片
add_salt_noise(im,0.2) #椒盐噪声
add_gauss_noise(im,3,10,3) #高斯噪声
cv.imshow("add_noise",im)
add_gauss_noise(b,3,5,3)
cv.imshow("add_gauss_noise",b)

ima=ave_filter(im,3) #3位均值滤波
imb=ave_function(im,3) #3位系统均值滤波
imc=ave_filter(im,5) #5位均值滤波
imd=gauss_filter_self(im,5) #5位高斯滤波
ime=gauss_function(im,3) #3位系统高斯滤波
imf=mid_filter(im,3)  #3位中值滤波
img=mid_function(im,3) #3位系统中值滤波
imh=mid_filter(im,5) #5位中值滤波
cv.imshow("ave_self_3*3",ima)
cv.imshow("ave_function_3*3",imb)
cv.imshow("ave_self_5*5",imc)
cv.imshow("gauss_self_3*3",imd)
cv.imshow("gauss_function_3*3",ime)
cv.imshow("mid_self_3*3",imf)
cv.imshow("mid_function_3*3",img)
cv.imshow("mid_self_5*5",imh)
im1=sobel_filter_self(im)#Sobel算子
im2=sobel_function(im)#系统Sobel算子
im3=laplacian_filter_self(im)#laplace算子
im4=laplacian_function(im)#系统laplace算子
cv.imshow("sobel_filter_self",im1)
cv.imshow("sobel_function",im2)
cv.imshow("laplacian_filter_self",im3)
cv.imshow("laplacian_function",im4)
print("init:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(a,a),psnr(a,a),ssim(a,a)))
print("add_salt:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(im,a),psnr(im,a),ssim(im,a)))
print("a:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(ima,a),psnr(ima,a),ssim(ima,a)))
print("b:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(imb,a),psnr(imb,a),ssim(imb,a)))
print("c:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(imc,a),psnr(imc,a),ssim(imc,a)))
print("d:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(imd,a),psnr(imd,a),ssim(imd,a)))
print("e:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(ime,a),psnr(ime,a),ssim(ime,a)))
print("f:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(imf,a),psnr(imf,a),ssim(imf,a)))
print("g:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(img,a),psnr(img,a),ssim(img,a)))
print("h:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(imh,a),psnr(imh,a),ssim(imh,a)))
print("im1:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(im1,a),psnr(im1,a),ssim(im1,a)))
print("im2:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(im2,a),psnr(im2,a),ssim(im2,a)))
print("im3:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(im3,a),psnr(im3,a),ssim(im3,a)))
print("im4:MSE=%.2f PSNR=%.2f SSIM=%.2f"%(mse(im4,a),psnr(im4,a),ssim(im4,a)))
cv.waitKey(0)

 
 

 
目录
 
1 原理
1.1 中值滤波
1.2 均值滤波
1.3 图像锐化
1.4 边缘检测
  
2 实现源代码（MATLAB）
2.1 中值滤波
2.2 均值滤波
2.3 锐化处理
2.3.0 说明
2.3.1 Laplace算子
2.3.2 Sobel算子
2.3.3 Prewitt算子
2.3.4 Roberts算子
   
2.4 边缘检测
  
3 实验结果
3.1 中值滤波（使用3*3大小模板）
3.2 均值滤波（使用3*3大小模板）
3.3 锐化处理（注：锐化效果未叠加到原图，展示图为轮廓图）
3.3.0原图
3.3.1 Laplace算子
3.3.2 Sobel算子
3.3.3 Prewitt算子
3.3.4 Roberts算子
   
3.4 边缘检测
 

 
1 原理
 
1.1 中值滤波
 
主要用来处理椒盐噪声。椒盐噪声的出现使得该点像素比周围的像素亮（暗）很多，而如果在某个模板中对像素由小到大进行重新排列，那么最亮的或者最暗的点一定被排在两侧。这时取模板中排在中间位置上的像素的灰度值替代待处理像素的值，就可以达到滤除噪声的目的。
 
将模板中心与像素位置重合
读取模板下各对应像素的灰度值
将这些像素从小到大排成1列
找出这些值里排在中间的1个
将这个中间值赋给模板中心位置像素。
 
1.2 均值滤波
 
采用邻域平均法。基本思想是用几个像素灰度的平均值来代替每个像素的灰度。如使用一个3*3大小的模板对图像进行滤波，每次移动模板后，对模板范围内的所有像素灰度值相加再除以系数9，四舍五入后作为模板中心像素的灰度值。
 
1.3 图像锐化
 
图像锐化的目的是为了突出图像的边缘信息，加强图像的轮廓特征，便于人眼或者机器的识别。图像中边缘和轮廓往往出现于图像中灰度图片的地方，而检查这些图片可以用微分实现，因而图像锐化主要是通过微分方法进行的。
 由于处理的对象是数字图像，其像素是离散的，对微分形式进行一些近似处理，得到离散的模板矩阵，称为算子。常见的算子诸如Laplacians、Roberts、Sobel等。用这些算子作为模板对图像进行处理，就能得到图像的轮廓图。将轮廓图与原图像相加即可得到原图像的锐化结果。
 
1.4 边缘检测
 
在图像锐化得到的轮廓图的基础上，根据轮廓图的灰度直方图设立阈值（一般在两峰一谷之间），对其进行二值化处理即可。
 
2 实现源代码（MATLAB）
 
2.1 中值滤波
 
function res = medianfilter(mat, width, height)
% res = medianfilter(img, width, height)
% -- res        the result of the median filter process
% -- img        the input image matrix(should be 2 dimensions)
% -- width      the width of template
% -- height     the height of template
% e.g.:
%   filename='noise.jpg';
%   img=imread(filename);
%   width=3;height=3;
%   res = medianfilter(img, width, height);
%   imshow(res);
%   mkdir('results/medianfilter');
%   imwrite(res,['./results/medianfilter/',filename]);
 
res=im2uint8(zeros(size(mat)));
[rows,cols]=size(mat);
 
for i = 1:(rows-height + 1)
    for j = 1:(cols-width + 1)
        temp = mat(i:i+height-1, j:j+width-1);
        temp = sort(temp(:));
        res((i+i+height-1)/2, (j+j+width-1)/2) = temp((length(temp)+1)/2);
    end
end
 
 
2.2 均值滤波
 
function res = meanfilter(mat, width, height)
% res = medianfilter(img, width, height)
% -- res        the result of the median filter process
% -- img        the input image matrix(should be 2 dimensions)
% -- width      the width of template
% -- height     the height of template
% e.g.:
%   filename='noise.jpg';
%   img=imread(filename);
%   width=3;height=3;
%   res = meanfilter(img, width, height);
%   imshow(res);
%   mkdir('results/meanfilter');
%   imwrite(res,['./results/meanfilter/',filename]);
 
res=im2uint8(zeros(size(mat)));
[rows,cols]=size(mat);
 
for i = 1:(rows-height + 1)
    for j = 1:(cols-width + 1)
        temp = mat(i:i+height-1, j:j+width-1);
        res((i+i+height-1)/2, (j+j+width-1)/2) = round(mean(temp(:)));
    end
end
 
 
2.3 锐化处理
 
2.3.0 说明
 
在分别实现使用Sobel、Prewitt、Roberts、Laplace算子的滤波函数后，希望实现一个“万能滤波函数”。
 搜索得知MATLAB内置的滤波函数为imfilter，支持输入图像、模板等参数，输出滤波后结果。其对“最外面一圈像素”的处理方式是在最外层“加一圈0像素”，并对全图像进行处理。
 仿照其功能实现myfilter函数，代码如下：
 
function res=myfilter(I,H)
% function res=myfilter(I,H)
% -- res    the result of the filter process
% -- I      the input image
% -- H      the input template matrix(should be 2 dimensions)
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   I=imread(filename);
%   H=[1,1,1;1,-8,1;1,1,1];
%   res=myfilter(I,H);
%   imshow(res);
%   mkdir('results/myfilter');
%   imwrite(res,['./results/myfilter/',filename]);
I=im2uint8(I);
H=int16(H);
[rows,cols,channels]=size(I);
res=uint8(zeros(rows,cols,channels));
[height,width]=size(H);
hh=(height+1)/2;
hw=(width+1)/2;
for k=1:channels
    chan=int16(zeros(rows+hh,cols+hw));
    chan(hh:rows+hh-1,hw:cols+hw-1)=I(:,:,k);
    for i=1:rows
        for j=1:cols
            temp=chan(i:i+hh,j:j+hw);
            res(i,j,k)=sum(sum(H.*temp));
        end
    end
end
 
通过如下代码实现两者效果对比：
 
I=imread('锐化及边缘检测用途.jpg');
H=[1,1,1;1,-8,1;1,1,1];
imres=imfilter(I,H);
myres=myfilter(I,H);
cmp=(imres==myres);
length(cmp(cmp==0))
 
得到结果为：
 
ans =

     0
 
于是复用myfilter函数，以另一种方式实现了Sobel、Prewitt、Roberts、Laplace算子滤波函数，在原命名后添加了“2”以区分。
 需要注意的是，myfilter最多算是imfilter的一个弱化版本，运行速度远不及imfilter。另外，经查阅，imfilter默认使用的是’corr’参数，而不是’conv’，即默认使用的不是卷积运算。但本实验遇到的算子（Sobel、Prewitt、Laplace）模板都可以认为是对称的，实验得知本实验中对这三类算子在’conv’和’corr’两种参数下结果是完全一致的，故myfilter也没有提供诸如’corr’和’conv’之类参数的支持（即不支持模板矩阵不对称的卷积运算，仅支持相关运算）。
 
2.3.1 Laplace算子
 
采用的laplacian operator为[1,1,1;1,-8,1;1,1,1]，故若要得到锐化图像，应该使用原图像减去该函数计算结果再输出。
 
function laplacian=laplacianfilter(image)
% function laplacian=laplacianfilter(image)
% -- laplacian      the result of the laplacian operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=laplacianfilter(img);
%   imshow(res);
%   mkdir('results/laplacianfilter');
%   imwrite(res,['./results/laplacianfilter/',filename]);
 
image=im2uint8(image);
[rows,cols,channels]=size(image);
laplacian=uint8(zeros(rows,cols,channels));
 
for k=1:channels
    chan=int16(zeros(rows+2,cols+2));
    chan(2:rows+1,2:cols+1)=image(:,:,k);
    for i=1:rows
        for j=1:cols
            temp=chan(i:i+2, j:j+2);
            laplacian(i,j,k)=sum(temp(:))-9*chan(i+1,j+1);
        end
    end
end
 
复用myfilter函数实现：
 
function laplacian=laplacianfilter2(image)
% function laplacian=laplacianfilter2(image)
% -- laplacian      the result of the laplacian operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=laplacianfilter2(img);
%   imshow(res);
%   mkdir('results/laplacianfilter');
%   imwrite(res,['./results/laplacianfilter/',filename]);
laplacianoperator=[1,1,1;1,-8,1;1,1,1];
laplacian=myfilter(image,laplacianoperator);
 
 
2.3.2 Sobel算子
 
function sobel=sobelfilter(image)
% function sobel=sobelfilter(image)
% -- sobel          the result of the sobel operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=sobelfilter(img);
%   imshow(res);
%   mkdir('results/sobelfilter');
%   imwrite(res,['./results/sobelfilter/',filename]);
sobeloperator1=int16([-1,0,1;-2,0,2;-1,0,1]);
sobeloperator2=int16([-1,-2,-1;0,0,0;1,2,1]);
image=im2uint8(image);
[rows,cols,channels]=size(image);
sobel1=uint8(zeros(rows,cols,channels));
sobel2=uint8(zeros(rows,cols,channels));
 
for k=1:channels
    chan=int16(zeros(rows+2,cols+2));
    chan(2:rows+1,2:cols+1)=image(:,:,k);
    for i=1:rows
        for j=1:cols
            temp=chan(i:i+2,j:j+2);
            sobel1(i,j,k)=sum(sum(sobeloperator1.*temp));
            sobel2(i,j,k)=sum(sum(sobeloperator2.*temp));
        end
    end
end
sobel=abs(sobel1)+abs(sobel2);
 
复用myfilter函数实现：
 
function sobel=sobelfilter2(image)
% function sobel=sobelfilter2(image)
% -- sobel          the result of the sobel operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=sobelfilter2(img);
%   imshow(res);
%   mkdir('results/sobelfilter');
%   imwrite(res,['./results/sobelfilter/',filename]);
sobeloperator1=[-1,0,1;-2,0,2;-1,0,1];
sobeloperator2=[-1,-2,-1;0,0,0;1,2,1];
sobel=myfilter(image,sobeloperator1)+myfilter(image,sobeloperator2);
 
 
2.3.3 Prewitt算子
 
function prewitt=prewittfilter(image)
% function prewitt=prewittfilter(image)
% -- prewitt        the result of the prewitt operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=prewittfilter(img);
%   imshow(res);
%   mkdir('results/prewittfilter');
%   imwrite(res,['./results/prewittfilter/',filename]);
image=im2uint8(image);
[rows,cols,channels]=size(image);
prewitt1=uint8(zeros(rows,cols,channels));
prewitt2=uint8(zeros(rows,cols,channels));
for k=1:channels
    chan=int16(zeros(rows+2,cols+2));
    chan(2:rows+1,2:cols+1)=image(:,:,k);
    for i=1:rows
        for j=1:cols
            temp=chan(i+2,j:j+2)-chan(i,j:j+2);
            prewitt1(i,j,k)=sum(temp(:));
            temp=chan(i:i+2,j+2)-chan(i:i+2,j);
            prewitt2(i,j,k)=sum(temp(:));
        end
    end
end
 
prewitt=prewitt1+prewitt2;
 
复用myfilter函数实现：
 
function prewitt=prewittfilter2(image)
% function prewitt=prewittfilter2(image)
% -- prewitt        the result of the prewitt operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=prewittfilter2(img);
%   imshow(res);
%   mkdir('results/prewittfilter');
%   imwrite(res,['./results/prewittfilter/',filename]);
prewittoperator1=[-1,-1,-1;0,0,0;1,1,1];
prewittoperator2=[-1,0,1;-1,0,1;-1,0,1];
prewitt=myfilter(image,prewittoperator1)+myfilter(image,prewittoperator2);
 
2.3.4 Roberts算子
 
为方便计算，将Roberts算子修改为3*3矩阵，复用myfilter函数实现：
 
function roberts=robertsfilter2(image)
% function roberts=robertsfilter2(image)
% -- roberts        the result of the laplacian operator filter process
% -- image          the input image
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   res=robertsfilter2(img);
%   imshow(res);
%   mkdir('results/robertsfilter2');
%   imwrite(res,['./results/robertsfilter2/',filename]);
robertsoperator1=[0,0,0;0,-1,0;0,0,1];
robertsoperator2=[0,0,0;0,0,-1;0,1,0];
roberts=myfilter(image,robertsoperator1)+myfilter(image,robertsoperator2);
 
 
2.4 边缘检测
 
同样地，参考MATLAB内置边缘检测函数edge实现了其一个简化版本myedge：支持直接输入模板矩阵，或者输入字符串以调用对应模板矩阵；没有自动阈值选择，需要手动输入）。代码如下：
 
function res=myedge(img,operator,threshold)
% function res=myedge(img,operator,threshold)
% -- res        the edge result of the input image
% -- image      the input image
% -- operator   the input operator(should be 2 dimensional odd-order matrix) or string like 'laplacian'/'sobel'/'prewitt'/'roberts'
% -- threshold  threshold to decide whther a pixel on the edge should be black or white 
% e.g.:
%   filename='锐化及边缘检测用途.jpg';
%   img=imread(filename);
%   operator='laplacian';
%   threshold=233;
%   res=myedge(img,operator,threshold);
%   imshow(res);
%   mkdir(['results/myedge/',operator]);
%   imwrite(res,['./results/myedge/',operator,'/',filename]);
if ischar(operator)
    switch operator
        case 'laplacian'
            res=laplacianfilter2(img);
        case 'sobel'
            res=sobelfilter2(img);
        case 'prewitt'
            res=prewittfilter2(img);
        case 'roberts'
            res=robertsfilter2(img);
    end
else
    res=myfilter(img,operator);
end
res(res<threshold)=0;res(res>=threshold)=255;
 
3 实验结果
 
3.1 中值滤波（使用3*3大小模板）
 
 
编写了通用函数，可以在输入参数中指定模板大小。
 
3.2 均值滤波（使用3*3大小模板）
 
 
编写了通用函数，可以在输入参数中指定模板大小。
 
3.3 锐化处理（注：锐化效果未叠加到原图，展示图为轮廓图）
 
3.3.0原图
 
 
3.3.1 Laplace算子
 
 
3.3.2 Sobel算子
 
 
3.3.3 Prewitt算子
 
 
3.3.4 Roberts算子
 
 
 
3.4 边缘检测
 
以laplacian operator为模板矩阵，阈值为233。其余结果类似，可参照help myedge中的例子，调用myedge函数生成。
 
 
注：遇到的一个问题是，边缘检测应当是只有0和255两种灰度的图像，在MATLAB中生成以后使用imshow查看没有问题，但保存为.jpg格式后出现了其他灰色。再次读取该.jpg图片，数值上与原来MATLAB生成图像并不一致，猜测可能是压缩编码为.jpg格式时出现了损失。

实验题目:
 
编程实现灰度图像的中值滤波平滑处理。滤波模板的大小自定（可为3×3、5×5、7×7、15×15等）。实验图像可从提供的实验图像集中的噪声图像中选取。
 思考题：（选做）
 编程实现灰度图像的均值滤波平滑处理；也可尝试实现灰度图像的锐化处理，包括Sobel、Prewitt、Roberts、Laplace、Canny边缘检测等。
 
源码:
 
均值滤波
 
function ava(n)
imNew=im2double(imread('noise.jpg'));
len=floor(n/2);
a(1:n,1:n)=1;   %a即n×n模板,元素全是1
%对原始图像进行扩展，复制边界的值
imNew_pad=padarray(imNew,[len,len],'symmetric');
[M,N]=size(imNew_pad);
New=zeros(size(imNew));
for i=1+len:M-len
    for j=1+len:N-len
        c=imNew_pad(i-(n-1)/2:i+(n-1)/2,j-(n-1)/2:j+(n-1)/2).*a; %取出x1中从(i,j)开始的n行n列元素与模板相乘  
        s=sum(sum(c));                 %求c矩阵中各元素之和  
        New(i-(n-1)/2,j-(n-1)/2)=s/(n*n); %将与模板运算后的各元素的均值赋给模板中心位
    end
end
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(New);
title(['灰度图像的均值滤波平滑处理',num2str(n),'X',num2str(n),'图像']);
end
 
中值滤波
 
// 中值函数
function mid=mid(A)
len=length(A);
for i = 1:len
    for j = 1:len-i
        if A(j) > A(j+1)
            temp = A(j);       %核心代码
            A(j) = A(j+1);
            A(j+1) = temp;
        end
    end
end
mid=A(floor(len/2)+1);
end

//中值滤波 
function mido(bianSize)
imNew=im2double(imread('noise.jpg'));
%扩展区域的行列数，floor取整
len=floor(bianSize/2);
%对原始图像进行扩展，没有参数默认0填充
imNew_pad=padarray(imNew,[len,len]);
[M,N]=size(imNew_pad);
New=zeros(size(imNew));
for i=1+len:M-len
    for j=1+len:N-len
        %从扩展图像中，取出局部图像
        Block=imNew_pad(i-len:i+len,j-len:j+len);
        %将多维矩阵转换为一维数组
        Block=Block(:);
        %取这组数的中值,赋值给输出图像        
        New(i-len,j-len)=mid(Block);
    end
end
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(New);
title(['灰度图像的中值滤波平滑处理',num2str(bianSize),'X',num2str(bianSize),'图像']);
end
 
Prewitt算子
 
function Prewitt
imNew=im2double(imread('锐化及边缘检测用途.jpg'));
[M , N, R]=size(imNew);
for i=2:M-1
    for j=2:N-1
        Dx=[imNew(i+1,j-1)-imNew(i-1,j-1)]+[imNew(i+1,j)-imNew(i-1,j)]+[imNew(i+1,j+1)-imNew(i-1,j+1)];
        Dy=[imNew(i-1,j+1)-imNew(i-1,j-1)]+[imNew(i,j+1)-imNew(i,j-1)]+[imNew(i+1,j+1)-imNew(i+1,j-1)];
       P(i,j)=sqrt(Dx^2+Dy^2);
     
    end
end
for i=1:M-1
    for j=1:N-1
        if (P(i,j)<0.5)
            P(i,j)=1;
        else P(i,j)=0;
        end
    end
end
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(P,[]);
title('Prewitt边缘检测');  %画出边缘检测后的图像
end

 
Roberts算子
 
function Roberts
imNew=im2double(imread('锐化及边缘检测用途.jpg'));
[M,N]=size(imNew);
newimNew=imNew;%为保留图像的边缘一个像素
for j=1:M-1 %进行边界提取
    for k=1:N-1
        robertsNum = abs(imNew(j,k)-imNew(j+1,k+1)) + abs(imNew(j+1,k)-imNew(j,k+1));
        newimNew(j,k)=robertsNum;
    end
end
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(newimNew);
title('Roberts边缘检测');  %画出边缘检测后的图像
end


 
Sobel算子
 
function Sobel
imNew=im2double(imread('锐化及边缘检测用途.jpg'));
[M,N] = size(imNew);   % 获得图像的高度和宽度
F2 = double(imNew);        
U = double(imNew);       
uSobel = imNew;
for i = 2:M - 1   %sobel边缘检测
    for j = 2:N - 1
        Gx = (U(i+1,j-1) + 2*U(i+1,j) + F2(i+1,j+1)) - (U(i-1,j-1) + 2*U(i-1,j) + F2(i-1,j+1));
        Gy = (U(i-1,j+1) + 2*U(i,j+1) + F2(i+1,j+1)) - (U(i-1,j-1) + 2*U(i,j-1) + F2(i+1,j-1));
        uSobel(i,j) = sqrt(Gx^2 + Gy^2); 
    end
end 
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(im2uint8(uSobel));
title('Sobel边缘检测');  %画出边缘检测后的图像
end
 
交互面板
 
function varargout = mmm2(varargin)
% MMM2 MATLAB code for mmm2.fig
%      MMM2, by itself, creates a new MMM2 or raises the existing
%      singleton*.
%
%      H = MMM2 returns the handle to a new MMM2 or the handle to
%      the existing singleton*.
%
%      MMM2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MMM2.M with the given input arguments.
%
%      MMM2('Property','Value',...) creates a new MMM2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before mmm2_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to mmm2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help mmm2

% Last Modified by GUIDE v2.5 23-Mar-2020 23:58:07

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @mmm2_OpeningFcn, ...
                   'gui_OutputFcn',  @mmm2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before mmm2 is made visible.
function mmm2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to mmm2 (see VARARGIN)

% Choose default command line output for mmm2
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes mmm2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = mmm2_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
clear;
mido(3)
% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
clear;
mido(5)
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
clear;
mido(7);
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
clear;
mido(9)
% --- Executes on button press in pushbutton7.
function pushbutton7_Callback(hObject, eventdata, handles)
clear;
ava(3)
% --- Executes on button press in pushbutton8.
function pushbutton8_Callback(hObject, eventdata, handles)
clear;
ava(5)
% --- Executes on button press in pushbutton9.
function pushbutton9_Callback(hObject, eventdata, handles)
clear;
ava(7)
% --- Executes on button press in pushbutton10.
function pushbutton10_Callback(hObject, eventdata, handles)
clear;
ava(9)
% --- Executes on button press in pushbutton12.
function pushbutton12_Callback(hObject, eventdata, handles)
clear;
Sobel();
% --- Executes on button press in pushbutton13.
function pushbutton13_Callback(hObject, eventdata, handles)
clear;
Prewitt();
% --- Executes on button press in pushbutton14.
function pushbutton14_Callback(hObject, eventdata, handles)
clear;
Roberts();
% --- Executes on button press in pushbutton16.
function pushbutton16_Callback(hObject, eventdata, handles)
clear;
imNew=im2double(imread('锐化及边缘检测用途.jpg'));
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(imNew);
title('原图2'); 
% --- Executes on button press in pushbutton17.
function pushbutton17_Callback(hObject, eventdata, handles)
clear;
imNew=im2double(imread('noise.jpg'));
figure('toolbar','none','menubar','none');
set (gca,'position',[0.05,0.03,0.90,0.90] );
imshow(imNew);
title('原图1'); 

 
实验效果
 

 原图
 
 
第一排是中值滤波，第二排是均值滤波，
 从左到右是3X3,5X5,7X7,9X9模板大小
 
 

 从左到右分别是Sobel,Prewitt,Roberts算子