1## Copyright (C) 2021 David Legland
2## All rights reserved.
3##
4## Redistribution and use in source and binary forms, with or without
5## modification, are permitted provided that the following conditions are met:
6##
7##     1 Redistributions of source code must retain the above copyright notice,
8##       this list of conditions and the following disclaimer.
9##     2 Redistributions in binary form must reproduce the above copyright
10##       notice, this list of conditions and the following disclaimer in the
11##       documentation and/or other materials provided with the distribution.
12##
13## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
14## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
16## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
17## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
18## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
19## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
20## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
21## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
22## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23##
24## The views and conclusions contained in the software and documentation are
25## those of the authors and should not be interpreted as representing official
26## policies, either expressed or implied, of the copyright holders.
27
28function b = isPointOnEdge(point, edge, varargin)
29%ISPOINTONEDGE Test if a point belongs to an edge.
30%
31%   Usage
32%   B = isPointOnEdge(POINT, EDGE)
33%   B = isPointOnEdge(POINT, EDGE, TOL)
34%
35%   Description
36%   B = isPointOnEdge(POINT, EDGE)
37%   with POINT being [xp yp], and EDGE being [x1 y1 x2 y2], returns TRUE if
38%   the point is located on the edge, and FALSE otherwise.
39%
40%   B = isPointOnEdge(POINT, EDGE, TOL)
41%   Specify an optilonal tolerance value TOL. The tolerance is given as a
42%   fraction of the norm of the edge direction vector. Default is 1e-14.
43%
44%   B = isPointOnEdge(POINTARRAY, EDGE)
45%   B = isPointOnEdge(POINT, EDGEARRAY)
46%   When one of the inputs has several rows, return the result of the test
47%   for each element of the array tested against the single parameter.
48%
49%   B = isPointOnEdge(POINTARRAY, EDGEARRAY)
50%   When both POINTARRAY and EDGEARRAY have the same number of rows,
51%   returns a column vector with the same number of rows.
52%   When the number of rows are different and both greater than 1, returns
53%   a Np-by-Ne matrix of booleans, containing the result for each couple of
54%   point and edge.
55%
56%   Examples
57%   % create a point array
58%   points = [10 10;15 10; 30 10];
59%   % create an edge array
60%   vertices = [10 10;20 10;20 20;10 20];
61%   edges = [vertices vertices([2:end 1], :)];
62%
63%   % Test one point and one edge
64%   isPointOnEdge(points(1,:), edges(1,:))
65%   ans =
66%       1
67%   isPointOnEdge(points(3,:), edges(1,:))
68%   ans =
69%       0
70%
71%   % Test one point and several edges
72%   isPointOnEdge(points(1,:), edges)'
73%   ans =
74%        1     0     0     1
75%
76%   % Test several points and one edge
77%   isPointOnEdge(points, edges(1,:))'
78%   ans =
79%        1     1     0
80%
81%   % Test N points and N edges
82%   isPointOnEdge(points, edges(1:3,:))'
83%   ans =
84%        1     0     0
85%
86%   % Test NP points and NE edges
87%   isPointOnEdge(points, edges)
88%   ans =
89%        1     0     0     1
90%        1     0     0     0
91%        0     0     0     0
92%
93%
94%   See also
95%   edges2d, points2d, isPointOnLine
96%
97
98%   ---------
99%   author : David Legland
100%   INRA - TPV URPOI - BIA IMASTE
101%   created the 31/10/2003.
102%
103
104%   HISTORY
105%   11/03/2004 change input format: edge is [x1 y1 x2 y2].
106%   17/01/2005 if test N edges with N points, return N boolean.
107%   21/01/2005 normalize test for colinearity, so enhance precision
108%   22/05/2009 rename to isPointOnEdge, add psb to specify tolerance
109%   26/01/2010 fix bug in precision computation
110%   04/10/2010 fix a bug, and clean up code
111%   28/10/2010 fix bug to have N results when input is N points and N
112%       edges, add support for arrays with different numbers of rows, and
113%       update doc.
114%   2011-06-15 rewrites by using less memory, and avoiding repmat when psb
115
116
117% extract computation tolerance
118tol = 1e-14;
119if ~isempty(varargin)
120    tol = varargin{1};
121end
122
123% number of edges and of points
124nPoints = size(point, 1);
125nEdges = size(edge, 1);
126
127% adapt size of inputs if needed, and extract elements for computation
128if nPoints == nEdges
129    % When the number of points and edges is the same, the one-to-one test
130    % will be computed, so there is no need to repeat matrices
131    dx = edge(:,3) - edge(:,1);
132    dy = edge(:,4) - edge(:,2);
133    lx = point(:,1) - edge(:,1);
134    ly = point(:,2) - edge(:,2);
135
136elseif nPoints == 1
137    % one point, several edges
138    dx = edge(:, 3) - edge(:, 1);
139    dy = edge(:, 4) - edge(:, 2);
140    lx = point(ones(nEdges, 1), 1) - edge(:, 1);
141    ly = point(ones(nEdges, 1), 2) - edge(:, 2);
142
143elseif nEdges == 1
144    % several points, one edge
145    dx = (edge(3) - edge(1)) * ones(nPoints, 1);
146    dy = (edge(4) - edge(2)) * ones(nPoints, 1);
147    lx = point(:, 1) - edge(1);
148    ly = point(:, 2) - edge(2);
149
150else
151    % Np points and Ne edges:
152    % Create an array for each parameter, so that the result will be a
153    % Np-by-Ne matrix of booleans (requires more memory, and uses repmat)
154
155    x0 = repmat(edge(:, 1)', nPoints, 1);
156    y0 = repmat(edge(:, 2)', nPoints, 1);
157    dx = repmat(edge(:, 3)', nPoints,  1) - x0;
158    dy = repmat(edge(:, 4)', nPoints,  1) - y0;
159
160    lx = repmat(point(:, 1), 1, nEdges) - x0;
161    ly = repmat(point(:, 2), 1, nEdges) - y0;
162end
163
164% test if point is located on supporting line
165b1 = abs(lx.*dy - ly.*dx) ./ (dx.*dx + dy.*dy) < tol;
166
167% compute position of point with respect to edge bounds
168% use different tests depending on line angle
169ind     = abs(dx) > abs(dy);
170t       = zeros(size(dx));
171t(ind)  = lx( ind) ./ dx( ind);
172t(~ind) = ly(~ind) ./ dy(~ind);
173
174% check if point is located between edge bounds
175b = t >- tol & t-1 < tol & b1;
176