目录
理论知识
代码步骤
读取图像
得到金字塔图像
提取底层特征
计算显著图
显著图综合
运行结果展示

理论知识
显著性检测 是指按照人类的视觉注意机制，判断出图像中的显著区域，并为该区域分配较高的显著值，通常认为显著区域更有可能包含目标，利用显著性检测的方法能够快速在图像中找到可能的目标区域，可以减少计算量。显著性检测主要是提取底层视觉特征来计算局部区域的差异性。



Itti模型是一种经典的自底向上基于底层显著特征计算的显著性检测模型，其实现过程大致分为高斯滤波、计算底层空间特征图和计算显著性图，结构如下图所示



篇幅原因，详细原理可以参考Itti算法
代码步骤
全部代码：github链接
读取图像
读取图像后，将其转化为double类型。
img1=im2double(img1);

得到金字塔图像
一共是九层金字塔，除去原图外，还需要八层。
%得到金字塔图像
w=fspecial('gaussian',[3 3]);
img2=imresize(imfilter(img1,w),[m/2 n/2]);
img3=imresize(imfilter(img2,w),[m/4 n/4]);
img4=imresize(imfilter(img3,w),[m/8 n/8]);
img5=imresize(imfilter(img4,w),[m/16 n/16]);
img6=imresize(imfilter(img5,w),[m/32 n/32]);
img7=imresize(imfilter(img6,w),[m/64 n/64]);
img8=imresize(imfilter(img7,w),[m/128 n/128]);
img9=imresize(imfilter(img8,w),[m/256 n/256]);

提取底层特征
分别提取亮度金字塔、色度金字塔、方向特征金字塔。

定义亮度金字塔提取函数
function I=extract_Ifeature(img)
% 输入：待提取图像
% 输出：亮度特征图
I=(img(:,:,1)+img(:,:,2)+img(:,:,3))./3;

定义色度金字塔提取函数，红绿蓝三个特征，以红色为例
function R=extract_Rfeature(img)
% 输入：待提取图像
% 输出：红色R特征
R=img(:,:,1)-((img(:,:,2)+img(:,:,3))./2);

%亮度金字塔
I1=extract_Ifeature(img1);
I2=extract_Ifeature(img2);
I3=extract_Ifeature(img3);
I4=extract_Ifeature(img4);
I5=extract_Ifeature(img5);
I6=extract_Ifeature(img6);
I7=extract_Ifeature(img7);
I8=extract_Ifeature(img8);
I9=extract_Ifeature(img9);
% figure,
% imshow(I2);


%色度金字塔
%红色特征
R1=extract_Rfeature(img1);
R2=extract_Rfeature(img2);
R3=extract_Rfeature(img3);
R4=extract_Rfeature(img4);
R5=extract_Rfeature(img5);
R6=extract_Rfeature(img6);
R7=extract_Rfeature(img7);
R8=extract_Rfeature(img8);
R9=extract_Rfeature(img9);
% figure,
% imshow(R1);

%绿色特征
G1=extract_Gfeature(img1);
G2=extract_Gfeature(img2);
G3=extract_Gfeature(img3);
G4=extract_Gfeature(img4);
G5=extract_Gfeature(img5);
G6=extract_Gfeature(img6);
G7=extract_Gfeature(img7);
G8=extract_Gfeature(img8);
G9=extract_Gfeature(img9);
% figure,
% imshow(G1);

%蓝色特征
B1=extract_Bfeature(img1);
B2=extract_Bfeature(img2);
B3=extract_Bfeature(img3);
B4=extract_Bfeature(img4);
B5=extract_Bfeature(img5);
B6=extract_Bfeature(img6);
B7=extract_Bfeature(img7);
B8=extract_Bfeature(img8);
B9=extract_Bfeature(img9);
% figure,
% imshow(B1);

