[interpolation] mesh2d gateway introduced
[scilab.git] / scilab / modules / interpolation / macros / mesh2d.sci
1 //
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 //
4 // Copyright (C) 2018 - St├ęphane MOTTELET
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // For more information, see the COPYING file which you should have received
8 // along with this program.
9 //
10
11 function [tri, hull]=mesh2d(varargin)
12     if argn(2) == 0 then
13         x = rand(10,1,"normal");
14         y = rand(10,1,"normal");
15         [tri, hull]=mesh2d(x,y);
16         tri = [tri;tri(1,:)];
17         plot(x,y,'o');
18         xstring(x,y,string(1:size(x,"*")));
19         xpolys(matrix(x(tri),size(tri)),matrix(y(tri),size(tri)));
20         plot(x(hull),y(hull),'-r');
21         return
22     end
23     if argn(2) < 2 || argn(2) > 3
24         error(msprintf(_("%s: Wrong number of input argument(s): %d or %d expected.\n"), "mesh2d", 2, 3));
25     end
26     if argn(2) == 3 && argn(1) <> 1
27         error(msprintf(_("%s: Wrong number of output argument(s): %d expected.\n"), "mesh2d", 1));
28     end
29     if argn(2) == 2 && argn(1) > 2
30         error(msprintf(_("%s: Wrong number of output argument(s): %d or %d expected.\n"), "mesh2d", 1, 2));
31     end
32     for i = 1:length(varargin)
33         if type(varargin(i)) <> 1 || ~isreal(varargin(i))
34             error(msprintf(_("%s: Wrong type for argument #%d: Real matrix expected.\n"), "mesh2d", i));
35         end
36     end
37     if (size(varargin(1),"*") <> size(varargin(2),"*"))
38         error(msprintf(_("%s: Arguments #%d and #%d must have the same sizes.\n"), "mesh2d", 1,2));
39     end
40     if (size(varargin(1),"*") < 3)
41         error(msprintf(_("%s: Wrong size for input argument #%d : At least %d elements expected.\n"), "mesh2d)", 1, 3));
42     end
43     if argn(2) == 2
44         [hull, area] = quickhull(varargin(1), varargin(2));
45         if area < %eps
46             error(msprintf(_("%s: all points are aligned.\n"),"mesh2d"))
47         end
48         tri = mesh2di(varargin(1), varargin(2), hull);
49     else
50         tri = mesh2di(varargin(1), varargin(2), varargin(3));
51     end
52 end
53
54 function [ind, area]=quickhull(x, y)
55     // https://fr.wikipedia.org/wiki/Quickhull
56     // Quickhull is a method of computing the convex hull of a finite set of points in the plane.
57     // It uses a divide and conquer approach similar to that of quicksort, from which its name derives.
58     // Its best case complexity is considered to be O(n * log(n))
59     // The following implementation supports all degenerates cases.
60
61     // Determine vertices of minimum and maximum abcissa
62     [minx,imin] = min(x(:));
63     [maxx,imax] = max(x(:));
64     if maxx-minx < %eps * (1 + max(maxx,minx))
65       [miny,imin] = min(y(:));
66       [maxy,imax] = max(y(:));
67       if maxy-miny < %eps * (1 + max(maxy,miny))
68         error(msprintf(_("%s: all points are coincident.\n"),"mesh2d"))
69       end
70     end  
71     
72     ind = 1:size(x,"*");
73     ind([imin;imax]) = [];
74     // Call recursive divide and conquer function
75     // 1-for vertices under segment [imin imax]
76     [ind1,remain] = conquer(x(:),y(:),ind,[imin;imax]);
77     if ~isempty(remain) then
78         // 2-for vertices (if any) above segment [imin imax]
79         ind2 = conquer(x(:),y(:),remain,[imax;imin]);
80         ind = [ind1(1:$-1);ind2];
81     else // if all vertices were under segment [imin imax]
82         ind=[ind1;imin];
83     end
84     area = sum(x(ind(1:$-1)).*y(ind(2:$))-x(ind(2:$)).*y(ind(1:$-1)))/2
85 end
86
87 function [out, remain]=conquer(x, y, ind, indH)
88     out = indH;
89     xh = x(indH);
90     yh = y(indH);
91     // compute signed distance w.r.t to segment [indH(1) indH(2)]
92     u=[xh(2)-xh(1);yh(2)-yh(1)];
93     d = (y(ind)-yh(1))*u(1)-(x(ind)-xh(1))*u(2);
94     // 1-divide
95     // ind(iplus) = vertices to the right of segment [indH(1) indH(2)]
96     iplus = find(d<=0);
97     remain = ind;
98     remain(iplus) = [];
99     if ~isempty(iplus) then
100         // #indnew = the furthest vertex to the right of segment [indH(1),indH(2)
101         [ma,ima]=min(d(iplus));
102         if ma > -%eps // all the righ-of-segment vertices are *on* the segment
103             iOnSeg = ind(iplus);
104             // we sort vertices by increasing L_1 distance from indH(1)
105             [sorted,iSorted] = gsort(abs(x(iOnSeg)-x(indH(1)))+...
106                                      abs(y(iOnSeg)-y(indH(1))),'g','i');
107             out = [indH(1);iOnSeg(iSorted)';indH(2)]
108             return
109         end
110         indnew = ind(iplus(ima));
111         iplus(ima) = [];
112         // 2-conquer
113         if ~isempty(iplus) then
114             // process ind(iplus) vertices w.r.t. segment [indH(1) indnew]
115             [ind1,remainplus]=conquer(x,y,ind(iplus),[indH(1);indnew]);
116             if ~isempty(remainplus) then
117                 // process remaining vertices w.r.t. segment [indnew indH(2)]
118                 ind2=conquer(x,y,remainplus,[indnew;indH(2)]);
119                 out = [ind1(1:$-1);ind2];
120             else
121                 out = [ind1;indH(2)];
122             end
123         else // #indnew vertex was alone to the right of [indH(1) indH(2)]
124             out = [indH(1);indnew;indH(2)];
125         end
126     end
127 end
128