1## Copyright (C) 2015 Carnë Draug <carandraug@octave.org> 2## 3## This program is free software; you can redistribute it and/or 4## modify it under the terms of the GNU General Public License as 5## published by the Free Software Foundation; either version 3 of the 6## License, or (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, but 9## WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 11## General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program; if not, see 15## <http:##www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {} hist3 (@var{X}) 19## @deftypefnx {Function File} {} hist3 (@var{X}, @var{nbins}) 20## @deftypefnx {Function File} {} hist3 (@var{X}, @qcode{"Nbins"}, @var{nbins}) 21## @deftypefnx {Function File} {} hist3 (@var{X}, @var{centers}) 22## @deftypefnx {Function File} {} hist3 (@var{X}, @qcode{"Ctrs"}, @var{centers}) 23## @deftypefnx {Function File} {} hist3 (@var{X}, @qcode{"Edges"}, @var{edges}) 24## @deftypefnx {Function File} {[@var{N}, @var{C}] =} hist3 (@dots{}) 25## @deftypefnx {Function File} {} hist3 (@dots{}, @var{prop}, @var{val}, @dots{}) 26## @deftypefnx {Function File} {} hist3 (@var{hax}, @dots{}) 27## Produce bivariate (2D) histogram counts or plots. 28## 29## The elements to produce the histogram are taken from the Nx2 matrix 30## @var{X}. Any row with NaN values are ignored. The actual bins can 31## be configured in 3 different: number, centers, or edges of the bins: 32## 33## @table @asis 34## @item Number of bins (default) 35## Produces equally spaced bins between the minimum and maximum values 36## of @var{X}. Defined as a 2 element vector, @var{nbins}, one for each 37## dimension. Defaults to @code{[10 10]}. 38## 39## @item Center of bins 40## Defined as a cell array of 2 monotonically increasing vectors, 41## @var{centers}. The width of each bin is determined from the adjacent 42## values in the vector with the initial and final bin, extending to Infinity. 43## 44## @item Edge of bins 45## Defined as a cell array of 2 monotonically increasing vectors, 46## @var{edges}. @code{@var{N}(i,j)} contains the number of elements 47## in @var{X} for which: 48## 49## @itemize @w{} 50## @item 51## @var{edges}@{1@}(i) <= @var{X}(:,1) < @var{edges}@{1@}(i+1) 52## @item 53## @var{edges}@{2@}(j) <= @var{X}(:,2) < @var{edges}@{2@}(j+1) 54## @end itemize 55## 56## The consequence of this definition is that values outside the initial 57## and final edge values are ignored, and that the final bin only contains 58## the number of elements exactly equal to the final edge. 59## 60## @end table 61## 62## The return values, @var{N} and @var{C}, are the bin counts and centers 63## respectively. These are specially useful to produce intensity maps: 64## 65## @example 66## [counts, centers] = hist3 (data); 67## imagesc (centers@{1@}, centers@{2@}, counts) 68## @end example 69## 70## If there is no output argument, or if the axes graphics handle 71## @var{hax} is defined, the function will plot a 3 dimensional bar 72## graph. Any extra property/value pairs are passed directly to the 73## underlying surface object. 74## 75## @seealso{hist, histc, lookup, mesh} 76## @end deftypefn 77 78function [N, C] = hist3 (X, varargin) 79 if (nargin < 1) 80 print_usage (); 81 endif 82 83 next_argin = 1; 84 should_draw = true; 85 if (isaxes (X)) 86 hax = X; 87 X = varargin{next_argin++}; 88 elseif (nargout == 0) 89 hax = gca (); 90 else 91 should_draw = false; 92 endif 93 94 if (! ismatrix (X) || columns (X) != 2) 95 error ("hist3: X must be a 2 columns matrix"); 96 endif 97 98 method = "nbins"; 99 val = [10 10]; 100 if (numel (varargin) >= next_argin) 101 this_arg = varargin{next_argin++}; 102 if (isnumeric (this_arg)) 103 method = "nbins"; 104 val = this_arg; 105 elseif (iscell (this_arg)) 106 method = "ctrs"; 107 val = this_arg; 108 elseif (numel (varargin) >= next_argin 109 && any (strcmpi ({"nbins", "ctrs", "edges"}, this_arg))) 110 method = tolower (this_arg); 111 val = varargin{next_argin++}; 112 else 113 next_argin--; 114 endif 115 endif 116 117 have_centers = false; 118 switch (tolower (method)) 119 case "nbins" 120 [r_edges, c_edges] = edges_from_nbins (X, val); 121 case "ctrs" 122 have_centers = true; 123 centers = val; 124 [r_edges, c_edges] = edges_from_centers (val); 125 case "centers" 126 ## This was supported until 1.2.4 when the Matlab compatible option 127 ## 'Ctrs' was added. 128 persistent warned = false; 129 if (! warned) 130 warning ("hist3: option `centers' is deprecated. Use `ctrs'"); 131 endif 132 have_centers = true; 133 centers = val; 134 [r_edges, c_edges] = edges_from_centers (val); 135 case "edges" 136 if (! iscell (val) || numel (val) != 2 137 || ! all (cellfun (@isvector, val))) 138 error ("hist3: EDGES must be a cell array with 2 vectors"); 139 endif 140 [r_edges] = vec (val{1}, 2); 141 [c_edges] = vec (val{2}, 2); 142 out_rows = any (X < [r_edges(1) c_edges(1)] 143 | X > [r_edges(end) c_edges(end)], 2); 144 X(out_rows,:) = []; 145 otherwise 146 ## we should never get here... 147 error ("hist3: invalid binning method `%s'", method); 148 endswitch 149 150 ## We only remove the NaN now, after having computed the bin edges, 151 ## because the extremes from each column that define the edges may 152 ## be paired with a NaN. While such values do not appear on the 153 ## histogram, they must still be used to compute the histogram 154 ## edges. 155 X(any (isnan (X), 2), :) = []; 156 157 r_idx = lookup (r_edges, X(:,1), "l"); 158 c_idx = lookup (c_edges, X(:,2), "l"); 159 160 counts_size = [numel(r_edges) numel(c_edges)]; 161 counts = accumarray ([r_idx, c_idx], 1, counts_size); 162 163 if (should_draw) 164 counts = counts.'; 165 z = zeros ((size (counts) +1) *2); 166 z(2:end-1,2:end-1) = kron (counts, ones (2, 2)); 167 ## Setting the values for the end of the histogram bin like this 168 ## seems straight wrong but that's hwo Matlab plots look. 169 y = [kron(c_edges, ones (1, 2)) (c_edges(end)*2-c_edges(end-1))([1 1])]; 170 x = [kron(r_edges, ones (1, 2)) (r_edges(end)*2-r_edges(end-1))([1 1])]; 171 mesh (hax, x, y, z, "facecolor", [.75 .85 .95], varargin{next_argin:end}); 172 else 173 N = counts; 174 if (isargout (2)) 175 if (! have_centers) 176 C = {(r_edges + [diff(r_edges)([1:end end])]/ 2) ... 177 (c_edges + [diff(c_edges)([1:end end])]/ 2)}; 178 else 179 C = centers(:)'; 180 C{1} = vec (C{1}, 2); 181 C{2} = vec (C{2}, 2); 182 endif 183 endif 184 endif 185 186endfunction 187 188function [r_edges, c_edges] = edges_from_nbins (X, nbins) 189 if (! isnumeric (nbins) || numel (nbins) != 2) 190 error ("hist3: NBINS must be a 2 element vector"); 191 endif 192 inits = min (X, [], 1); 193 ends = max (X, [], 1); 194 ends -= (ends - inits) ./ vec (nbins, 2); 195 196 ## If any histogram side has an empty range, then still make NBINS 197 ## but then place that value at the centre of the centre bin so that 198 ## they appear in the centre in the plot. 199 single_bins = inits == ends; 200 if (any (single_bins)) 201 inits(single_bins) -= (floor (nbins(single_bins) ./2)) + 0.5; 202 ends(single_bins) = inits(single_bins) + nbins(single_bins) -1; 203 endif 204 205 r_edges = linspace (inits(1), ends(1), nbins(1)); 206 c_edges = linspace (inits(2), ends(2), nbins(2)); 207endfunction 208 209function [r_edges, c_edges] = edges_from_centers (ctrs) 210 if (! iscell (ctrs) || numel (ctrs) != 2 || ! all (cellfun (@isvector, ctrs))) 211 error ("hist3: CTRS must be a cell array with 2 vectors"); 212 endif 213 r_edges = vec (ctrs{1}, 2); 214 c_edges = vec (ctrs{2}, 2); 215 r_edges(2:end) -= diff (r_edges) / 2; 216 c_edges(2:end) -= diff (c_edges) / 2; 217endfunction 218 219%!demo 220%! X = [ 221%! 1 1 222%! 1 1 223%! 1 10 224%! 1 10 225%! 5 5 226%! 5 5 227%! 5 5 228%! 5 5 229%! 5 5 230%! 7 3 231%! 7 3 232%! 7 3 233%! 10 10 234%! 10 10]; 235%! hist3 (X) 236 237%!test 238%! N_exp = [ 0 0 0 5 20 239%! 0 0 10 15 0 240%! 0 15 10 0 0 241%! 20 5 0 0 0]; 242%! 243%! n = 100; 244%! x = [1:n]'; 245%! y = [n:-1:1]'; 246%! D = [x y]; 247%! N = hist3 (D, [4 5]); 248%! assert (N, N_exp); 249 250%!test 251%! N_exp = [0 0 0 0 1 252%! 0 0 0 0 1 253%! 0 0 0 0 1 254%! 1 1 1 1 93]; 255%! 256%! n = 100; 257%! x = [1:n]'; 258%! y = [n:-1:1]'; 259%! D = [x y]; 260%! C{1} = [1 1.7 3 4]; 261%! C{2} = [1:5]; 262%! N = hist3 (D, C); 263%! assert (N, N_exp); 264 265## bug 44987 266%!test 267%! D = [1 1; 3 1; 3 3; 3 1]; 268%! [c, nn] = hist3 (D, {0:4, 0:4}); 269%! exp_c = zeros (5); 270%! exp_c([7 9 19]) = [1 2 1]; 271%! assert (c, exp_c); 272%! assert (nn, {0:4, 0:4}); 273 274%!test 275%! for i = 10 276%! assert (size (hist3 (rand (9, 2), "Edges", {[0:.2:1]; [0:.2:1]})), [6 6]) 277%! endfor 278 279%!test 280%! edge_1 = linspace (0, 10, 10); 281%! edge_2 = linspace (0, 50, 10); 282%! [c, nn] = hist3 ([1:10; 1:5:50]', "Edges", {edge_1, edge_2}); 283%! exp_c = zeros (10, 10); 284%! exp_c([1 12 13 24 35 46 57 68 79 90]) = 1; 285%! assert (c, exp_c); 286%! 287%! assert (nn{1}, edge_1 + edge_1(2)/2, eps*10^4) 288%! assert (nn{2}, edge_2 + edge_2(2)/2, eps*10^4) 289 290%!shared X 291%! X = [ 292%! 5 2 293%! 5 3 294%! 1 4 295%! 5 3 296%! 4 4 297%! 1 2 298%! 2 3 299%! 3 3 300%! 5 4 301%! 5 3]; 302 303%!test 304%! N = zeros (10); 305%! N([1 10 53 56 60 91 98 100]) = [1 1 1 1 3 1 1 1]; 306%! C = {(1.2:0.4:4.8), (2.1:0.2:3.9)}; 307%! assert (nthargout ([1 2], @hist3, X), {N C}, eps*10^3) 308 309%!test 310%! N = zeros (5, 7); 311%! N([1 5 17 18 20 31 34 35]) = [1 1 1 1 3 1 1 1]; 312%! C = {(1.4:0.8:4.6), ((2+(1/7)):(2/7):(4-(1/7)))}; 313%! assert (nthargout ([1 2], @hist3, X, [5 7]), {N C}, eps*10^3) 314%! assert (nthargout ([1 2], @hist3, X, "Nbins", [5 7]), {N C}, eps*10^3) 315 316%!test 317%! N = [0 1 0; 0 1 0; 0 0 1; 0 0 0]; 318%! C = {(2:5), (2.5:1:4.5)}; 319%! assert (nthargout ([1 2], @hist3, X, "Edges", {(1.5:4.5), (2:4)}), {N C}) 320 321%!test 322%! N = [0 0 1 0 1 0; 0 0 0 1 0 0; 0 0 1 4 2 0]; 323%! C = {(1.2:3.2), (0:5)}; 324%! assert (nthargout ([1 2], @hist3, X, "Ctrs", C), {N C}) 325%! assert (nthargout ([1 2], @hist3, X, C), {N C}) 326 327%!test 328%! [~, C] = hist3 (rand (10, 2), "Edges", {[0 .05 .15 .35 .55 .95], 329%! [-1 .05 .07 .2 .3 .5 .89 1.2]}); 330%! C_exp = {[ 0.025 0.1 0.25 0.45 0.75 1.15], ... 331%! [-0.475 0.06 0.135 0.25 0.4 0.695 1.045 1.355]}; 332%! assert (C, C_exp, eps*10^2) 333 334## Test how handling of out of borders is different whether we are 335## defining Centers or Edges. 336%!test 337%! Xv = repmat ([1:10]', [1 2]); 338%! 339%! ## Test Centers 340%! assert (hist3 (Xv, "Ctrs", {1:10, 1:10}), eye (10)) 341%! 342%! N_exp = eye (6); 343%! N_exp([1 end]) = 3; 344%! assert (hist3 (Xv, "Ctrs", {3:8, 3:8}), N_exp) 345%! 346%! N_exp = zeros (8, 6); 347%! N_exp([1 2 11 20 29 38 47 48]) = [2 1 1 1 1 1 1 2]; 348%! assert (hist3 (Xv, "Ctrs", {2:9, 3:8}), N_exp) 349%! 350%! ## Test Edges 351%! assert (hist3 (Xv, "Edges", {1:10, 1:10}), eye (10)) 352%! assert (hist3 (Xv, "Edges", {3:8, 3:8}), eye (6)) 353%! assert (hist3 (Xv, "Edges", {2:9, 3:8}), [zeros(1, 6); eye(6); zeros(1, 6)]) 354%! 355%! N_exp = zeros (14); 356%! N_exp(3:12, 3:12) = eye (10); 357%! assert (hist3 (Xv, "Edges", {-1:12, -1:12}), N_exp) 358%! 359%! ## Test for Nbins 360%! assert (hist3 (Xv), eye (10)) 361%! assert (hist3 (Xv, [10 10]), eye (10)) 362%! assert (hist3 (Xv, "nbins", [10 10]), eye (10)) 363%! assert (hist3 (Xv, [5 5]), eye (5) * 2) 364%! 365%! N_exp = zeros (7, 5); 366%! N_exp([1 9 10 18 26 27 35]) = [2 1 1 2 1 1 2]; 367%! assert (hist3 (Xv, [7 5]), N_exp) 368 369%!test # bug #51059 370%! D = [1 1; NaN 2; 3 1; 3 3; 1 NaN; 3 1]; 371%! [c, nn] = hist3 (D, {0:4, 0:4}); 372%! exp_c = zeros (5); 373%! exp_c([7 9 19]) = [1 2 1]; 374%! assert (c, exp_c) 375%! assert (nn, {0:4, 0:4}) 376 377## Single row of data or cases where all elements have the same value 378## on one side of the histogram. 379%!test 380%! [c, nn] = hist3 ([1 8]); 381%! exp_c = zeros (10, 10); 382%! exp_c(6, 6) = 1; 383%! exp_nn = {-4:5, 3:12}; 384%! assert (c, exp_c) 385%! assert (nn, exp_nn, eps) 386%! 387%! [c, nn] = hist3 ([1 8], [10 11]); 388%! exp_c = zeros (10, 11); 389%! exp_c(6, 6) = 1; 390%! exp_nn = {-4:5, 3:13}; 391%! assert (c, exp_c) 392%! assert (nn, exp_nn, eps) 393 394## NaNs paired with values defining the histogram edges. 395%!test 396%! [c, nn] = hist3 ([1 NaN; 2 3; 6 9; 8 NaN]); 397%! exp_c = zeros (10, 10); 398%! exp_c(2, 1) = 1; 399%! exp_c(8, 10) = 1; 400%! exp_nn = {linspace(1.35, 7.65, 10) linspace(3.3, 8.7, 10)}; 401%! assert (c, exp_c) 402%! assert (nn, exp_nn, eps*100) 403 404## Columns full of NaNs (recent Matlab versions seem to throw an error 405## but this did work like this on R2010b at least). 406%!test 407%! [c, nn] = hist3 ([1 NaN; 2 NaN; 6 NaN; 8 NaN]); 408%! exp_c = zeros (10, 10); 409%! exp_nn = {linspace(1.35, 7.65, 10) NaN(1, 10)}; 410%! assert (c, exp_c) 411%! assert (nn, exp_nn, eps*100) 412 413## Behaviour of an empty X after removal of rows with NaN. 414%!test 415%! [c, nn] = hist3 ([1 NaN; NaN 3; NaN 9; 8 NaN]); 416%! exp_c = zeros (10, 10); 417%! exp_nn = {linspace(1.35, 7.65, 10) linspace(3.3, 8.7, 10)}; 418%! assert (c, exp_c) 419%! assert (nn, exp_nn, eps*100) 420