meanshift图像分割实例matlab程序(支持高斯核和Epanechnikov核)

本来不想写这篇文章的，因为代码只是对CSDN上一篇博客（https://blog.csdn.net/full_speed_turbo/article/details/83181982）上的程序做了小小的改动，不敢掠人之美。不过原博客大概由于没有给出中文注释和处理结果图像而似乎并没有受到应有的关注。故而纠结了几天，最后还是决定在这里给出改动后的版本以助其发扬光大，同时希望能帮助到有需要的人。

基于meanshift的图像分割函数：

function l=meanshift(im,hs,hr,err,kernelType) 
% MEANSHSEGM Mean shift segmentation
% CMP Vision Algorithms http://visionbook.felk.cvut.cz
%  
% Segment a grayscale or color image using mean shift image segmentation
% Warning: Segmenting a medium-size image (200x 300 pixels)
% takes a couple of minutes and the time increases with
% the image and kernel size.  
%
% Usage: l = meanshsegm(im,hs,hr)
%  im   [m x n x d]  Scalar (d=1) or color (d=3) input image.
%  hs  (default 10)
%      Spatial kernel size, in pixels.
%  hr  (default 20 for gaussian kernel and default 40 for Epanechnikov kernel)
%      Range kernel size.
% err  (default 0.2)
% norm of error which determine when to stop iteration
% kernelType (default 'gaussian')
% 'gaussian' for gaussian kernel, 'Epan' for Epanechnikov kernel
% l  [m x n]  Output labeling. Each pixel position contains an integer
%  1... N corresponding to an assigned region number; N is the
%  number of regions.
  
% 设置默认参数
if ~exist('kernelType','var')
    kernelType='gaussian'
end

if ~exist('hs','var')
    hs=15;
end

if ~exist('hr','var')
    if strcmp(kernelType,'gaussian')==1
        hr=20;
    else
        hr=40;
    end
end

if ~exist('err','var')
    err=0.2;
end

% 构造空域-值域联合尺度因子
[m,n,d] = size(im);
h  = [hs hs repmat(hr,1,d)];
z  = zeros( m, n, d+2 );
im = double(im);

% 对每个像素点进行均值移位
for ix = 1:n
  for iy = 1:m
    % 构造当前点空域-值域联合特征向量     
    y = double( [ix iy reshape(im(iy,ix,:),1,d)] );
    % 构造邻域像素点的空域-值域联合特征向量(组合成矩阵)  
    xl = max(ix-hs,1);  xh = min(ix+hs,n);  nw = xh-xl+1;
    yl = max(iy-hs,1);  yh = min(iy+hs,m);  mw = yh-yl+1;
    nw = nw*mw;  iw = (0:(nw-1))';
    fw = [fix(iw/mw+xl) mod(iw,mw)+yl reshape(im(yl:yh,xl:xh,:),[],d)];

    while true
      r = (fw-repmat(y,nw,1)) ./ repmat(h,nw,1);
      rs = sum( r.*r, 2 );

      if strcmp(kernelType,'gaussian')==1
          w=exp(-rs/2);
      else
          w=1-rs;
          w(w<0)=0;
      end
      y0 = y;
      y = w'*fw/sum(w);
      if norm(y-y0)<err, break; end
    end
    z(iy,ix,:) = y;
  end % for iy
end % for ix

% 使用超像素技术根据像素对的连接关系构造连通域
% 对应原始像素点的网格点值设为1，其余设为0
% 若水平、垂直像素对的特征量各维度均小于对应尺度值，则将对应网格点的值设为1
s = ones( 2*m+1, 2*n+1, 'int8' );
s(1:2:(2*m+1),:) = zeros( m+1, 2*n+1, 'int8' );
s(:,1:2:(2*n+1)) = zeros( 2*m+1, n+1, 'int8' );

s(2:2:2*m,3:2:(2*n-1)) = all(cat(3, ...                % horizontal edges
  abs(z(:,2:end,1:2)-z(:,1:(end-1),1:2)) < hs, ...
  abs(z(:,2:end,3:end)-z(:,1:(end-1),3:end)) < hr ),3);
s(3:2:(2*m-1),2:2:2*n) = all(cat(3, ...                % vertical edges
  abs(z(2:end,:,1:2)-z(1:(end-1),:,1:2)) < hs, ...
  abs(z(2:end,:,3:end)-z(1:(end-1),:,3:end)) < hr ),3);

% 连通域标记
l = bwlabel( s, 4 );                                   % find connected regions
l = l( 2:2:2*m, 2:2:2*n );                             % extract labeling

脚本文件：

% MEANSHSEGM_DEMO Demo showing the usage of meanshsegm 
% CMP Vision algorithms http://visionbook.felk.cvut.cz
% Example

clear all;close all;clc;
ImageDir='images/';%directory containing the images
% addpath('..') ;
% cmpviapath('..') ;
if (exist('output_images')~=7)
  mkdir('output_images');
end

if 1
img=imresize(imread([ImageDir 'onion.png']),0.5) ;
l=meanshift(img,15,20,0.1,'gaussian') ;
% l=meanshift(img,15,40,0.1,'Epan') ;

figure;
imagesc(img); % title('input image');
axis image ; axis off ; 

figure; 
imagesc(label2rgb(l-1,'jet','w','shuffle')) ; 
axis image ; axis off ; 
end

F=getframe(gcf);
I=frame2im(F);
[I,map]=rgb2ind(I,256);
imwrite(I,map,'r.png')

原图：

使用gaussian核采用默认值处理结果：

使用Epanechnikov核采用默认值处理结果：