%黄色特征
Y1=extract_Yfeature(img1);
Y2=extract_Yfeature(img2);
Y3=extract_Yfeature(img3);
Y4=extract_Yfeature(img4);
Y5=extract_Yfeature(img5);
Y6=extract_Yfeature(img6);
Y7=extract_Yfeature(img7);
Y8=extract_Yfeature(img8);
Y9=extract_Yfeature(img9);
% figure,
%imshow(Y1);

提取方向特征，使用Matlab自带的gabor滤波器
%方向特征
%使用内置gabor滤波器
gaborArray = gabor(3,[0 45 90 135]);%滤波器数组，波长为10，四个方向
img1_gaborMag=imgaborfilt(rgb2gray(img1),gaborArray);%结果包含四个方向的特征图，在第三维度
img2_gaborMag=imgaborfilt(rgb2gray(img2),gaborArray);
img3_gaborMag=imgaborfilt(rgb2gray(img3),gaborArray);
img4_gaborMag=imgaborfilt(rgb2gray(img4),gaborArray);
img5_gaborMag=imgaborfilt(rgb2gray(img5),gaborArray);
img6_gaborMag=imgaborfilt(rgb2gray(img6),gaborArray);
img7_gaborMag=imgaborfilt(rgb2gray(img7),gaborArray);
img8_gaborMag=imgaborfilt(rgb2gray(img8),gaborArray);
img9_gaborMag=imgaborfilt(rgb2gray(img9),gaborArray);
figure('Name',"原图的 Gabor 特征图（四个方向）")
for i=1:4
    subplot(2,2,i),imshow(img1_gaborMag(:,:,i),[]);
    theta = gaborArray(i).Orientation;
    lambda = gaborArray(i).Wavelength;
    title(sprintf('方向=%d, 波长=%d',theta,lambda));
end

计算显著图
根据公式计算显著图，这里分辨率不同的特征图合并时以函数中最大的分辨率为标准。

定义亮度特征取差函数
function I_cs=Ifeature_diff(I1,I2)
% 输入：I1为高分辨率图像，I2为低分辨率图像
% 输出：亮度特征取差结果
[m,n]=size(I1);
I2=imresize(I2,[m,n],'nearest');
I_cs=abs(double(I1)-double(I2));

定义色度特征取差函数
function RGBY=RGBfeature_diff(Rc,Gc,Rs,Gs)
%输入：高分辨率图像红色、绿色特征图，低分辨率图像红色、绿色特征图
%输出：红绿特征取差结果
[m,n]=size(Rc);
Rs=imresize(Rs,[m,n],'nearest');
Gs=imresize(Gs,[m,n],'nearest');
RGBY=abs((double(Rc)-double(Gc))-(double(Gs)-double(Rs)));

定义方向特征取差函数
function ocs=directionfeature_diff(oc,os)
[m,n]=size(oc);
os=imresize(os,[m,n],'nearest');
ocs=abs(oc-os);

%%
%%计算显著图,分辨率以函数中最大的分辨率为标准
%亮度特征取差
I_2_5=double(Ifeature_diff(I2,I5));
I_2_6=double(Ifeature_diff(I2,I6));
I_3_6=double(Ifeature_diff(I3,I6));
I_3_7=double(Ifeature_diff(I3,I7));
I_4_7=double(Ifeature_diff(I4,I7));
I_4_8=double(Ifeature_diff(I4,I8));

%色度特征取差
%RG
RG_2_5=double(RGBfeature_diff(R2,G2,R5,G5));
RG_2_6=double(RGBfeature_diff(R2,G2,R6,G6));
RG_3_6=double(RGBfeature_diff(R3,G3,R6,G6));
RG_3_7=double(RGBfeature_diff(R3,G3,R7,G7));
RG_4_7=double(RGBfeature_diff(R4,G4,R7,G7));
RG_4_8=double(RGBfeature_diff(R4,G4,R8,G8));

