/
compute_histogram.m
98 lines (85 loc) · 2.08 KB
/
compute_histogram.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
89
90
91
92
93
94
95
96
97
98
function [h,x] = compute_histogram(M, options)
% compute_histogram - compute the (symmetric) histogram of a vector after threshold.
%
% [h,x] = compute_histogram(M,options);
%
% If M contains only entry in 1...p, where p=max(M), then h is of size p
% h(k) = #{i \ M(i)=k} / length(M)
% If M contains only entry in -p...p, where p=max(abs(M)), then h is of size 2*p+1
% h(k) = #{i \ M(i)=k-p-1} / length(M)
% (you can also force this by setting options.symmetrize=1).
%
% For real valued signals, you can add thresholding by setting options.T
%
% Copyright (c) 2006 Gabriel Peyré
if nargin>1 && ~isstruct(options)
error('Valid syntax is "h = compute_histogram(M,options);".');
end
options.null = 0;
if isfield(options, 'avoid_zeros')
avoid_zeros = options.avoid_zeros;
else
avoid_zeros = 1;
end
if isfield(options, 'T')
T = options.T;
else
T = -1;
end
if length(T)>1
% build a cell array of histograms
h = {};
for i=1:length(T)
h{end+1} = compute_histogram(M,T(i));
end
return;
end
M = M(:);
if isempty(M)
h = []; x = [];
return;
end
% quantize
if T>0
[Mw, M] = perform_quantization(M, T);
M = double(M);
end
if isfield(options, 'symmetrize')
symmetrize = options.symmetrize;
if ~symmetrize && ~isempty(find(M<=0))
error('Negative entry do not work with options.symmetrize=0');
end
else
symmetrize = ~isempty(find(M<=0));
end
% bounds of histogram
hp = max(abs(M(:)));
hpmax = 10000;
if hp>hpmax
error(['Number of entries in histogram is > ' num2str(hpmax) '.']);
end
hp = min(double(hp),hpmax); % avoid too big histogramms
% compute histrogram
if symmetrize
h = zeros(2*hp+1,1);
x = -hp:hp;
for p = x
h(p+hp+1) = sum(M==p);
end
else
h = zeros(hp,1);
x = 1:hp;
for p = x
h(p) = sum(M==p);
end
end
% avoid zero value
if avoid_zeros
I = find(h==0);
J = find(h>0);
h(I) = mmin(h(J))/1e12;
end
% normalize to sum to 1
h = h/sum(h);
function y = mmin(x)
y = min(x(:));