Mercurial > hgweb > octave
changeset 32648:475939906cbb stable
rose.m: Overhaul function to get correct results for corner cases (bug #48889)
* rose.m: Add Programming Note to describe algorithm approach. Replace atan2()
with mod() for a more accurate way to map input to range [0, 2*pi). Use ranges
in degrees (integers), rather than radians (floating point), for improved
accuracy. Remove warning about uncentered bins. Rename input to "theta" and
outputs to "th" and "r" for clarity. Tighten test used to warn about bad
plotting results. Add specific BIST tests various corner cases.
| author | Rik <rik@octave.org> |
|---|---|
| date | Tue, 26 Dec 2023 10:22:53 -0800 |
| parents | f75f2abe189c |
| children | 6ec0b7516c5b fb59e77bf356 |
| files | scripts/plot/draw/rose.m |
| diffstat | 1 files changed, 82 insertions(+), 45 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/plot/draw/rose.m +++ b/scripts/plot/draw/rose.m @@ -24,12 +24,12 @@ ######################################################################## ## -*- texinfo -*- -## @deftypefn {} {} rose (@var{th}) -## @deftypefnx {} {} rose (@var{th}, @var{nbins}) -## @deftypefnx {} {} rose (@var{th}, @var{bins}) +## @deftypefn {} {} rose (@var{theta}) +## @deftypefnx {} {} rose (@var{theta}, @var{nbins}) +## @deftypefnx {} {} rose (@var{theta}, @var{bins}) ## @deftypefnx {} {} rose (@var{hax}, @dots{}) ## @deftypefnx {} {@var{h} =} rose (@dots{}) -## @deftypefnx {} {[@var{thout} @var{rout}] =} rose (@dots{}) +## @deftypefnx {} {[@var{th} @var{r}] =} rose (@dots{}) ## Plot an angular histogram. ## ## With one vector argument, @var{th}, plot the histogram with 20 angular bins. @@ -59,16 +59,19 @@ ## @end group ## @end example ## -## Programming Note: When specifying bin centers with the @var{bins} input, -## the edges for bins 2 to N-1 are spaced so that @code{@var{bins}(i)} is -## centered between the edges. The final edge is drawn halfway between bin N -## and bin 1. This guarantees that all input @var{th} will be placed into one -## of the bins, but also means that for some combinations bin 1 and bin N may -## not be centered on the user's given values. ## @seealso{hist, polar} ## @end deftypefn -function [thout, rout] = rose (varargin) +## Programming note: Ranges are calculated in degrees and then converted to +## radians because the use of integers prevents accumulation of small errors +## that result when using floating point directly. +## The histogram counts are calculated using histc(). See the documentation. +## The final count from histc() contains any values *exactly* equal to the +## last bin edge which is always 2*pi. Because the input mapping of +## "mod (th, 2*pi)" changes any 2*pi values to 0, this last bin should always +## be zero and can be safely deleted. + +function [th, r] = rose (varargin) [hax, varargin, nargin] = __plt_get_axis_arg__ ("rose", varargin{:}); @@ -76,49 +79,62 @@ print_usage (); endif - th = varargin{1}; - ## Force theta to [0,2*pi] range - th = atan2 (sin (th), cos (th)); - th(th < 0) += 2*pi; + ## Force theta to range [0,2*pi) + th = mod (varargin{1}, 2*pi); custom_bins = false; if (nargin == 1) - bins = [1/40 : 1/20 : 1] * 2*pi; + bins = [9 : 18 : 360] / 180 * pi; else bins = varargin{2}; if (isscalar (bins)) - bins = [0.5/bins : 1/bins : 1] * 2*pi; + bins = [180/bins : 360/bins : 360] / 180 * pi; else custom_bins = true; - ## Force angles to [0,2*pi] range - bins = atan2 (sin (bins), cos (bins)); - bins(bins < 0) += 2*pi; - bins = unique (bins); + ## Force custom bins to [0,2*pi) range + bins = mod (bins, 2*pi); + bins = unique (bins); # de-duplicate and sort bins + bins = bins(:).'; # Force row vector endif endif - if (numel (bins) < 3) - warning ("rose: bin sizes >= pi will not plot correctly"); + + binedge = bins(1:end-1) + diff (bins) / 2; # halfway between bin centers + if (! custom_bins) + counts = histc (th, [0, binedge, 2*pi]); # Add implicit edges at 0, 2*pi + if (isrow (counts)) + counts = counts(:); + endif + ## FIXME: Remove in Octave 11 if no bug reports filed + if (any (counts(end,:))) + error ("rose: internal error, histc returned count for theta == 2*pi, please file a bug report"); + endif + counts(end,:) = []; # remove temporary bin + else + last_binedge = bins(end) + diff ([bins(end), 2*pi+bins(1)]) / 2; + if (last_binedge >= 2*pi) + counts = histc (th, [0, last_binedge - 2*pi, binedge, 2*pi]); + else + counts = histc (th, [0, binedge, last_binedge, 2*pi]); + endif + if (isrow (counts)) + counts = counts(:); + endif + counts(end-1,:) += counts(1,:); # Combine counts for first, last bin + ## FIXME: Remove in Octave 11 if no bug reports filed + if (any (counts(end,:))) + error ("rose: internal error, histc returned count for theta == 2*pi, please file a bug report"); + endif + counts([1,end], :) = []; # remove temporary bins endif - [counts, binctr] = hist (th, bins); - binctr = binctr(:).'; # Force row vector - if (isvector (counts)) - counts = counts(:); - endif - - binedge = binctr(1:end-1) + diff (binctr) / 2; binedge = [binedge ; zeros(size(binedge)); zeros(size(binedge)); binedge]; binedge = binedge(:); if (! custom_bins) - ## Add in implicit edges at 0 and 2*pi + ## Add implicit edges at 0 and 2*pi th = [0; 0; binedge; 2*pi ; 0]; else - ## Add in final edge - last_bin_edge = binctr(end) + diff ([binctr(end), (2*pi+binctr(1))])/2; - if ((binedge(end) + last_bin_edge)/2 != binctr(end)) - warning ("rose: bin 1 and bin %d are not centered", numel (binctr)); - endif - th = [0; last_bin_edge; binedge; last_bin_edge; 0]; + ## Add final edge for custom bin + th = [0; last_binedge; binedge; last_binedge; 0]; endif r = zeros (4 * rows (counts), columns (counts)); @@ -126,6 +142,10 @@ r(3:4:end, :) = counts; if (nargout < 2) + if (any (diff (bins) >= pi)) + warning ("rose: bin sizes >= pi will not plot correctly"); + endif + oldfig = []; if (! isempty (hax)) oldfig = get (0, "currentfigure"); @@ -140,11 +160,8 @@ end_unwind_protect if (nargout > 0) - thout = htmp; + th = htmp; endif - else - thout = th; - rout = r; endif endfunction @@ -165,10 +182,30 @@ %! rose ([0, 2, 3, 5], [0, pi/2, pi, 3*pi/2]); %! title ("rose() angular histogram plot with specified bins"); +## Test mapping inputs to [0, 2*pi), 2*pi mapped to bin 1. +%!test +%! [t, r] = rose ([1:1:360]/180*pi + 2*pi); +%! assert (diff (t(2:4:end)), 2*pi/20 * ones (19, 1)); +%! assert (r(2:4:end), 18*ones (20, 1)); + +## Custom # of bins, values exactly at 0 and 2*pi go to bin 1 +%!test +%! [t,r] = rose ([0, 2*pi], 4); +%! assert (size (t), [16, 1]); +%! assert (size (r), [16, 1]); +%! assert ([t(2); t(3:4:end)], [0; pi/2; pi; 3*pi/2; 2*pi]); +%! assert (r(2:4:end), [2; 0; 0; 0]); + +## Custom bins, synthesized bin1 cut-off is exactly 36 degrees +%!test +%! [t,r] = rose (deg2rad ([35, 36]), pi * [1/2, 1, 1.5, 1.9]); +%! assert (r(2:4:end), [1; 0; 0; 1]); + +## Custom bins, synthesized bin1 cut-off is exactly -36 degrees +%!test +%! [t,r] = rose (deg2rad ([-36, -37, 360]), pi * [1/10, 1/2, 1, 1.5]); +%! assert (r(2:4:end), [0; 0; 1; 2]); + ## Test input validation %!error <Invalid call> rose () %!error <Invalid call> rose (1,2,3) -%!warning <bin sizes .= pi will not plot correctly> -%! [th, r] = rose ([1 2 2 4 4 4], 2); -%!warning <bin 1 and bin 3 are not centered> -%! [th, r] = rose ([1 2 2 4 4 4], [1 2 3]);