%BY
BY_2_5=double(RGBfeature_diff(B2,B2,Y5,Y5));
BY_2_6=double(RGBfeature_diff(B2,B2,Y6,Y6));
BY_3_6=double(RGBfeature_diff(B3,B3,Y6,Y6));
BY_3_7=double(RGBfeature_diff(B3,B3,Y7,Y7));
BY_4_7=double(RGBfeature_diff(B4,B4,Y7,Y7));
BY_4_8=double(RGBfeature_diff(B4,B4,Y8,Y8));

%方向特征取差
%0度
O_2_5_0=double(directionfeature_diff(img2_gaborMag(:,:,1),img5_gaborMag(:,:,1)));
O_2_6_0=double(directionfeature_diff(img2_gaborMag(:,:,1),img6_gaborMag(:,:,1)));
O_3_6_0=double(directionfeature_diff(img3_gaborMag(:,:,1),img6_gaborMag(:,:,1)));
O_3_7_0=double(directionfeature_diff(img3_gaborMag(:,:,1),img7_gaborMag(:,:,1)));
O_4_7_0=double(directionfeature_diff(img4_gaborMag(:,:,1),img7_gaborMag(:,:,1)));
O_4_8_0=double(directionfeature_diff(img4_gaborMag(:,:,1),img8_gaborMag(:,:,1)));

%45度
O_2_5_45=double(directionfeature_diff(img2_gaborMag(:,:,2),img5_gaborMag(:,:,2)));
O_2_6_45=double(directionfeature_diff(img2_gaborMag(:,:,2),img6_gaborMag(:,:,2)));
O_3_6_45=double(directionfeature_diff(img3_gaborMag(:,:,2),img6_gaborMag(:,:,2)));
O_3_7_45=double(directionfeature_diff(img3_gaborMag(:,:,2),img7_gaborMag(:,:,2)));
O_4_7_45=double(directionfeature_diff(img4_gaborMag(:,:,2),img7_gaborMag(:,:,2)));
O_4_8_45=double(directionfeature_diff(img4_gaborMag(:,:,2),img8_gaborMag(:,:,2)));

%90度
O_2_5_90=double(directionfeature_diff(img2_gaborMag(:,:,3),img5_gaborMag(:,:,3)));
O_2_6_90=double(directionfeature_diff(img2_gaborMag(:,:,3),img6_gaborMag(:,:,3)));
O_3_6_90=double(directionfeature_diff(img3_gaborMag(:,:,3),img6_gaborMag(:,:,3)));
O_3_7_90=double(directionfeature_diff(img3_gaborMag(:,:,3),img7_gaborMag(:,:,3)));
O_4_7_90=double(directionfeature_diff(img4_gaborMag(:,:,3),img7_gaborMag(:,:,3)));
O_4_8_90=double(directionfeature_diff(img4_gaborMag(:,:,3),img8_gaborMag(:,:,3)));

%135度
O_2_5_135=double(directionfeature_diff(img2_gaborMag(:,:,4),img5_gaborMag(:,:,4)));
O_2_6_135=double(directionfeature_diff(img2_gaborMag(:,:,4),img6_gaborMag(:,:,4)));
O_3_6_135=double(directionfeature_diff(img3_gaborMag(:,:,4),img6_gaborMag(:,:,4)));
O_3_7_135=double(directionfeature_diff(img3_gaborMag(:,:,4),img7_gaborMag(:,:,4)));
O_4_7_135=double(directionfeature_diff(img4_gaborMag(:,:,4),img7_gaborMag(:,:,4)));
O_4_8_135=double(directionfeature_diff(img4_gaborMag(:,:,4),img8_gaborMag(:,:,4)));

显著图综合
定义归一化函数，归一化函数很关键
function img1 = normalizeImg(img)
%输入：待归一化特征图
%输出：归一化特征图
[r,c]=size(img);
m = 0;%动态存储局部最小值
M=1;%全局最大值

img = mat2gray(img);%全局最大值M归一化,并且数值类型转化为double

%寻找局部最大值
step=16;%搜索步长
for i=1:step:r-step
    for j=1:step:c-step
        tempImg=img(i:i+step-1,j:j+step-1);%取块
        m=[m,max(max(tempImg))];
    end
