Merge pull request #82 from fastalgorithms/issue69-target-regions

Issue69 target regions closes #69

chunkie/@chunkgraph/chunkgraph.m

@@ -107,6 +107,12 @@
                 if (numel(fchnks)<i || isempty(fchnks{i}))
                     i1 = edgesendverts(1,i);
                     i2 = edgesendverts(2,i);
+                    if (i1 == i2)
+                        msg = "chunkgraph constructor: connecting a " + ...
+                            "vertex to itself by a line " + ...
+                            "may have unexpected behavior";
+                        warning(msg);
+                    end
                     v1 = verts(:,i1);
                     v2 = verts(:,i2);
                     fcurve = @(t) chnk.curves.linefunc(t,v1,v2);
@@ -115,24 +121,37 @@
                     %chnkr.vert = [v1,v2];
                     echnks(i) = chnkr;
                 elseif (~isempty(fchnks{i}) && isa(fchnks{i},'function_handle'))
-                    [vs,~,~] =fchnks{i}([0,1]);
+                    vs =fchnks{i}([0,1]);
                     chnkr = chunkerfunc(fchnks{i},cploc,pref);
                     chnkr = sort(chnkr);
                     i1 = edgesendverts(1,i);
                     i2 = edgesendverts(2,i);
-                    vfin0 = verts(:,i1);
-                    vfin1 = verts(:,i2);
-                    r0 = vs(:,1);
-                    r1 = vfin0;
-                    scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro');
-                    xdfin = vfin1(1)-vfin0(1);
-                    ydfin = vfin1(2)-vfin0(2);
-                    tfin = atan2(ydfin,xdfin);
-                    xdini = vs(1,2)-vs(1,1);
-                    ydini = vs(2,2)-vs(2,1);
-                    tini = atan2(ydini,xdini);
-                    trotat = tfin - tini;
-                    chnkr = move(chnkr,r0,r1,trotat,scale);
+                    if i1 ~= i2
+                        vfin0 = verts(:,i1);
+                        vfin1 = verts(:,i2);
+                        r0 = vs(:,1);
+                        r1 = vfin0;
+                        scale = norm(vfin1-vfin0,'fro')/norm(vs(:,2)-vs(:,1),'fro');
+                        xdfin = vfin1(1)-vfin0(1);
+                        ydfin = vfin1(2)-vfin0(2);
+                        tfin = atan2(ydfin,xdfin);
+                        xdini = vs(1,2)-vs(1,1);
+                        ydini = vs(2,2)-vs(2,1);
+                        tini = atan2(ydini,xdini);
+                        trotat = tfin - tini;
+                        chnkr = move(chnkr,r0,r1,trotat,scale);
+                    else
+                        vfin = verts(:,i1);
+                        if norm(vs(:,1)-vs(:,2))/max(abs(chnkr.r(:))) > 1e-14
+                            msg = "chunkgraph constructor: edge " + ...
+                            num2str(i) + " is defined as a loop but " + ...
+                            "curve defining the edge does not reconnect " + ...
+                            "to high precision";
+                            warning(msg);
+                        end
+                        chnkr = move(chnkr,zeros(size(vs(:,1))),vfin-vs(:,1),0,1);
+                    end
+
                     echnks(i) = chnkr;
                 end
             end

chunkie/chunkerinterior.m

@@ -61,12 +61,14 @@
     grid = true;
     [xx,yy] = meshgrid(x,y);
     pts = [xx(:).'; yy(:).'];
-elseif isa(ptsobj, "chunker")
+elseif isa(ptsobj, "chunker") || isa(ptsobj, "chunkgraph") || ...
+        (isstruct(ptsobj) && isfield(ptsobj,"r"))
     pts = ptsobj.r(:,:);
-elseif isa(ptsobj, "chunkgraph")
-    pts = ptsobj.r(:,:);
-else
+elseif isnumeric(ptsobj)
     pts = ptsobj;
+else
+    msg = "chunkerinterior: input 2 not a recognized type";
+    error(msg);
 end
 
 

chunkie/chunkgraphinregion.m

