-
Notifications
You must be signed in to change notification settings - Fork 5
/
FastCMeans.m
89 lines (75 loc) · 2.51 KB
/
FastCMeans.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [C,LUT,H]=FastCMeans(im,c)
% Partition N-dimensional grayscale image into c classes using a memory
% efficient implementation of the c-means (aka k-means) clustering
% algorithm. Computational efficiency is achieved by using the histogram
% of image intensities during clustering instead of the raw image data.
%
% INPUT:
% - im : N-dimensional grayscale image in integer format.
% - c : positive integer greater than 1 specifying the number of
% clusters. c=2 is the default setting. Alternatively, c can be
% specified as a k-by-1 array of initial cluster (aka prototype)
% centroids.
%
% OUTPUT :
% - C : 1-by-k array of cluster centroids.
% - LUT : L-by-1 array that specifies the intensity-class relations,
% where L is the dynamic intensity range of the input image.
% Specifically, LUT(1) corresponds to the (class) label assigned to
% min(im(:)) and LUT(L) corresponds to the label assigned
% to max(im(:)). LUT is used as input to 'apply_LUT' function to
% create a label image.
% - H : image histogram. If I=min(im(:)):max(im(:)) are the intensities
% present in the input image, then H(i) is the number of
% pixels/voxels with intensity I(i).
%
% AUTHOR : Anton Semechko (a.semechko@gmail.com)
%
% Default input arguments
if nargin<2 || isempty(c), c=2; end
% Basic error checking
if nargin<1 || isempty(im)
error('Insufficient number of input arguments')
end
msg='Revise variable used to specify class centroids. See function documentaion for more info.';
if ~isnumeric(c) || ~isvector(c)
error(msg)
end
if numel(c)==1 && (~isnumeric(c) || round(c)~=c || c<2)
error(msg)
end
% Check image format
if isempty(strfind(class(im),'int'))
error('Input image must be specified in integer format (e.g. uint8, int16)')
end
% Intensity range
Imin=double(min(im(:)));
Imax=double(max(im(:)));
I=(Imin:Imax)';
% Compute intensity histogram
H=hist(double(im(:)),I);
H=H(:);
% Initialize cluster centroids
if numel(c)>1 % uder-defined
C=c;
c=numel(c);
else % automatic
dI=(Imax-Imin)/c;
C=Imin+dI/2:dI:Imax;
end
% Update cluster centroids
IH=I.*H;
dC=Inf;
while dC>1E-3
C0=C;
% Distance to centroids
D=abs(bsxfun(@minus,I,C));
% Classify by proximity
[~,LUT]=min(D,[],2);
for j=1:c
C(j)=sum(IH(LUT==j))/sum(H(LUT==j));
end
C=sort(C,'ascend'); % enforce natural order
% Change in centroids
dC=max(abs(C-C0));
end