end

%取平均值
[~,cc]=size(m);
sum=0;
for i=1:cc
    sum=sum+m(i);
end
m_=sum/cc;
img1=img+((M-m_)^2);

%%综合显著图
%归一化图像
I_2_5=normalizeImg(I_2_5);
I_2_6=normalizeImg(I_2_6);
I_3_6=normalizeImg(I_3_6);
I_3_7=normalizeImg(I_3_7);
I_4_7=normalizeImg(I_4_7);
I_4_8=normalizeImg(I_4_8);

RG_2_5=normalizeImg(RG_2_5);
RG_2_6=normalizeImg(RG_2_6);
RG_3_6=normalizeImg(RG_3_6);
RG_3_7=normalizeImg(RG_3_7);
RG_4_7=normalizeImg(RG_4_7);
RG_4_8=normalizeImg(RG_4_8);

BY_2_5=normalizeImg(BY_2_5);
BY_2_6=normalizeImg(BY_2_6);
BY_3_6=normalizeImg(BY_3_6);
BY_3_7=normalizeImg(BY_3_7);
BY_4_7=normalizeImg(BY_4_7);
BY_4_8=normalizeImg(BY_4_8);

O_2_5_0=normalizeImg(O_2_5_0);
O_2_6_0=normalizeImg(O_2_6_0);
O_3_6_0=normalizeImg(O_3_6_0);
O_3_7_0=normalizeImg(O_3_7_0);
O_4_7_0=normalizeImg(O_4_7_0);
O_4_8_0=normalizeImg(O_4_8_0);

O_2_5_45=normalizeImg(O_2_5_45);
O_2_6_45=normalizeImg(O_2_6_45);
O_3_6_45=normalizeImg(O_3_6_45);
O_3_7_45=normalizeImg(O_3_7_45);
O_4_7_45=normalizeImg(O_4_7_45);
O_4_8_45=normalizeImg(O_4_8_45);

O_2_5_90=normalizeImg(O_2_5_90);
O_2_6_90=normalizeImg(O_2_6_90);
O_3_6_90=normalizeImg(O_3_6_90);
O_3_7_90=normalizeImg(O_3_7_90);
O_4_7_90=normalizeImg(O_4_7_90);
O_4_8_90=normalizeImg(O_4_8_90);

O_2_5_135=normalizeImg(O_2_5_135);
O_2_6_135=normalizeImg(O_2_6_135);
O_3_6_135=normalizeImg(O_3_6_135);
O_3_7_135=normalizeImg(O_3_7_135);
O_4_7_135=normalizeImg(O_4_7_135);
O_4_8_135=normalizeImg(O_4_8_135);

%综合计算
%亮度
I_c2=I_2_5+I_2_6;
I_c3=I_3_7+I_3_6;
I_c4=I_4_7+I_4_8;

I_c3=imresize(I_c3,size(I_c2),'nearest');
I_c4=imresize(I_c4,size(I_c2),'nearest');

I=(I_c2+I_c3+I_c4);

%色度
c2=RG_2_5+BY_2_5+RG_2_6+BY_2_6;
c3=RG_3_6+BY_3_6+RG_3_7+BY_3_7;
c4=RG_4_7+BY_4_7+RG_4_8+BY_4_8;

c3=imresize(c3,size(c2),'nearest');
c4=imresize(c4,size(c2),'nearest');

C=(c2+c3+c4);

%方向
O2=O_2_5_0+O_2_5_45+O_2_5_90+O_2_5_135+O_2_6_0+O_2_6_45+O_2_6_90+O_2_6_135;
O3=O_3_6_0+O_3_6_45+O_3_6_90+O_3_6_135+O_3_7_0+O_3_7_45+O_3_7_90+O_3_7_135;
O4=O_4_7_0+O_4_7_45+O_4_7_90+O_4_7_135+O_4_8_0+O_4_8_45+O_4_8_90+O_4_8_135;

