www.gusucode.com > 《matlab图像处理与界面编程宝典》秦襄培 编著，每章的MATLAB源代码程序 > 第23章/代码23-9.txt

    
he = imread('hestain.png');                                       % 读入图像
imshow(he), title('H&E image');                                  
text(size(he,2),size(he,1)+15,...                                 
     'Image courtesy of Alan Partin, Johns Hopkins University', ...
     'FontSize',7,'HorizontalAlignment','right');
cform = makecform('srgb2lab');         % 色彩空间转换
lab_he = applycform(he,cform);
ab = double(lab_he(:,:,2:3));          % 数据类型转换
nrows = size(ab,1);                    % 求矩阵尺寸
ncols = size(ab,2);                    % 求矩阵尺寸
ab = reshape(ab,nrows*ncols,2);        % 矩阵形状变换
nColors = 3;
% 重复聚类3次，以避免局部最小值
[cluster_idx cluster_center] = kmeans(ab,nColors,'distance','sqEuclidean', ...
                                      'Replicates',3);
pixel_labels = reshape(cluster_idx,nrows,ncols);   % 矩阵形状改变
imshow(pixel_labels,[]);                           % 显示图像
title('image labeled by cluster index');           % 设置图像标题
segmented_images = cell(1,3);                              % 细胞型数组
rgb_label = repmat(pixel_labels,[1 1 3]);                      % 矩阵平铺
for k = 1:nColors
    color = he;
    color(rgb_label ~= k) = 0;
    segmented_images{k} = color;
end
imshow(segmented_images{1}),                          % 显示处理后的图像
title('objects in cluster 1');                                 % 设置图像标题
imshow(segmented_images{2}),                           % 显示处理后的图像
title('objects in cluster 2');                                 % 设置图像标题
imshow(segmented_images{3}),                           % 显示处理后的图像
title('objects in cluster 3');                                 % 设置图像标题
mean_cluster_value = mean(cluster_center,2);
[tmp, idx] = sort(mean_cluster_value);
blue_cluster_num = idx(1);
L = lab_he(:,:,1);
blue_idx = find(pixel_labels == blue_cluster_num);
L_blue = L(blue_idx);
is_light_blue = im2bw(L_blue,graythresh(L_blue));         % 图像黑白转换
nuclei_labels = repmat(uint8(0),[nrows ncols]);           % 矩阵平铺
nuclei_labels(blue_idx(is_light_blue==false)) = 1;
nuclei_labels = repmat(nuclei_labels,[1 1 3]);            % 矩阵平铺
blue_nuclei = he;
blue_nuclei(nuclei_labels ~= 1) = 0;
imshow(blue_nuclei), title('blue nuclei');                 % 显示处理后的图像