图像处理之细胞计数器

很久以前写的一段代码。医学院某个实验室有一系列电子显微镜拍摄的细胞照片，他们需要人工统计每张图片里细胞的数目，对比前后照片来分析药物对组织产生的影响。人工计数的工作量太庞大了，于是他们找我们实验室合作。

这个问题本质上就是一个图像分割的问题，把照片中的每个细胞分割开，然后求一下连通域的个数就好。

% input image ---> tophat/bottomhat filter ---> image enhancement--->binarization by a histogram based threshold
%  ---> calculate seed points ---> segmentation or plot edges of blobs ---> count connected components

整个流程比较简单，图像预处理，调整亮度和对比度，区域分割，然后是连通域处理。用Matlab实现的，代码全都贴在这里了。这是实验阶段，后来又用OpenCVS和Qt写了一个简单的包含图形界面的程序给他们使用，如果有需要的同学可以邮件联系。

function [outputImage, number_of_nuclei] = nuclei_counter(sourceImage, threshold)
% This script accept an input nuclei image taken under microscope,
% then preprocess the image with tophat and bottomhat filter, enhance its 
% contrast by multiplying gamma. Convert it to binary image with a histogram 
% based threshold, count the number of connected components and detect edges
% of nucleus with LoG filter, then plot edges. Or calculate distance transform 
% and detect seed points,use watershed algorithm for segmentation.
% 
% flow chart:
% input image ---> tophat/bottomhat filter ---> image enhancement--->binarization by a histogram based threshold
%  ---> calculate seed points ---> segmentation or plot edges of blobs ---> count connected components
%  
% refrence:
% * Yousef Al-Kofahi, et al. Improved Automatic Detection and Segmentation of Cell Nuclei in Histopathology Images. 
%   IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 4, APRIL 2010
% * Jos B.T.M, et al, The Watershed Transform: Defnitions, Algorithms and Parallelization Strategies
%   Fundamenta Informaticae 41 (2001) 187{228
% 
% Input
%     sourceImage: 8 bits grayscale image;
%     threshold: threshold [0 1] for binarization.
%     
% Output:
%     outputImage: source image with overlayed contour
%     number_of_nuclei: number of nucleis in image
%     
% If you have any ideas or improvements about this method, welcome to contact me 'zhyh8341@gmail.com'.


% display original image
image = sourceImage;
figure;imshow(image);
title('input image');

total = numel(image);


% apply top hat and bottom hat filter
se = strel('disk',30);
tophat = imtophat(image,se);
bottomhat = imbothat(image,se);
filterImage = image + (tophat - bottomhat);
se = strel('disk',15);
tophat = imtophat(filterImage,se);
bottomhat = imbothat(filterImage,se);
filterImage = filterImage + (tophat - bottomhat);

% calculate histogram of filtered image
% estimate more than 78.5% area is background (pi/4 = .785)
[counts,x] = imhist(filterImage);
ssum = cumsum(counts);
bg = .215*total;
fg = .99*total;
low = find(ssum>bg, 1, 'first');
high = find(ssum>fg, 1, 'first');
adjustedImage = imadjust(filterImage, [low/255 high/255],[0 1],1.8);


% image binarization, threshold is choosen based on experience
if(nargin < 2)
    matrix = reshape(adjustedImage,total,1);
    matrix = sort(matrix);
    threshold = graythresh(matrix(total*.5:end));
end
binarization = im2bw(adjustedImage,threshold);


% open image and then detect edge using laplacian of gaussian
se2 = strel('disk',5);
afterOpening = imopen(binarization,se2);
nsize = 5; sigma = 3;
h = fspecial('log',nsize,sigma);
afterLoG = uint8(imfilter(double(afterOpening)*255,h,'same').*(sigma^2)); 


se2 = strel('disk',5);
afterOpening = imopen(binarization,se2);
number_of_nuclei = bwconncomp(afterOpening);


% % you can either use watershed method to do segmentation
% D = -bwdist(~afterOpening);
% D(~afterOpening) = -Inf;
% L = watershed(D);


outputImage = sourceImage + afterLoG*5;
figure;imshow(outputImage);
title('output image');

一个简单的例子：

% this is just an example to show how to use this function

% read an image
I = imread('img01.JPG');

% convert it to grayscale of just take one channel
% Igray = rgb2gray(I);
Iblue = I(:,:,3);

% call the nuclei counter function
[outputImage, number_of_nuclei] = nuclei_counter(Iblue);

效果如下：