O3=imresize(O3,size(O2),'nearest');
O4=imresize(O4,size(O2),'nearest');


O=normalizeImg(O2+O3+O4);

S=(double(I)+double(O)+double(C))./3;

figure,
imshow(S,[]);

运行结果展示
原图



金字塔图像



方向特征



显著图

主要对图像边缘进行检测，然后输出图像中某一个连通域的边缘坐标。
测试图像
从左到右从上到下分别为灰度处理、二值化、扣除小连通域、边缘识别
得到各个连通域的标号
绘制连通域2的坐标（坐标原点为原图左上角）
从第六个图形边缘提取结果中可以看到图形内部包含多个子连通域，目前还没有想到什么好办法能够把这个图形进行很好的处理，后续有时间再改进。
程序大部分内容用到了matlab内部包含的命令。
clc
clear all;
close all;
I_orin= imread('test.jpg');
m=size(I_orin,1);
n=size(I_orin,2);
I=rgb2gray(I_orin);
t=graythresh(I);
I_bw=imbinarize(I,t);
I_open=bwareaopen(~I_bw,1000);
I_edge=edge(I_open,'canny',0.09);
[B,L]=bwboundaries(I_edge,'noholes');
I_compare=zeros(m,n,3);
for i=1:3
I_compare(:,:,i)=I;
end
for i=1:m
    for j=1:n
        if I_edge(i,j)==1
            I_compare(i,j,1)=255;
            I_compare(i,j,2)=0;
            I_compare(i,j,3)=0;
        end
    end
end
figure(1)
subplot(221)
imshow(I)
subplot(222)
imshow(I_bw)
subplot(223)
imshow(I_open)
subplot(224)
imshow(I_edge)
figure(2)
imshow(I_edge); hold on;
colors=['b' 'g' 'r' 'c' 'm' 'y'];
for k=1:length(B)
  boundary = B{k};
  cidx = mod(k,length(colors))+1;
  plot(boundary(:,2), boundary(:,1),...
       colors(cidx),'LineWidth',2);

  %randomize text position for better visibility
  rndRow = ceil(length(boundary)/(mod(rand*k,7)+1));
  col = boundary(rndRow,2); row = boundary(rndRow,1);
  h = text(col+1, row-1, num2str(L(row,col)));
  set(h,'Color',colors(cidx),'FontSize',14,'FontWeight','bold');
end
hold off;
figure(3)
%连通域边界坐标
location=B{2};
plot(location(:,2),location(:,1),'bd-')

1.《一种简单的图像显著性计算模型 》，提供Matlab & opencv代码
http://www.cnblogs.com/ccbb/archive/2011/05/19/2051442.html


2.《图像特征：三个图像显著性区域特征提取方法》，提供三种算法的Matlab代码 
http://www.cnblogs.com/skyseraph/archive/2011/08/11/2135291.html


3.《国外的一个图像显著区域检测代码及其效果图 saliency region detection 》，提供C++源代码
http://blog.csdn.net/onezeros/article/details/6299745


4.《Global contrast based salient region detection》，提供C++源代码，并提供了其他算法的链接！！！
http://mmcheng.net/salobj/





LOG算子和CANNY算子边缘提取


边缘的含义：

在数字图像中，边缘是指图像局部变化最显著的部分，边缘主要存在于目标与目标，目标与背景之间，是图像局部特性的不连续性，如灰度的突变、纹理结构的突变、颜色的突变等。尽管图像的边缘点产生的原因各不相同，但他们都是图形上灰度不连续或灰度急剧变化的点，图像边缘分为阶跃状、斜坡状和屋顶状。


