Functions for line intersections

8756e0b0-4786-49cc-b282-8d7a293f7675 65 Functions for line intersections Marco Bernacca 4 5 decimal

&[DATE] &[TIME] - &[FILENAME]

These functions have been written to implement in Smath some features that are natively provided by either Scilab or Matlab in order to make available in Smath a function originally written in Matlab to find intersections between lines. The function in subject is P = InterX(L1, L2), written by NS. See below link: https://de.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections I made a translation for Scilab and (tried to) translate for Smath. To simplify the task, I have written some helper functions that might be useful for other tasks, so I provide them here to the community.

Function mdiff(V) Input: a vector (either a row or a column) Output: a vector of the same nature of input (a row or a column) containing the differences V(i+1)-V(i). Number of output elements are lenght(V)-1. The first line could appear a bit weird, but it is the quickest way I found to mantain the original nature of the input vector. Another possible implementation is given below.

V mdiff R V rows V rows 1 > (- V cols V cols 1 > (- matrix : ii 1 R length range R ii el V ii 1 + el V ii el - : for R 31 line :

Examples:

z2 4273151 mat :

V mdiff2 lR V length 1 - : V rows 1 > i 1 lR range R i el V i 1 + el V i el - : for i 1 lR range R 1 i el V i 1 + el V i el - : for if R 31 line :

z1 13589416 mat :

z1 mdiff 22315 - 15 mat

z2 mdiff 2 - 54 - 2 - 41 mat

Function map(x; y; fun(2)) Input: x, y: vectors or matrices with compatible sizes. fun: a function of 2 argumens Output: vector or matrix whose elements are fun(el(x);el(y)) Apply a given function to corresponding elements of compatible arrays. For info on the Matlab concept of "arrays of compatible sizes" please refer to: https://it.mathworks.com/help/matlab/matlab_prog/compatible-array-sizes-for-basic-operations.html

a b plus a b + :

a b times a b * :

a b minus a b - :

a b divide a b / :

x y 2 fun map xh x cols : xv x rows : yh y cols : yv y rows : okx xh yh ≡ (xh 1 ≡ (yh 1 ≡ (| (| (: oky xv yv ≡ (xv 1 ≡ (yv 1 ≡ (| (| (: okx oky & xx x : yy y : xh yh < ii 2 yh range xx xx x augment : for xh yh > ii 2 xh range yy yy y augment : for 1 if if xv yv < ii 2 yv range xx xx x stack : for xv yv > ii 2 xv range yy yy y stack : for 1 if if ii 1 xv yv Max range jj 1 xh yh Max range zz ii jj el xx ii jj el yy ii jj el fun : for for zz 61 line Incompatible sizes error if 71 line :

z1 z2 A B times map 4122032361626101618872135566328391524271213589456 mat

Examples:

Function index(x; cond) Input: x: a vector (row or column). cond: expression evaluating to an array of T/F values. Output: a vector of same nature of input array whose components are those corresponding to a TRUE value for the given condition.

V cond index R V rows V cols matrix : j 1 : i 1 V length range cond vectorize i el R j el V i el : j j 1 + : 21 line continue if for R rows 1 > (R 1 j 1 - 11 submatrix R 111 j 1 - submatrix if 41 line :

Examples:

z1 z1 4 > index 58913 mat

z2 z1 mdiff 0 > index 427341 mat

Function plotG(x; y; char; size; color) Thanks to Omorr!! :-)

x y char size color plotG n x length : plot x 1 el y 1 el char size color augment : i 2 n range plot plot x i el y i el char size color augment stack : for plot 41 line :

Function D(x; y) Helper function for InterX. See below. Sorry for the long line, but I do like functional programming... :-) At the end this is simpler than an average lisp expression!!

x y D x 1 x rows 1 x cols 1 - submatrix y A B minus map x 1 x rows 2 x cols submatrix y A B minus map A B times map :

function P = InterX(L1,L2) //INTERX Intersection of curves // P = INTERX(L1,L2) returns the intersection points of two curves L1 // and L2. The curves L1,L2 can be either closed or open and are described // by two-row-matrices, where each row contains its x- and y- coordinates. // The intersection of groups of curves (e.g. contour lines, multiply // connected regions etc) can also be computed by separating them with a // column of NaNs as for example // // L = [x11 x12 x13 ... NaN x21 x22 x23 ...; // y11 y12 y13 ... NaN y21 y22 y23 ...] // // P has the same structure as L1 and L2, and its rows correspond to the // x- and y- coordinates of the intersection points of L1 and L2. If no // intersections are found, the returned P is empty. // // P = INTERX(L1) returns the self-intersection points of L1. To keep // the code simple, the points at which the curve is tangent to itself are // not included. P = INTERX(L1,L1) returns all the points of the curve // together with any self-intersection points. // // Author : NS // Version: 3.0, 21 Sept. 2010 // Modified for Scilab by Marco Bernacca 07 March 2017 // Two words about the algorithm: Most of the code is self-explanatory. // The only trick lies in the calculation of C1 and C2. To be brief, this // is essentially the two-dimensional analog of the condition that needs // to be satisfied by a function F(x) that has a zero in the interval // [a,b], namely // F(a)*F(b) <= 0 // C1 and C2 exactly do this for each segment of curves 1 and 2 // respectively. If this condition is satisfied simultaneously for two // segments then we know that they will cross at some point. // Each factor of the 'C' arrays is essentially a matrix containing // the numerators of the signed distances between points of one curve // and line segments of the other. /* //...Argument checks and assignment of L2 error(nargchk(1,2,nargin)); if nargin == 1, L2 = L1; hF = @lt; //...Avoid the inclusion of common points else L2 = varargin{1}; hF = @le; end Note by MBernacca: this is not implemented in the Scilab version. L1 and L2 must be both specified and sintax is P = InterX(L1, L2)) */ //...Preliminary stuff x1 = L1(1,:)'; x2 = L2(1,:); y1 = L1(2,:)'; y2 = L2(2,:); dx1 = diff(x1); dy1 = diff(y1); dx2 = diff(x2); dy2 = diff(y2); //...Determine 'signed distances' S1 = dx1.*y1(1:$-1) - dy1.*x1(1:$-1); S2 = dx2.*y2(1:$-1) - dy2.*x2(1:$-1); // Original Matlab statements // C1 = feval(hF,D(bsxfun(@times,dx1,y2)-bsxfun(@times,dy1,x2),S1),0); // C2 = feval(hF,D((bsxfun(@times,y1,dx2)-bsxfun(@times,x1,dy2))',S2'),0)'; // Modified for Scilab C1 = D(bsrfun(dx1,y2,times)-bsrfun(dy1,x2,times),S1) <= 0; C2 = (D((bsrfun(y1,dx2,times)-bsrfun(x1,dy2,times))',S2') <= 0)'; //...Obtain the segments where an intersection is expected [i,j] = find(C1 & C2); if isempty(i) then P = zeros(2,0);return; end; //...Transpose and prepare for output i=i'; dx2=dx2'; dy2=dy2'; S2 = S2'; L = dy2(j).*dx1(i) - dy1(i).*dx2(j); i = i(L~=0); j=j(L~=0); L=L(L~=0); //...Avoid divisions by 0 //...Solve system of eqs to get the common points P = unique([dx2(j).*S1(i) - dx1(i).*S2(j), dy2(j).*S1(i) - dy1(i).*S2(j)]./[L L],'r')'; endfunction /* function u = D(x,y) u = bsxfun(@minus,x(:,1:end-1),y).*bsxfun(@minus,x(:,2:end),y); end Note by MBer: Matlab function u = bsxfun(fun, a, b) substituted by function u = bsrfun(a, b, fun) written to support "compatible sizes" feature not available in Scilab */ //========================================================== // Below functions written by Marco Bernacca 07 March 2017 //========================================================== function u = D(a, b) u = bsrfun(a(:, 1:$-1), b, minus) .* bsrfun(a(:, 2:$), b, minus); endfunction //_____________________________________________________________________________ // Function: Apply a function to a set of elements // Developer: Marco Bernacca // Origination: 2017-03-02 // Edited: 2017-03-02 //_____________________________________________________________________________ function u = bsrfun(a, b, fun) [chk, ra, rb, ms] = chkarg(a,b); // Replicates Matlab's "Compatible sizes" feature if chk then a = repmat(a, ra); b = repmat(b, rb); u = fun(a, b); else disp ("bsrfun: Incompatible sizes of arguments.") abort end endfunction //========================================================== function [c,ra,rb,d] = chkarg(a,b) // Helper function to support Matlab's "Compatible sizes" feature aa = size(a); bb = size(b); na = ndims(a); nb = ndims(b); if na > nb then bb = [bb, ones(na-nb)]; elseif na < nb then aa = [aa, ones(nb-na)]; end cc = min(aa,bb)(aa<>bb); c = and(cc == [] || cc == 1); d = max(aa,bb); ra = ones(1, length(d)); rb = ones(1, length(d)); ra(aa==1) = d(aa==1); rb(bb==1) = d(bb==1); endfunction //========================================================== function u = plus(a, b) u = a + b; endfunction //========================================================== function u = minus(a, b) u = a - b; endfunction //========================================================== function u = times(a, b) u = a .* b; endfunction //========================================================== function u = divide(a, b) u = a ./ b; endfunction //========================================================== // Scilab users: just copy/paste this entire text region to a Scilab console. // Example: timer(); t = linspace(0,2*%pi); r1 = sin(4*t)+2; x1 = r1.*cos(t); y1 = r1.*sin(t); r2 = sin(8*t)+2; x2 = r2.*cos(t); y2 = r2.*sin(t); P = InterX([x1;y1],[x2;y2]); plot(x1,y1,x2,y2,P(1,:),P(2,:),'ro') timer() //End of File_____________________________________________________________________________

Scilab graphic output

Now let's try to do the same thing in Smath. Some plugins are required and in particular: - Matlab C++ Math Library Attention: to my knowledge,this plugin only works with 32 bit version of Smath. I am using 0.98 build 6484 dated 02 Oct 2017.

Function InterX(L1; L2) Input: L1; L2 : 2-rows matrices of coords (Row1 = x-coord; Row2 = y-coord) Output: 2-rows matricx of coords of intersecting points

L1 L2 InterX x1 L1 1 row transpose : x2 L2 1 row : y1 L1 2 row transpose : y2 L2 2 row : dx1 x1 mdiff : dy1 y1 mdiff : dx2 x2 mdiff : dy2 y2 mdiff : ii 1 dx1 length range S1 ii el dx1 ii el y1 ii el * dy1 ii el x1 ii el * - : S2 ii el dx2 ii el y2 ii el * dy2 ii el x2 ii el * - : 21 line for S2 S2 transpose eval : C1 dx1 y2 A B times map dy1 x2 A B times map - S1 D 0 mwle : C2 y1 dx2 A B times map x1 dy2 A B times map - (transpose S2 transpose D 0 mwle transpose : ff C1 C2 mwand mwfind : ii ff 1 el eval : jj ff 2 el eval : ii length 0 ≡ 21 Zeros ii ii transpose : dx2 dx2 transpose : dy2 dy2 transpose : S2 S2 transpose : L dy2 jj el dx1 ii el A B times map dy1 ii el dx2 jj el A B times map - : ii ii L 0 ≠ index : jj jj L 0 ≠ index : L L L 0 ≠ index : PP1 dx2 jj el S1 ii el A B times map dx1 ii el S2 jj el A B times map - : PP2 dy2 jj el S1 ii el A B times map dy1 ii el S2 jj el A B times map - : PP PP1 PP2 augment L L augment A B divide map : PP 121 line if 161 line :

Let's test it with a simple example:

L1 123456890112203 0.5 28 mat :

L2 345665790133202128 mat :

tt L1 L2 InterX : 4.5 2 5.5 1 5.6667 0.6667 7.25 1.875 8.75 1.125 52 mat

plot1 tt 1 col tt 2 col o 12 Green plotG :

L1 transpose L2 transpose plot1 3 1 sys

and a bit more complicate one:

t 02 π * (2 π 100 / * range :

r1 4 t * sin vectorize 2 + :

r2 8 t * sin vectorize 2 + :

x1 r1 t cos vectorize A B times map :

y1 r1 t sin vectorize A B times map :

x2 r2 t cos vectorize A B times map :

y2 r2 t sin vectorize A B times map :

arg1 x1 transpose y1 transpose stack : 2 2.2443 2.4622 2.637 2.755 2.8066 2.7875 2.6985 2.5455 0 0.1412 0.311 0.503 0.7074 0.9119 1.1036 1.2698 1.3994 29 mat

arg2 x2 transpose y2 transpose stack : 2 2.4769 2.8219 2.9449 2.8136 2.4611 1.9761 1.4766 1.0774 0 0.1558 0.3565 0.5618 0.7224 0.7997 0.7824 0.6948 0.5923 29 mat

aa arg1 arg2 InterX :

plot2 aa 1 col aa 2 col o 12 Green plotG :

arg1 transpose arg2 transpose plot2 3 1 sys

Great! It works! So we have Smath performing like Scilab! Well, almost... The timing of Smath on my hardware says:

Whilst, on same hardware, Scilab says:

and this is expressed in seconds... But anyway Smath is a great tool!!

In the hope that these brief notes could be of help to somebody, enjoy! Marco