Mercurial > octave-antonio

comparison scripts/geometry/rectint.m @ 20116:33e706b6b7be

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Carnë Draug <carandraug@octave.org>
date Wed, 22 Apr 2015 20:54:14 +0100
parents 4197fc428c7d
children 7503499a252b
comparison
equal deleted inserted replaced
20115:7e0e8fb16201 20116:33e706b6b7be
1 ## Copyright (C) 2008-2015 Bill Denney 1 ## Copyright (C) 2015 Carnë Draug
2 ## 2 ##
3 ## This file is part of Octave. 3 ## This file is part of Octave.
4 ## 4 ##
5 ## Octave is free software; you can redistribute it and/or modify it 5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by 6 ## under the terms of the GNU General Public License as published by
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b}) 20 ## @deftypefn {Function File} {@var{area} =} rectint (@var{a}, @var{b})
21 ## Compute area or volume of intersection of rectangles or ND boxes.
21 ## 22 ##
22 ## Compute the area of intersection of rectangles in @var{a} and 23 ## Compute the area of intersection of rectangles in @var{a} and
23 ## rectangles in @var{b}. Rectangles are defined as [x y width height] 24 ## rectangles in @var{b}. N dimensional boxes are supported in which
24 ## where x and y are the minimum values of the two orthogonal 25 ## case the volume, or hypervolume is computed according to the number
25 ## dimensions. 26 ## of dimensions.
26 ## 27 ##
27 ## If @var{a} or @var{b} are matrices, then the output, @var{area}, is a 28 ## 2 dimensional rectangles are defined as @code{[xpos ypos width height]}
28 ## matrix where the i-th row corresponds to the i-th row of a and the j-th 29 ## where xpos and ypos are the position of the bottom left corner.
29 ## column corresponds to the j-th row of b. 30 ## Higher dimensions are supported where the coordinates for the minimum
31 ## value of each dimension follow the length of the box in that dimension,
32 ## e.g., @code{[xpos ypos zpos kpos @dots{} width height depth k_length @dots{}]}.
33 ##
34 ## Each row of @var{a} and @var{b} define a rectangle, and if both define
35 ## multiple rectangles, then the output, @var{area}, is a matrix where
36 ## the i-th row corresponds to the i-th row of a and the j-th column
37 ## corresponds to the j-th row of b.
30 ## 38 ##
31 ## @seealso{polyarea} 39 ## @seealso{polyarea}
32 ## @end deftypefn 40 ## @end deftypefn
33 41
34 ## Author: Bill Denney <bill@denney.ws> 42 ## Author: 2015 Carnë Draug <carandraug@octave.org>
35 43
36 function area = rectint (a, b) 44 function dists = rectint (a, b)
37 45
38 if (nargin != 2) 46 if (nargin != 2)
39 print_usage (); 47 print_usage ();
40 elseif (ndims (a) != 2 || ndims (b) != 2) 48 elseif (columns (a) != columns (b))
41 error ("rectint: expecting arguments to be 2-d arrays"); 49 error ("rectint: A and B must have same number of columns");
42 elseif (columns (a) != 4) 50 elseif (ndims (a) > 2)
43 error ("rectint: A must have 4 columns"); 51 error ("rectint: A and B must be 2-d arrays");
44 elseif (columns (b) != 4) 52 elseif (mod (columns (a), 2))
45 error ("rectint: B must have 4 columns"); 53 error ("rectint: number of columns of A and B must be a multiple of two");
46 elseif (any ([a(:,3:4);b(:,3:4)](:) < 0))
47 error ("rectint: all widths and heights must be > 0");
48 endif 54 endif
49 55
50 ## This runs faster if the number of rows of a is greater than the 56 nd = columns (a) / 2;
51 ## number of rows of b. Swap them and transpose to make it run 57 na = rows (a);
52 ## faster. 58 nb = rows (b);
53 swapinputs = false ();
54 if (rows (a) > rows (b))
55 tmp = a;
56 a = b;
57 b = tmp;
58 swapinputs = true ();
59 endif
60 59
61 area = zeros (rows (a), rows (b)); 60 a_start = a(:,1:nd);
62 r1 = [a(:,1:2) a(:,1:2)+a(:,3:4)]; 61 b_start = b(:,1:nd);
63 r2 = [b(:,1:2) b(:,1:2)+b(:,3:4)];
64 for i = 1:columns (area)
65 ## Find the location of each point relative to the other points.
66 r1x1small = r1(:,1) < r2(i,1);
67 r1x1large = r1(:,1) > r2(i,3);
68 r1x1mid = (r1(:,1) >= r2(i,1)) & (r1(:,1) <= r2(i,3));
69 r1x2small = r1(:,3) < r2(i,1);
70 r1x2large = r1(:,3) > r2(i,3);
71 r1x2mid = (r1(:,3) >= r2(i,1)) & (r1(:,3) <= r2(i,3));
72 62
73 r1y1small = r1(:,2) < r2(i,2); 63 a_end = a_start + a(:,nd+1:end);
74 r1y1large = r1(:,2) > r2(i,4); 64 b_end = b_start + b(:,nd+1:end);
75 r1y1mid = (r1(:,2) >= r2(i,2)) & (r1(:,2) <= r2(i,4));
76 r1y2small = r1(:,4) < r2(i,2);
77 r1y2large = r1(:,4) > r2(i,4);
78 r1y2mid = (r1(:,4) >= r2(i,2)) & (r1(:,4) <= r2(i,4));
79 65
80 ## determine the width of the rectangle 66 a_start = reshape (a_start, [na 1 nd]);
81 ## r1 completely encloses r2 67 b_start = reshape (b_start, [1 nb nd]);
82 area(r1x1small & r1x2large,i) = r2(i,3) - r2(i,1);
83 ## the range goes from r2x min to r1x max
84 mask = r1x1small & r1x2mid;
85 area(mask,i) = r1(mask,3) - r2(i,1);
86 ## the range goes from r1x min to r2x max
87 mask = r1x1mid & r1x2large;
88 area(mask,i) = r2(i,3) - r1(mask,1);
89 ## the range goes from r1x min to r1x max
90 mask = r1x1mid & r1x2mid;
91 area(mask,i) = r1(mask,3) - r1(mask,1);
92 68
93 ## determine the height of the rectangle 69 a_end = reshape (a_end, [na 1 nd]);
94 ## r1 completely encloses r2 70 b_end = reshape (b_end, [1 nb nd]);
95 area(r1y1small & r1y2large,i) .*= r2(i,4) - r2(i,2);
96 ## the range goes from r2y min to r1y max
97 mask = r1y1small & r1y2mid;
98 area(mask,i) .*= r1(mask,4) - r2(i,2);
99 ## the range goes from r1y min to r2y max
100 mask = r1y1mid & r1y2large;
101 area(mask,i) .*= r2(i,4) - r1(mask,2);
102 ## the range goes from r1x min to r1x max
103 mask = r1y1mid & r1y2mid;
104 area(mask,i) .*= r1(mask,4) - r1(mask,2);
105 71
106 endfor 72 ## We get a 3d matrix where each dimension is in the 3rd dimension
107 73 dists = bsxfun (@min , a_end, b_end) - bsxfun (@max, a_start, b_start);
108 if (swapinputs) 74 dists(dists < 0) = 0;
109 area = area'; 75 dists = prod (dists, 3);
110 endif
111 76
112 endfunction 77 endfunction
113
114 78
115 ## Exactly overlapping 79 ## Exactly overlapping
116 %!assert (rectint ([0 0 1 1], [0 0 1 1]), 1) 80 %!assert (rectint ([0 0 1 1], [0 0 1 1]), 1)
117 ## rect2 completely enclosed by rect1 81 ## rect2 completely enclosed by rect1
118 %!assert (rectint ([-1 -1 3 3], [0 0 1 1]), 1) 82 %!assert (rectint ([-1 -1 3 3], [0 0 1 1]), 1)
128 %!assert (rectint ([0 0 1 1], [0 1 2 2]), 0) 92 %!assert (rectint ([0 0 1 1], [0 1 2 2]), 0)
129 ## Correct orientation of output 93 ## Correct orientation of output
130 %!assert (rectint ([0 0 1 1;0.5 0.5 1 1;-1 -1 2 2], [1 1 2 2]), [0;0.25;0]) 94 %!assert (rectint ([0 0 1 1;0.5 0.5 1 1;-1 -1 2 2], [1 1 2 2]), [0;0.25;0])
131 %!assert (rectint ([1 1 2 2], [0 0 1 1;0.5 0.5 1 1;-1 -1 2 2]), [0 0.25 0]) 95 %!assert (rectint ([1 1 2 2], [0 0 1 1;0.5 0.5 1 1;-1 -1 2 2]), [0 0.25 0])
132 96
97 ## bug #44904
98 %!assert (rectint ([0 0 5 5], [6 6 5 5]), 0)
99 %!assert (rectint ([0 0 5 5], [0 6 5 5]), 0)
100 %!assert (rectint ([0 0 5 5], [6 0 5 5]), 0)
101 %!assert (rectint ([0 0 0 5 5 5], [0 0 6 5 5 5]), 0)
102
103 ## Test volumes
104 %!shared r1, r2, r3, r4, r5
105 %! r1 = [ 5 3 0 7 5 2];
106 %! r2 = [ 2 5 0 4 2 2];
107 %! r3 = [ 10 7 0 10 3 2];
108 %! r4 = [ 10 -5 0 5 7 2];
109 %! r5 = [-10 0 0 40 11 2];
110
111 %!assert (rectint (r5, r1), 70)
112 %!assert (rectint (r5, r4), 20)
113 %!assert (rectint (r5, [r1; r2; r3; r4]), [70 16 60 20])
114
115 ## Test multiple volumes in both A and B
116 %!assert (rectint ([r2; r5], [r1; r3; r4]), [4 0 0; 70 60 20])