一般图像边缘检测方法主要有如下四个步骤：

	图像滤波：传统边缘检测算法主要是基于图像强度的一阶和二阶导数，但导数的计算对噪声很敏感，因此必须使用滤波器来改善与噪声有关的边缘检测器的性能。需要指出的是，大多数滤波器在降低噪声的同时也造成了边缘强度的损失，因此，在增强边缘和降低噪声之间需要一个折衷的选择。

	图像增强：增强边缘的基础是确定图像各点邻域强度的变化值。增强算法可以将邻域(或局部)强度值有显著变化的点突显出来。边缘增强一般是通过计算梯度的幅值来完成的。

	图像检测：在图像中有许多点的梯度幅值比较大，而这些点在特定的应用领域中并不都是边缘，所以应该用某种方法来确定哪些点是边缘点。最简单的边缘检测判断依据是梯度幅值。

	图像定位：如果某一应用场合要求确定边缘位置，则边缘的位置可在子像素分辨率上来估计，边缘的方位也可以被估计出来。


LOG边缘检测算子

在20世纪70年代，Marr理论根据神经生理学实验得出了以下结论：物体的边界是将亮度图像与其解释连接起来的最重要线索。边缘检测技术在当时是基于很小邻域的卷积，只对特殊图像效果好。这些边缘检测子的主要缺点是它们依赖物体的大小且对噪声敏感。

基于二阶导数过零点的边缘检测技术探究了阶跃边缘对应于图像函数陡峭的变化这一事实。图像函数的一阶导数在对应于图像边缘的位置上应该取得极值，因此二阶导数在同一位置应该为0；而寻找过零点位置比起极值来得更容易和更准确。关键的问题是如何稳定地计算二阶导数，一种可能性是首先平滑图像（减小噪声），再计算二阶导数。在选择平滑滤波器时，需要满足两个标准：

	滤波器应该是平滑的且在邻域中大致上是有限带宽的，以便减少会导致函数变化的可能频率数。

	空间定位的约束要求滤波器的响应来自于图像中邻近的点。这两个标准矛盾的，但是可以通过使用高斯分布同时得到优化。在实践中，需要准确地考虑优化的含义。

选择Laplacian算子的基础有两个基本概念。第一，算子的高斯部分会模糊图像，从而在尺寸上将结构的灰度(包括噪声)降低到远小于σ的程度，而且高斯函数能在空间和频率两个域平滑图像，因而在原图像中引入不存在的人为干扰(如振铃)的可能性很小。第二，拉普拉斯有各向同性(旋转不变)的重要优点，符合人的视觉系统特性，而且对任何模板方向的灰度变化有相等的响应，从而避免了使用多个模板去计算图像中任何点处的最强响应。

LOG算法步骤如下：

	取样得到的高斯低通滤波器对输入图像滤波。

	计算第一步得到图像的拉普拉斯。

	找到步骤2所得图像的零交叉。


Canny边缘检测算子介绍

1986年，JOHN CANNY 提出一个很好的边缘检测算法，被称为Canny边缘检测器，

Canny方法基于三个基本的目标：

	检测标准: 低错误率。所有的边缘都应该被找到，并且没有假边缘。

	行为标准: 边缘点应被很好的定位。测到的边缘点与真实边缘的中心之间距离应该最小。

	单响应标准: 单一的边缘点效应。对于真实的边缘点，检测器应该仅仅返回一个点。也就是真实边缘的局部最大数应该是最小的。意味着仅存在单一边缘点的位置，检测器不应指出多个边缘像素。

Canny边缘检测算法由一下基本步骤组成：

	用一个高斯滤波器平滑输入图像。

	计算梯度幅值图像和角度图像。

	对梯度幅值图像应用非最大抑制。

	用双阈值处理和连接分析来检测并连接边缘。