@@ -0,0 +1,112 @@
+function [ids] = chunkgraphinregion(cg,ptsobj,opts)
+%CHUNKGRAPHINREGION returns an array indicating the region number for
+% each point specified by pts.
+%
+% Syntax: ids = chunkerinterior(cg,pts,opts)
+%         ids = chunkerinterior(cg,{x,y},opts) % meshgrid version
+%
+% Input:
+%   cg - chunkgraph object describing geometry
+%   ptsobj - object describing the target points, can be specified as
+%       * (cg.echnks(1).dim,:) array of points to test
+%       * {x,y} - length 2 cell array. the points checked then have the
+%           coordinates of a mesh grid [xx,yy] = meshgrid(x,y)
+%       * chunker object, in which case it uses ptsobj.r(:,:) 
+%       * chunkgraph object, in which case it uses ptsobj.r(:,:) 
+%
+% Optional input:
+%   opts - options structure with entries:
+%       opts.fmm = boolean, use FMM 
+%       opts.flam = boolean, use FLAM routines
+%       opts.axissym = boolean, chunker is axissymmetric
+%  Note on the default behavior: 
+%    by default it tries to use the fmm if it exists, if it doesn't
+%    then unless explicitly set to false, it tries to use flam
+%
+% Output:
+%   ids - integer array, pts(:,i) is in region ids(i)
+%
+% Examples:
+%   verts = [1 0 -1 0; 0 1 0 -1]; edgesendverts = [1:3, 3, 4; 2:3, 1, 4, 1];
+%   cg = chunkgraph(verts,edgesendverts);
+%   pts = 2*randn(2,100);
+%   ids = chunkgraphinregion(cg,pts);
+%
+% see also CHUNKGRAPH, CHUNKERINTERIOR
+
+% author: Travis Askham ([email protected])
+
+if nargin < 3
+    opts = [];
+end
+
+% Assign appropriate object to chnkr
+msg = "chunkgraphinregion: input 1 must be chunkgraph";
+assert(class(cg) == "chunkgraph",msg);
+
+% Figure out size of ids array based on ptsobj
+if isa(ptsobj, "cell")
+    assert(length(ptsobj)==2,'second input should be either 2xnpts array or length 2 cell array');
+    x = ptsobj{1};
+    y = ptsobj{2};
+    npts = length(x)*length(y);
+elseif isa(ptsobj, "chunker") || isa(ptsobj, "chunkgraph") || ...
+        (isstruct(ptsobj) && isfield(ptsobj,"r"))
+    npts = size(ptsobj.r(:,:),2);
+elseif isnumeric(ptsobj)
+    npts = size(ptsobj(:,:),2);
+else
+    msg = "chunkgraphinregion: input 2 not a recognized type";
+    error(msg);
+end
+
+ids = nan(npts,1);
+
+% loop over regions and use chunkerinterior to label
+% TODO: make a more efficient version 
+
+nedge0 = numel(cg.echnks);
+if nedge0 > 0
+    k = cg.echnks(1).k;
+    t = cg.echnks(1).tstor;
+    w = cg.echnks(1).wstor;
+    p = struct("k",k);
+else
+    ids = [];
+    return
+end
+
+nr = numel(cg.regions);
+for ir = 1:nr
+    ncomp = numel(cg.regions{ir});
+    intmp = zeros(npts,1);
+    for ic = 1:ncomp
+        edgelist = cg.regions{ir}{ic};
+        nedge = numel(edgelist);
+        chnkrs(nedge) = chunker(p,t,w);
+        for ie = 1:nedge
+            eid = edgelist(ie);
+            if eid > 0
+                if ir == 1
+                    chnkrs(ie) = cg.echnks(eid);
+                else
+                    chnkrs(ie) = reverse(cg.echnks(eid));
+                end
+            else
+                if ir == 1
+                    chnkrs(ie) = reverse(cg.echnks(-eid));
+                else
+                    chnkrs(ie) = cg.echnks(-eid);
+                end
+            end
+        end
+        intmp = intmp + reshape(chunkerinterior(merge(chnkrs(1:nedge)),ptsobj,opts),npts,1);
+    end
+
+    intmp = intmp > 0;
+    if ir == 1
+        ids(~intmp) = ir;
+    else
+        ids(intmp) = ir;
+    end
+end

devtools/test/chunkgraph_basicTest.m

@@ -38,11 +38,6 @@
 
 cgrph1 = balance(cgrph1);
 
