Index: /issm/trunk/externalpackages/canos/InterX.m
===================================================================
--- /issm/trunk/externalpackages/canos/InterX.m	(revision 10039)
+++ /issm/trunk/externalpackages/canos/InterX.m	(revision 10039)
@@ -0,0 +1,82 @@
+function P = InterX(L1,varargin)
+%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.
+%   
+%   Example:
+%       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')
+
+%   Author : NS
+%   Version: 3.0, 21 Sept. 2010
+
+%   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
+       
+    %...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:end-1) - dy1.*x1(1:end-1);
+    S2 = dx2.*y2(1:end-1) - dy2.*x2(1:end-1);
+    
+    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)';
+
+    %...Obtain the segments where an intersection is expected
+    [i,j] = find(C1 & C2); 
+    if isempty(i),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],'rows')';
+              
+    function u = D(x,y)
+        u = bsxfun(@minus,x(:,1:end-1),y).*bsxfun(@minus,x(:,2:end),y);
+    end
+end
Index: /issm/trunk/externalpackages/canos/licenseNS.txt
===================================================================
--- /issm/trunk/externalpackages/canos/licenseNS.txt	(revision 10039)
+++ /issm/trunk/externalpackages/canos/licenseNS.txt	(revision 10039)
@@ -0,0 +1,24 @@
+Copyright (c) 2009, NS
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without 
+modification, are permitted provided that the following conditions are 
+met:
+
+    * Redistributions of source code must retain the above copyright 
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright 
+      notice, this list of conditions and the following disclaimer in 
+      the documentation and/or other materials provided with the distribution
+      
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+POSSIBILITY OF SUCH DAMAGE.