From d1b8c8b3dd2ef1714c075eeaad14518caac72b36 Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Fri, 31 May 2024 11:25:49 -0400 Subject: [PATCH 1/2] chunkgraphinregion routine for labeling targets by region. fixes some minor bugs in chunkgraph constructor and adds appropriate tests --- chunkie/@chunkgraph/chunkgraph.m | 47 ++++++++---- chunkie/chunkgraphinregion.m | 111 +++++++++++++++++++++++++++ devtools/test/chunkgraph_basicTest.m | 108 ++++++++++++++++++++++++-- 3 files changed, 247 insertions(+), 19 deletions(-) create mode 100644 chunkie/chunkgraphinregion.m 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/chunkgraphinregion.m b/chunkie/chunkgraphinregion.m new file mode 100644 index 0000000..bd49c07 --- /dev/null +++ b/chunkie/chunkgraphinregion.m @@ -0,0 +1,111 @@ +function [ids] = chunkgraphinregion(chnkobj,ptsobj,opts) +%CHUNKERINTERIOR returns an array indicating whether each point specified +% % by pts is inside the domain. Assumes the domain is closed. +% +% Syntax: in = chunkerinterior(chnkobj,pts,opts) +% in = chunkerinterior(chnkobj,{x,y},opts) % meshgrid version +% +% Input: +% chnkobj - chunker object or chunkgraph object describing geometry +% ptsobj - object describing the target points, can be specified as +% * (chnkr.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 chunker.r(:,:) +% * chunkgraph object, in which case it uses chunkgraph.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: +% in - logical array, if in(i) is true, then pts(:,i) is inside the +% domain or for a mesh grid [xx(i); yy(i)] is inside the domain. +% +% Examples: +% chnkr = chunkerfunc(@(t) starfish(t)); +% pts = 2*randn(2,100); +% in = chunkerinterior(chnkr,pts); +% + +% 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(chnkobj) == "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(chnkobj.echnks); +if nedge0 > 0 + k = chnkobj.echnks(1).k; + t = chnkobj.echnks(1).tstor; + w = chnkobj.echnks(1).wstor; + p = struct("k",k); +else + ids = []; + return +end + +nr = numel(chnkobj.regions); +for ir = 1:nr + ncomp = numel(chnkobj.regions{ir}); + intmp = zeros(npts,1); + for ic = 1:ncomp + edgelist = chnkobj.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) = chnkobj.echnks(eid); + else + chnkrs(ie) = reverse(chnkobj.echnks(eid)); + end + else + if ir == 1 + chnkrs(ie) = reverse(chnkobj.echnks(-eid)); + else + chnkrs(ie) = chnkobj.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..dfc2016 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,109 @@ 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); + +idstrue = polygonids(cg,xx,yy); + +assert(nnz(ids(:)-idstrue(:)) == 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) From 6d207a0b926c56f0644690247c9d890765afc4f9 Mon Sep 17 00:00:00 2001 From: Travis Askham Date: Fri, 31 May 2024 11:35:01 -0400 Subject: [PATCH 2/2] update docs, error handling in chunkerinterior --- chunkie/chunkerinterior.m | 10 +++--- chunkie/chunkgraphinregion.m | 51 ++++++++++++++-------------- devtools/test/chunkgraph_basicTest.m | 2 ++ 3 files changed, 34 insertions(+), 29 deletions(-) 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 index bd49c07..bce5060 100644 --- a/chunkie/chunkgraphinregion.m +++ b/chunkie/chunkgraphinregion.m @@ -1,18 +1,18 @@ -function [ids] = chunkgraphinregion(chnkobj,ptsobj,opts) -%CHUNKERINTERIOR returns an array indicating whether each point specified -% % by pts is inside the domain. Assumes the domain is closed. +function [ids] = chunkgraphinregion(cg,ptsobj,opts) +%CHUNKGRAPHINREGION returns an array indicating the region number for +% each point specified by pts. % -% Syntax: in = chunkerinterior(chnkobj,pts,opts) -% in = chunkerinterior(chnkobj,{x,y},opts) % meshgrid version +% Syntax: ids = chunkerinterior(cg,pts,opts) +% ids = chunkerinterior(cg,{x,y},opts) % meshgrid version % % Input: -% chnkobj - chunker object or chunkgraph object describing geometry +% cg - chunkgraph object describing geometry % ptsobj - object describing the target points, can be specified as -% * (chnkr.dim,:) array of points to test +% * (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 chunker.r(:,:) -% * chunkgraph object, in which case it uses chunkgraph.r(:,:) +% * 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: @@ -24,14 +24,15 @@ % then unless explicitly set to false, it tries to use flam % % Output: -% in - logical array, if in(i) is true, then pts(:,i) is inside the -% domain or for a mesh grid [xx(i); yy(i)] is inside the domain. +% ids - integer array, pts(:,i) is in region ids(i) % % Examples: -% chnkr = chunkerfunc(@(t) starfish(t)); +% 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); -% in = chunkerinterior(chnkr,pts); +% ids = chunkgraphinregion(cg,pts); % +% see also CHUNKGRAPH, CHUNKERINTERIOR % author: Travis Askham (askhamwhat@gmail.com) @@ -41,7 +42,7 @@ % Assign appropriate object to chnkr msg = "chunkgraphinregion: input 1 must be chunkgraph"; -assert(class(chnkobj) == "chunkgraph",msg); +assert(class(cg) == "chunkgraph",msg); % Figure out size of ids array based on ptsobj if isa(ptsobj, "cell") @@ -64,38 +65,38 @@ % loop over regions and use chunkerinterior to label % TODO: make a more efficient version -nedge0 = numel(chnkobj.echnks); +nedge0 = numel(cg.echnks); if nedge0 > 0 - k = chnkobj.echnks(1).k; - t = chnkobj.echnks(1).tstor; - w = chnkobj.echnks(1).wstor; + 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(chnkobj.regions); +nr = numel(cg.regions); for ir = 1:nr - ncomp = numel(chnkobj.regions{ir}); + ncomp = numel(cg.regions{ir}); intmp = zeros(npts,1); for ic = 1:ncomp - edgelist = chnkobj.regions{ir}{ic}; + 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) = chnkobj.echnks(eid); + chnkrs(ie) = cg.echnks(eid); else - chnkrs(ie) = reverse(chnkobj.echnks(eid)); + chnkrs(ie) = reverse(cg.echnks(eid)); end else if ir == 1 - chnkrs(ie) = reverse(chnkobj.echnks(-eid)); + chnkrs(ie) = reverse(cg.echnks(-eid)); else - chnkrs(ie) = chnkobj.echnks(-eid); + chnkrs(ie) = cg.echnks(-eid); end end end diff --git a/devtools/test/chunkgraph_basicTest.m b/devtools/test/chunkgraph_basicTest.m index dfc2016..0bb2bfc 100644 --- a/devtools/test/chunkgraph_basicTest.m +++ b/devtools/test/chunkgraph_basicTest.m @@ -94,10 +94,12 @@ 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