比较分析

	相同尺度下的LOG算子总是能比Canny算子检测出更多的细节，Canny算子能对真正比较显著地边缘给出检测。LOG算子容易受尺度的影响，不同尺度下的边缘点要用不同尺度的LOG算子检测，Canny 算子受尺度的影响不太明显，不同尺度下，边缘点的位置都有偏差，但几乎相同；

	从对噪声的敏感程度来看，LOG边缘检测子是采用二阶导数过零点的检测方法，故对噪声更敏感一些，因此从抑制噪声方面来讲，Canny边缘检测子不容易受到噪声的干扰，而相同尺度下LOG算子却容易受到噪声的干扰，抑制噪声的能力要弱一些。LOG算子对噪声的抑制能力随着尺度的增加而增加，相同尺度下的Canny算子比LOG算子的抗噪声能力强，而LOG算子比Canny算子的边缘点准确；

	对于弱边缘，由于Canny边缘检测子采用两种不同的阈值分别检测强边缘和弱边缘，并且当弱边缘和强边缘相连时，才将弱边缘包含在输出图像中，故而Canny算子更能检测出真正的弱边缘，但是Canny边缘检测子检测出的边缘的位置会有一定范围的误差，

	对于假边缘，LOG边缘检测子相对比较容易受到噪声干扰，会检测出更多的细节，也容易检测出一些由于噪声引起的假边缘，但是LOG边缘检测子对边缘位置的检测还是很准确的。


参考文献

[1]	D.Marr and E. Hildreth, Theory of Edge Detection, Proc. R. Soc. Lond. B207:187-217.

[2]	John Canny, A Computational Approach to Edge Detection, IEEE Trans. PAMI, 8(6):679-698.

[3]	James J. Clark, Authenticating Edges Produced by Zero-crossing Algorithms, IEEE T. PAMI, 11(1), 1989, pp.43-57.


附录：实验代码
code：独立加性高斯白噪声
function [distImg,endtImg] = add_noise(origImg,u)
origImg = im2double(origImg);
distImg = imnoise(origImg, 'gaussian', 0, u);
endtImg = distImg.*255;
end

code：对图像进行高斯滤波
function g=my_gaussian(src,n,k)  
n1=(n+1)/2;                                         
[m,l]=size(src);  
b=zeros(n,n);                                      
I=zeros(m,l);  
g=zeros(m,l);  
img=zeros(m,l);  
[Ay, Ax ,dim ] = size(src);
if dim>1
    src = rgb2gray(src);
end
for i=1:n  
    for j=1:n  
          b(i,j)=(exp(-((i-n1)^2 +(j-n1)^2)/(2*k)))/(2*pi*k);  
    end  
end  
b = b/sum(b(:));        
I=double(src);  
img=conv2(I,b,'same');              
g=uint8(img);    
end

code：LOG算子实现图像处理以及结果展示
function [thrr] =log_img(noise,n,k,thr)
[filename, pathname] = uigetfile({'*.jpg'; '*.bmp'; '*.gif'}, '选择图片');
if filename == 0
    return;
end
imgsrc = imread([pathname, filename]);
[y, x, dim] = size(imgsrc);  %转换为灰度图像
if dim>1
    imgsrc = rgb2gray(imgsrc);
end
if (x~=256 && y~=256)  %转换为8bit 256*256图像
    imgsrc = im2uint8(imgsrc);
    imgsrc=imresize(imgsrc,[256,256]);
end
[showing_noise, img_noise]= add_noise(imgsrc,noise);
img= my_gaussian(img_noise,n,k); 
[img_log,thrr]=edge(img,'log',thr);
figure(1)    
        imshow(imgsrc);%原图
figure(2)  
        imshow(showing_noise);
figure(3)  
        imshow(img);%高斯滤波后
figure(4)  
        imshow(img_log);
end  

code：CANNY算子实现
function [ m, theta, sector, canny1,  canny2, bin] = canny_step( src,  lowTh)
[Ay, Ax, dim ] = size(src);
if dim>1
    src = rgb2gray(src);
end
src = double(src);
m = zeros(Ay, Ax); 
theta = zeros(Ay, Ax);
sector = zeros(Ay, Ax);
canny1 = zeros(Ay, Ax);
canny2 = zeros(Ay, Ax);
bin = zeros(Ay, Ax);
 
