# HG changeset patch # User Carnë Draug # Date 1429732454 -3600 # Node ID 33e706b6b7bed76544be750fd79b7dbc77eebb6f # Parent 7e0e8fb16201bd7a5888b4bd60add61d864a6a71 rectint: fix for non-overlapping rectangle, and support ND boxes (bug #44904) * scripts/geometry/rectint.m: it was giving incorrect results (not zero) when the coordinates for X overlapped but the Y did not (added tests for this bug). This is a complete rewrite of the function. The new version performs approximately 1000 times faster and adds supports N dimensional boxes in which case it measures the volume or hypervolume of the hypervolume intersections. * NEWS: make note that rectint() now supports ND. diff -r 7e0e8fb16201 -r 33e706b6b7be NEWS --- a/NEWS Wed Apr 22 08:41:50 2015 -0700 +++ b/NEWS Wed Apr 22 20:54:14 2015 +0100 @@ -179,7 +179,7 @@ ** The following functions now support N-dimensional arrays: - fliplr flipud rot90 + fliplr flipud rot90 rectint ** The new warning ID "Octave:data-file-in-path" replaces the three previous separate warning IDs "Octave:fopen-file-in-path", diff -r 7e0e8fb16201 -r 33e706b6b7be scripts/geometry/rectint.m --- a/scripts/geometry/rectint.m Wed Apr 22 08:41:50 2015 -0700 +++ b/scripts/geometry/rectint.m Wed Apr 22 20:54:14 2015 +0100 @@ -1,4 +1,4 @@ -## Copyright (C) 2008-2015 Bill Denney +## Copyright (C) 2015 Carnë Draug ## ## This file is part of Octave. ## @@ -18,100 +18,64 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b}) +## Compute area or volume of intersection of rectangles or ND boxes. ## ## Compute the area of intersection of rectangles in @var{a} and -## rectangles in @var{b}. Rectangles are defined as [x y width height] -## where x and y are the minimum values of the two orthogonal -## dimensions. +## rectangles in @var{b}. N dimensional boxes are supported in which +## case the volume, or hypervolume is computed according to the number +## of dimensions. ## -## If @var{a} or @var{b} are matrices, then the output, @var{area}, is a -## matrix where the i-th row corresponds to the i-th row of a and the j-th -## column corresponds to the j-th row of b. +## 2 dimensional rectangles are defined as @code{[xpos ypos width height]} +## where xpos and ypos are the position of the bottom left corner. +## Higher dimensions are supported where the coordinates for the minimum +## value of each dimension follow the length of the box in that dimension, +## e.g., @code{[xpos ypos zpos kpos @dots{} width height depth k_length @dots{}]}. +## +## Each row of @var{a} and @var{b} define a rectangle, and if both define +## multiple rectangles, then the output, @var{area}, is a matrix where +## the i-th row corresponds to the i-th row of a and the j-th column +## corresponds to the j-th row of b. ## ## @seealso{polyarea} ## @end deftypefn -## Author: Bill Denney +## Author: 2015 Carnë Draug -function area = rectint (a, b) +function dists = rectint (a, b) if (nargin != 2) print_usage (); - elseif (ndims (a) != 2 || ndims (b) != 2) - error ("rectint: expecting arguments to be 2-d arrays"); - elseif (columns (a) != 4) - error ("rectint: A must have 4 columns"); - elseif (columns (b) != 4) - error ("rectint: B must have 4 columns"); - elseif (any ([a(:,3:4);b(:,3:4)](:) < 0)) - error ("rectint: all widths and heights must be > 0"); - endif - - ## This runs faster if the number of rows of a is greater than the - ## number of rows of b. Swap them and transpose to make it run - ## faster. - swapinputs = false (); - if (rows (a) > rows (b)) - tmp = a; - a = b; - b = tmp; - swapinputs = true (); + elseif (columns (a) != columns (b)) + error ("rectint: A and B must have same number of columns"); + elseif (ndims (a) > 2) + error ("rectint: A and B must be 2-d arrays"); + elseif (mod (columns (a), 2)) + error ("rectint: number of columns of A and B must be a multiple of two"); endif - area = zeros (rows (a), rows (b)); - r1 = [a(:,1:2) a(:,1:2)+a(:,3:4)]; - r2 = [b(:,1:2) b(:,1:2)+b(:,3:4)]; - for i = 1:columns (area) - ## Find the location of each point relative to the other points. - r1x1small = r1(:,1) < r2(i,1); - r1x1large = r1(:,1) > r2(i,3); - r1x1mid = (r1(:,1) >= r2(i,1)) & (r1(:,1) <= r2(i,3)); - r1x2small = r1(:,3) < r2(i,1); - r1x2large = r1(:,3) > r2(i,3); - r1x2mid = (r1(:,3) >= r2(i,1)) & (r1(:,3) <= r2(i,3)); + nd = columns (a) / 2; + na = rows (a); + nb = rows (b); - r1y1small = r1(:,2) < r2(i,2); - r1y1large = r1(:,2) > r2(i,4); - r1y1mid = (r1(:,2) >= r2(i,2)) & (r1(:,2) <= r2(i,4)); - r1y2small = r1(:,4) < r2(i,2); - r1y2large = r1(:,4) > r2(i,4); - r1y2mid = (r1(:,4) >= r2(i,2)) & (r1(:,4) <= r2(i,4)); + a_start = a(:,1:nd); + b_start = b(:,1:nd); + + a_end = a_start + a(:,nd+1:end); + b_end = b_start + b(:,nd+1:end); - ## determine the width of the rectangle - ## r1 completely encloses r2 - area(r1x1small & r1x2large,i) = r2(i,3) - r2(i,1); - ## the range goes from r2x min to r1x max - mask = r1x1small & r1x2mid; - area(mask,i) = r1(mask,3) - r2(i,1); - ## the range goes from r1x min to r2x max - mask = r1x1mid & r1x2large; - area(mask,i) = r2(i,3) - r1(mask,1); - ## the range goes from r1x min to r1x max - mask = r1x1mid & r1x2mid; - area(mask,i) = r1(mask,3) - r1(mask,1); + a_start = reshape (a_start, [na 1 nd]); + b_start = reshape (b_start, [1 nb nd]); - ## determine the height of the rectangle - ## r1 completely encloses r2 - area(r1y1small & r1y2large,i) .*= r2(i,4) - r2(i,2); - ## the range goes from r2y min to r1y max - mask = r1y1small & r1y2mid; - area(mask,i) .*= r1(mask,4) - r2(i,2); - ## the range goes from r1y min to r2y max - mask = r1y1mid & r1y2large; - area(mask,i) .*= r2(i,4) - r1(mask,2); - ## the range goes from r1x min to r1x max - mask = r1y1mid & r1y2mid; - area(mask,i) .*= r1(mask,4) - r1(mask,2); + a_end = reshape (a_end, [na 1 nd]); + b_end = reshape (b_end, [1 nb nd]); - endfor - - if (swapinputs) - area = area'; - endif + ## We get a 3d matrix where each dimension is in the 3rd dimension + dists = bsxfun (@min , a_end, b_end) - bsxfun (@max, a_start, b_start); + dists(dists < 0) = 0; + dists = prod (dists, 3); endfunction - ## Exactly overlapping %!assert (rectint ([0 0 1 1], [0 0 1 1]), 1) ## rect2 completely enclosed by rect1 @@ -130,3 +94,23 @@ %!assert (rectint ([0 0 1 1;0.5 0.5 1 1;-1 -1 2 2], [1 1 2 2]), [0;0.25;0]) %!assert (rectint ([1 1 2 2], [0 0 1 1;0.5 0.5 1 1;-1 -1 2 2]), [0 0.25 0]) +## bug #44904 +%!assert (rectint ([0 0 5 5], [6 6 5 5]), 0) +%!assert (rectint ([0 0 5 5], [0 6 5 5]), 0) +%!assert (rectint ([0 0 5 5], [6 0 5 5]), 0) +%!assert (rectint ([0 0 0 5 5 5], [0 0 6 5 5 5]), 0) + +## Test volumes +%!shared r1, r2, r3, r4, r5 +%! r1 = [ 5 3 0 7 5 2]; +%! r2 = [ 2 5 0 4 2 2]; +%! r3 = [ 10 7 0 10 3 2]; +%! r4 = [ 10 -5 0 5 7 2]; +%! r5 = [-10 0 0 40 11 2]; + +%!assert (rectint (r5, r1), 70) +%!assert (rectint (r5, r4), 20) +%!assert (rectint (r5, [r1; r2; r3; r4]), [70 16 60 20]) + +## Test multiple volumes in both A and B +%!assert (rectint ([r2; r5], [r1; r3; r4]), [4 0 0; 70 60 20])