Matlab二维Otsu阈值分割

% 二维Otsu
clear all;  % 清除所有变量  在function里不能使用
% close all;  % 关闭所有figure   在function里不能使用

clc;

[filenamepathname]=uigetfile（{‘*.jpg;*.bmp;*.tif;*.png;*.gif‘‘All Image Files‘;‘*.*‘‘All Files‘}）; 

tic;
t1 = clock;

% 读取图像
I = imread（[pathnamefilename]）;
% I = imread（‘F:\matlab2016\bin\pic\lena.jpg‘）;
figure（）;
subplot（231）;
imshow（I）title（‘原图‘）;

numel = numel（size（I））;
% 图像灰度化
if numel>2
    I_gray = rgb2gray（I）;
else
    I_gray = I;
end
subplot（232）;
imshow（I_gray）title（‘灰度图‘）;

% 图像大小
[mn] = size（I_gray）;

% NI=filter2（fspecial（‘average‘3）I_gray）; % 第一行  第一列  最后一行  最后一列  不满足八邻域  补0计算

NI = I_gray;

% 3×3模板
for i = 2:m-1
    for j = 2:n-1
        pixsum=sum（sum（I_gray（i-1:i+1j-1:j+1）））;
        pix=round（pixsum/9）;
        NI（ij）=pix;
    end
end
subplot（233）;
imshow（NI）title（‘邻域平均灰度图‘）;

Hist=zeros（256256）;
for i = 1:m
    for j = 1:n
        pix1=I_gray（ij）;
        pix2=NI（ij）;
        Hist（pix1+1pix2+1）=Hist（pix1+1pix2+1）+1;           
    end
end
subplot（234）;
mesh（double（Hist））;
title（‘二维灰度直方图‘）;
xlim（[0 256]）;% X轴范围
ylim（[0 256]）;% Y轴范围
xlabel（‘灰度级‘）;
ylabel（‘邻域平均灰度级‘）;
zlabel（‘像素数‘）;

% 灰度级
len =