-
Notifications
You must be signed in to change notification settings - Fork 0
/
st_immidist.m
77 lines (72 loc) · 1.86 KB
/
st_immidist.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
function [D] = st_immidist(ste)
[imarray] = st_sampleimg(ste.img, ste.xy, 10);
D = zeros(length(imarray));
for k = 1:length(imarray)
for l = 1:length(imarray)
D(k, l) = mi(imarray{k}, imarray{l});
end
end
end
function n = hist2(A, B, L)
%HIST2 Calculates the joint histogram of two images or signals
%
% n=hist2(A,B,L) is the joint histogram of matrices A and B, using L
% bins for each matrix.
%
% See also MI, HIST.
% jfd, 15-11-2006, working
% 27-11-2006, memory usage reduced (sub2ind)
% 22-10-2008, added support for 1D matrices
% 01-09-2009, commented specific code for sensorimotor signals
% 24-08-2011, speed improvements by Andrew Hill
ma = min(A(:));
MA = max(A(:));
mb = min(B(:));
MB = max(B(:));
% For sensorimotor variables, in [-pi,pi]
% ma=-pi;
% MA=pi;
% mb=-pi;
% MB=pi;
% Scale and round to fit in {0,...,L-1}
A = round((A - ma)*(L - 1)/(MA - ma + eps));
B = round((B - mb)*(L - 1)/(MB - mb + eps));
n = zeros(L);
x = 0:L - 1;
for i = 0:L - 1
n(i+1, :) = histc(B(A == i), x, 1);
end
end
function I = mi(A, B, varargin)
%MI Determines the mutual information of two images or signals
%
% I=mi(A,B) Mutual information of A and B, using 256 bins for
% histograms
% I=mi(A,B,L) Mutual information of A and B, using L bins for histograms
%
% Assumption: 0*log(0)=0
%
% See also ENTROPY.
% jfd, 15-11-2006
% 01-09-2009, added case of non-double images
% 24-08-2011, speed improvements by Andrew Hill
if nargin >= 3
L = varargin{1};
else
L = 256;
end
A = double(A);
B = double(B);
na = hist(A(:), L);
na = na / sum(na);
nb = hist(B(:), L);
nb = nb / sum(nb);
n2 = hist2(A, B, L);
n2 = n2 / sum(n2(:));
I = sum(minf(n2, na'*nb));
end
% -----------------------
function y = minf(pab, papb)
I = find(papb(:) > 1e-12 & pab(:) > 1e-12); % function support
y = pab(I) .* log2(pab(I)./papb(I));
end