%梯度计算
for y = 1:(Ay-1)
    for x = 1:(Ax-1)
        gx =  src(y, x) + src(y+1, x) - src(y, x+1)  - src(y+1, x+1);
        gy = -src(y, x) + src(y+1, x) - src(y, x+1) + src(y+1, x+1);
        m(y,x) = (gx^2+gy^2)^0.5 ;
        theta(y,x) = atand(gx/gy)  ;
        tem = theta(y,x); 
        if (tem<67.5)&&(tem>22.5)
            sector(y,x) =  0;    
        elseif (tem<22.5)&&(tem>-22.5)
            sector(y,x) =  3;    
        elseif (tem<-22.5)&&(tem>-67.5)
            sector(y,x) =   2;    
        else
            sector(y,x) =   1;    
        end
            
    end    
end

%非极大值抑制
for y = 2:(Ay-1)
    for x = 2:(Ax-1)        
        if 0 == sector(y,x) %右上 - 左下
            if ( m(y,x)>m(y-1,x+1) )&&( m(y,x)>m(y+1,x-1)  )
                canny1(y,x) = m(y,x);
            else
                canny1(y,x) = 0;
            end
        elseif 1 == sector(y,x) %竖直方向
            if ( m(y,x)>m(y-1,x) )&&( m(y,x)>m(y+1,x)  )
                canny1(y,x) = m(y,x);
            else
                canny1(y,x) = 0;
            end
        elseif 2 == sector(y,x) %左上 - 右下
            if ( m(y,x)>m(y-1,x-1) )&&( m(y,x)>m(y+1,x+1)  )
                canny1(y,x) = m(y,x);
            else
                canny1(y,x) = 0;
            end
        elseif 3 == sector(y,x) %横方向
            if ( m(y,x)>m(y,x+1) )&&( m(y,x)>m(y,x-1)  )
                canny1(y,x) = m(y,x);
            else
                canny1(y,x) = 0;
            end
        end        
    end
end

%双阈值检测
ratio = 2;
for y = 2:(Ay-1)
    for x = 2:(Ax-1)        
        if canny1(y,x)<lowTh 
            canny2(y,x) = 0;
            bin(y,x) = 0;
            continue;
        elseif canny1(y,x)>ratio*lowTh 
            canny2(y,x) = canny1(y,x);
            bin(y,x) = 1;
            continue;
        else 
            tem =[canny1(y-1,x-1), canny1(y-1,x), canny1(y-1,x+1);
                       canny1(y,x-1),    canny1(y,x),   canny1(y,x+1);
                       canny1(y+1,x-1), canny1(y+1,x), canny1(y+1,x+1)];
            temMax = max(tem);
            if temMax(1) > ratio*lowTh
                canny2(y,x) = temMax(1);
                bin(y,x) = 1;
                continue;
            else
                canny2(y,x) = 0;
                bin(y,x) = 0;
                continue;
            end
        end
    end
end
end

code：CANNY图像处理及结果展示
function  canny_plotimg(noise,n,thr,k)
[filename, pathname] = uigetfile({'*.jpg'; '*.bmp'; '*.gif'}, '选择图片');
if filename == 0
    return;
end
imgsrc = imread([pathname, filename]);
[y, x, dim] = size(imgsrc);
if dim>1
    imgsrc = rgb2gray(imgsrc);
end
if (x~=256 && y~=256)   %转换为8bit 256*256图像
    imgsrc = im2uint8(imgsrc);
    imgsrc=imresize(imgsrc,[256,256]);
end 
[showing_noise, img_noise]= add_noise(imgsrc,noise);
img= my_gaussian(img_noise,n,k); 
[m theta sector canny1  canny2 bin] = canny1step(img, thr);
figure(1)    
        imshow(imgsrc);%原图
figure(2)  
        imshow(showing_noise);
figure(3)  
        imshow(img);%高斯滤波后
figure(4)  
        imshow(uint8(m));%导数
figure(5)  
        imshow(uint8(canny1));%非极大值抑制
figure(6)  
        imshow(uint8(canny2));%双阈值