第一步：将彩色图像从RGB转化到lab色彩空间
close all; clear all; clc; I_rgb = imread('a.jpg'); figure(1); %将彩色图像从RGB转化到lab色彩空间 C = makecform('srgb2lab');%设置转换格式 I_lab = applycform(I_rgb, C);
第二步：取出lab空间的a分量和b分量
ab = double(I_lab(:,:,2:3)); nrows = size(ab,1); ncols = size(ab,2); ab = reshape(ab,nrows*ncols,2);
第三步：设定分割的区域个数进行分割
nColors = 2; %分割的区域个数,可进行修改 [cluster_center,cluster_idx,mindist,q2,quality] = kmeanss(ab,nColors,2); pixel_labels = reshape(cluster_idx,nrows,ncols); figure(1); imshow(pixel_labels,[],'border','tight','InitialMagnification','fit'); title('聚类结果');
其中的kmeanss如下：
function [centers,mincenter,mindist,q2,quality] = kmeanss(data,initcenters,method) tic if nargin < 3 method = 2; end [n,dim] = size(data); if max(size(initcenters)) == 1 k = initcenters; [centers, mincenter, mindist, lower, computed] = anchors(mean(data),k,data); total = computed; skipestep = 1; else centers = initcenters; mincenter = zeros(n,1); total = 0; skipestep = 0; [k,dim2] = size(centers); if dim ~= dim2 error('dim(data) ~= dim(centers)'); end; end nchanged = n; iteration = 0; oldmincenter = zeros(n,1); while nchanged > 0 % do one E step, then one M step computed = 0; if method == 0 & ~skipestep for i = 1:n for j = 1:k distmat(i,j) = calcdist(data(i,:),centers(j,:)); end end [mindist,mincenter] = min(distmat,[],2); computed = k*n; elseif (method == 1 | (method == 2 & iteration == 0)) & ~skipestep mindist = Inf*ones(n,1); lower = zeros(n,k); for j = 1:k jdist = calcdist(data,centers(j,:)); lower(:,j) = jdist; track = find(jdist < mindist); mindist(track) = jdist(track); mincenter(track) = j; end computed = k*n; elseif method == 2 & ~skipestep computed = 0; % for each center, nndist is half the distance to the nearest center % if d(x,center) < nndist then x cannot belong to any other center % mindist is an upper bound on the distance of each point to its nearest center nndist = min(centdist,[],2); % the following usually is not faster % ldist = min(lower,[],2); % mobile = find(mindist > max(nndist(mincenter),ldist)); mobile = find(mindist > nndist(mincenter)); % recompute distances for point i and center j % only if j can possibly be the new nearest center % for speed, the first check has been optimized by modifying centdist % swapping the order of the checks is slower for data with natural clusters mdm = mindist(mobile); mcm = mincenter(mobile); for j = 1:k % the following is incorrect: for j = unique(mcm)' track = find(mdm > centdist(mcm,j)); if isempty(track) continue; end alt = find(mdm(track) > lower(mobile(track),j)); if isempty(alt) continue; end track1 = mobile(track(alt)); % calculate exact distances to the mincenter % recalculate separately for each jj to avoid copying too much of data % redo may be empty, but we don't need to check this redo = find(~recalculated(track1)); redo = track1(redo); c = mincenter(redo); computed = computed + size(redo,1); for jj = unique(c)' rp = redo(find(c == jj)); udist = calcdist(data(rp,:),centers(jj,:)); lower(rp,jj) = udist; mindist(rp) = udist; end recalculated(redo) = 1; track2 = find(mindist(track1) > centdist(mincenter(track1),j)); track1 = track1(track2); if isempty(track1) continue; end % calculate exact distances to center j track4 = find(lower(track1,j) < mindist(track1)); if isempty(track4) continue; end track5 = track1(track4); jdist = calcdist(data(track5,:),centers(j,:)); computed = computed + size(track5,1); lower(track5,j) = jdist; % find which points really are assigned to center j track2 = find(jdist < mindist(track5)); track3 = track5(track2); mindist(track3) = jdist(track2); mincenter(track3) = j; end % for j=1:k end % if method oldcenters = centers; diff = find(mincenter ~= oldmincenter); diffj = unique([mincenter(diff);oldmincenter(diff)])'; diffj = diffj(find(diffj > 0)); if size(diff,1) < n/3 & iteration > 0 for j = diffj plus = find(mincenter(diff) == j); minus = find(oldmincenter(diff) == j); oldpop = pop(j); pop(j) = pop(j) + size(plus,1) - size(minus,1); if pop(j) == 0 continue; end centers(j,:) = (centers(j,:)*oldpop + sum(data(diff(plus),:),1) - sum(data(diff(minus),:),1))/pop(j); end else for j = diffj track = find(mincenter == j); pop(j) = size(track,1); if pop(j) == 0 continue; end % it's correct to have mean(data(track,:),1) but this can make answer worse! centers(j,:) = mean(data(track,:),1); end end if method == 2 for j = diffj offset = calcdist(centers(j,:),oldcenters(j,:)); computed = computed + 1; if offset == 0 continue; end track = find(mincenter == j); mindist(track) = mindist(track) + offset; lower(:,j) = max(lower(:,j) - offset,0); end recalculated = zeros(n,1); realdist = alldist(centers); centdist = 0.5*realdist + diag(Inf*ones(k,1)); computed = computed + k + k*(k-1)/2; end nchanged = size(diff,1) + skipestep; iteration = iteration+1; skipestep = 0; oldmincenter = mincenter; [iteration toc nchanged computed] total = total + computed; end % while nchanged > 0 udist = calcdist(data,centers(mincenter,:)); quality = mean(udist); q2 = mean(udist.^2); [iteration toc quality q2 total]
在kmeanss 中的calcdist如下所示：
function distances = calcdist(data,center) % input: vector of data points, single center or multiple centers % output: vector of distances [n,dim] = size(data); [n2,dim2] = size(center); % Using repmat is slower than using ones(n,1) % delta = data - repmat(center,n,1); % delta = data - center(ones(n,1),:); % The following is fastest: not duplicating the center at all if n2 == 1 distances = sum(data.^2, 2) - 2*data*center' + center*center'; elseif n2 == n distances = sum( (data - center).^2 ,2); else error('bad number of centers'); end % Euclidean 2-norm distance: distances = sqrt(distances); % Inf-norm distance: % distances = max(abs(distances),[],2);
在kmeanss 中的anchors如下所示：
function [centers, mincenter, mindist, lower, computed] = anchors(firstcenter,k,data) % choose k centers by the furthest-first method [n,dim] = size(data); centers = zeros(k,dim); lower = zeros(n,k); mindist = Inf*ones(n,1); mincenter = ones(n,1); computed = 0; centdist = zeros(k,k); for j = 1:k if j == 1 newcenter = firstcenter; else [maxradius,i] = max(mindist); newcenter = data(i,:); end centers(j,:) = newcenter; centdist(1:j-1,j) = calcdist(centers(1:j-1,:),newcenter); centdist(j,1:j-1) = centdist(1:j-1,j)'; computed = computed + j-1; inplay = find(mindist > centdist(mincenter,j)/2); newdist = calcdist(data(inplay,:),newcenter); computed = computed + size(inplay,1); lower(inplay,j) = newdist; move = find(newdist < mindist(inplay)); shift = inplay(move); mincenter(shift) = j; mindist(shift) = newdist(move); end
在kmeanss 中的alldist如下所示：
function centdist = alldist(centers) % output: matrix of all pairwise distances % input: data points (centers) k = size(centers,1); centdist = zeros(k,k); for j = 1:k centdist(1:j-1,j) = calcdist(centers(1:j-1,:),centers(j,:)); end centdist = centdist+centdist';
第四步：显示分割后的各个区域
dst = cell(1,nColors); rgb_label = repmat(pixel_labels,[1 1 3]); for k = 1:nColors color = I_rgb; color(rgb_label ~= k) = 0; dst{k} = color; end for i=1:nColors figure(i+2); imshow(dst{i}); title('分割结果'); end