-% adjacent triangles
-
-verts = [1 0 -1 0; 0 1 0 -1]; edgesendverts = [1:3, 3, 4; 2:3, 1, 4, 1];
-cg = chunkgraph(verts,edgesendverts);
-
 % a graph with two edges in the old format
 
 verts = [1 0 1; -1 0 1]; edge2verts = [-1 1 0; 0 -1 1];
@@ -53,6 +48,111 @@
 
 assert(nnz(cg1.v2emat-cg2.v2emat) == 0);
 
+% testing multiply connected
+
+verts = [1 0 -1 4 3 2; 0 1 0 0 1 0]; 
+edgends = [1 2 3 4 5 6; 2 3 1 5 6 4];
+cg1 = chunkgraph(verts,edgends);
+
+assert(numel(cg1.regions) == 3);
+
+% testing connected by bridge
+
+verts = [1 0 -1 4 3 2; 0 1 0 0 1 0]; 
+edgends = [1 2 3 4 5 6 1; 3 1 2 6 4 5 6];
+cg1 = chunkgraph(verts,edgends);
+
+assert(numel(cg1.regions) == 3)
+
+% testing a loop
+
+verts = [2;1]; edgends = [1;1];
+cfun = @(t) [cos(2*pi*t(:).'); sin(2*pi*t(:).').*sin(pi*t(:).')];
+fchnks = {cfun};
+cg1 = chunkgraph(verts,edgends,fchnks);
+
+assert(numel(cg1.regions) == 2)
+
+% testing nested 
+
+verts = [1 0 -1 2 0 -2; 0 1 0 -1 2 -1];
+edgends = [1 2 3 4 5 6; 2 3 1 5 6 4];
+cg1 = chunkgraph(verts,edgends);
+
+assert(numel(cg1.regions) == 3)
+
+% region id for targets: compare to inpolygon routine for polygonal domains
+
+% adjacent triangles
+
+verts = [1 0 -1 0; 0 1 0 -1]; edgesendverts = [1:3, 3, 4; 2:3, 1, 4, 1];
+cg = chunkgraph(verts,edgesendverts);
+
+x1 = linspace(-pi,pi); 
+[xx,yy] = meshgrid(x1,x1);
+
+targs = [xx(:).'; yy(:).'];
+
+ids = chunkgraphinregion(cg,targs);
+ids2 = chunkgraphinregion(cg,{x1,x1});
+
+idstrue = polygonids(cg,xx,yy);
+
+assert(nnz(ids(:)-idstrue(:)) == 0)
+assert(nnz(ids(:)-ids2(:)) == 0)
+
+% nested triangles
+
+verts = [1 0 -1 2 0 -2; 0 1 0 -1 2 -1];
+edgends = [1 2 3 4 5 6; 2 3 1 5 6 4];
+cg = chunkgraph(verts,edgends);
+
+x1 = linspace(-pi,pi); % avoid edge cases...
+[xx,yy] = meshgrid(x1,x1);
+
+targs = [xx(:).'; yy(:).'];
+
+ids = chunkgraphinregion(cg,targs);
+
+idstrue = polygonids(cg,xx,yy);
+
+assert(nnz(ids(:)-idstrue(:)) == 0)
+
+function idstrue = polygonids(cg,xx,yy)
+
+verts = cg.verts;
+edgesendverts = cg.edgesendverts;
+
+idstrue = nan(size(xx(:)));
+regions = cg.regions;
+for j = 1:numel(regions)
+    vlist = [];
+    for ic = 1:numel(regions{j})
+        elist = regions{j}{ic};
+        for ie = 1:numel(elist)
+            eid = elist(ie);
+            if eid > 0
+                v1 = edgesendverts(1,eid);
+                v2 = edgesendverts(2,eid);
+            else
+                v1 = edgesendverts(2,-eid);
+                v2 = edgesendverts(1,-eid);
+            end
+            vlist = [vlist, v1];
+            if ie == numel(elist)
+                vlist = [vlist, v2];
+            end
+        end
+    end
+    intmp = inpolygon(xx(:),yy(:),verts(1,vlist),verts(2,vlist));
+    if j == 1
+        idstrue(intmp == 0) = j;
+    else
+        idstrue(intmp > 0) = j;
+    end
+end
+
+end
 
 function [r,d,d2] = sinearc(t,amp,frq)