diff --git a/chunkie/@chunkgraph/chunkgraph.m b/chunkie/@chunkgraph/chunkgraph.m index 6ec0053..b4a2679 100644 --- a/chunkie/@chunkgraph/chunkgraph.m +++ b/chunkie/@chunkgraph/chunkgraph.m @@ -107,6 +107,12 @@ if (numel(fchnks) 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 diff --git a/chunkie/chunkerinterior.m b/chunkie/chunkerinterior.m index f0c29fe..f33ea59 100644 --- a/chunkie/chunkerinterior.m +++ b/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 diff --git a/chunkie/chunkgraphinregion.m b/chunkie/chunkgraphinregion.m new file mode 100644 index 0000000..bce5060 --- /dev/null +++ b/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 (askhamwhat@gmail.com) + +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 diff --git a/devtools/test/chunkgraph_basicTest.m b/devtools/test/chunkgraph_basicTest.m index d771081..0bb2bfc 100644 --- a/devtools/test/chunkgraph_basicTest.m +++ b/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)