Merge pull request #117 from fastalgorithms/patch-chunkerinterior

Patch chunkerinterior

Author: mrachh

chunkie/@chunker/centroids.m

@@ -6,10 +6,10 @@
 k = chnkr.k;
 nch = chnkr.nch;
 dim = chnkr.dim;
-[~,w] = lege.exps(k);
+w = chnkr.wstor;
 ww = repmat( (w(:)).',dim,1);
 wall = reshape(ww(:)*ones(1,nch),dim,k,nch);
 
 ctrs = reshape(sum( (chnkr.r).*wall , 2),dim,nch)/2;
 
-end
+end

chunkie/@chunker/flagnear_rectangle.m

@@ -4,7 +4,7 @@
 % logical array. If the (i,j) entry is non zero then pts(:,i) is close
 % to chunk j in chnkr.
 %
-% Syntax: flag = flagnear(chnkr,pts,opts)
+% Syntax: flag = flagnear_rectangle(chnkr,pts,opts)
 %
 % Input:
 %   chnkr - chunker object describing curve

chunkie/chunkerinterior.m

@@ -124,27 +124,16 @@
    useflam_final = useflam;
 end
 
-% use bernstein ellipses and rectangles to flag problematic points
-%
 
 eps_local = 1e-3;
-
-rho = 1.2;
-optsflag = [];  optsflag.rho = rho; optsflag.occ = 5;
-if grid
-    flag = flagnear_rectangle_grid(chnkr,x,y,optsflag);
-else
-    flag = flagnear_rectangle(chnkr,pts,optsflag);
-end
-
+rho = 1.6;
 npoly = chnkr.k;
 nlegnew = chnk.ellipse_oversample(rho,npoly,eps_local);
 nlegnew = max(nlegnew,chnkr.k);
 
 [chnkr2] = upsample(chnkr,nlegnew);
 
 
-
 icont = false;
 if usefmm_final
    try
@@ -185,14 +174,23 @@
 % for points where the integral might be inaccurate:
 % find close boundary point and check normal direction
 
+iffy = min(abs(vals1+1),abs(vals1)) > 1e-2;
+
+ipt = find(iffy(:));
+pts_iffy = pts(:,iffy);
+
+optsflag = [];  optsflag.rho = rho; optsflag.occ = 5;
+flag = flagnear_rectangle(chnkr,pts_iffy,opts);
+flag2 = flagnear_rectangle(chnkr,chnkr.r(:,:),opts);
+flag2 = ((flag2.'*kron(speye(chnkr.nch),ones(chnkr.k,1))) > 1).';
 
-nnzpt = sum(flag~=0,2);
-ipt = find(nnzpt);
+flag = (flag2*flag.').';
 
-npts = numel(pts)/2;
-distmins = inf(npts,1);
-dss = zeros(2,npts);
-rss = zeros(2,npts);
+npts_iffy = numel(pts_iffy)/2;
+assert(npts_iffy == length(ipt));
+distmins = inf(npts_iffy,1);
+dss = zeros(2,npts_iffy);
+rss = zeros(2,npts_iffy);
 
 k = chnkr.k;
 [t,~,u] = lege.exps(k);
@@ -204,7 +202,8 @@
     dval = chnkr.d(:,:,i);
     nval = chnkr.n(:,:,i);
     [ji] = find(flag(:,i));
-    ptsi = pts(:,ji);
+
+    ptsi = pts_iffy(:,ji);
     nptsi = size(ptsi,2);
     dist2all = reshape(sum( abs(reshape(ptsi,2,1,nptsi) ...
                     - reshape(rval,2,k,1)).^2, 1),k,nptsi);
@@ -223,9 +222,9 @@
     diffs = rval(:,ipti)-ptsi;
     dots = sum(ptsn.*diffs,1);
 
-    jsus = abs(dots) < 2e-1*sqrt(dist2all);
+    jsus = abs(dots) < 0.9*sqrt(dist2all);
     jii = ji(jsus);
-    [~,rs,ds,~,dist2s] = chnk.chunk_nearparam(rval,pts(:,jii),[],t,u);
+    [~,rs,ds,~,dist2s] = chnk.chunk_nearparam(rval,pts_iffy(:,jii),[],t,u);
     for j = 1:length(jii)
         jj = jii(j);
         if dist2s(j) < distmins(jj)
@@ -237,6 +236,8 @@
 end
 
 for i = 1:length(ipt)
-    jj = ipt(i);
-    in(jj) = (rss(:,jj)-pts(:,jj)).'*[dss(2,jj);-dss(1,jj)] > 0;
+    if distmins(i) < inf
+        jj = ipt(i);
+        in(jj) = (rss(:,i)-pts_iffy(:,i)).'*[dss(2,i);-dss(1,i)] > 0;
+    end
 end

devtools/test/chunkerinteriorTest.m

@@ -2,10 +2,8 @@
 % inside a domain or not
 % 
 
-clearvars; close all;
 seed = 8675309;
 rng(seed);
-addpaths_loc();
 
 doadap = false;
 
@@ -89,4 +87,46 @@
 targs = starfish(ttarg).*scal;
 
 in = chunkerinterior(chnkr,targs,struct('axissym',true));
-assert(all(in(:) == (scal(:) <= 1)));
+assert(all(in(:) == (scal(:) <= 1)));
+
+%% 
+
+% a stress test. previously this triggered a bug for the old version which
+% used the normals to determine interior points
+
+amp = 0.25;
+scale = .3;
+
+ctr = [-2;-1.6];
+chnkr_int= chunkerfunc(@(t) starfish(t,3,amp,ctr,pi/4,scale)); 
+chnkr_int = sort(reverse(chnkr_int));
+
+a = max(vecnorm(chnkr_int.r(:,:)))*1.01;
+chnkr_ext = chunkerfunc(@(t) [a*cos(t(:).'); a*sin(t(:).')]); 
+chnkr = merge([chnkr_ext,chnkr_int]);
+
+% return
+
+L = max(abs(chnkr.r),[],"all");
+x1 = linspace(-L,L,100);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).']; 
+xv=[chnkr_ext.r(1,:),nan,chnkr_int.r(1,:)];
+yv=[chnkr_ext.r(2,:),nan,chnkr_int.r(2,:)];
+%tic; in0 = inpolygon(xx(:),yy(:),xv,yv); toc
+tt = atan2(yy-ctr(2),xx-ctr(1)) + 2*pi;
+st = starfish(tt(:),3,amp,[0;0],pi/4,scale);
+ss2 = reshape(st(1,:).^2 + st(2,:).^2,size(xx));
+in0 = and((xx.^2+yy.^2)<a^2, ...
+    (xx-ctr(1)).^2 + (yy-ctr(2)).^2 > ss2);
+in0 = in0(:);
+tic; in = chunkerinterior(chnkr,{x1,x1}); toc
+
+% clf
+% plot(chnkr,'g-o')
+% hold on
+% scatter(xx(in(:)),yy(in(:)),'bo')
+% scatter(xx(~in(:)),yy(~in(:)),'rx')
+
+assert(all(in0 == in));
+

devtools/test/mixedbcTest.m

@@ -93,14 +93,9 @@
 targets = zeros(2,length(xxtarg(:)));
 targets(1,:) = xxtarg(:); targets(2,:) = yytarg(:);
 
-start = tic; in1 = chunkerinterior(merge(cgrph.echnks(edir)),{xtarg,ytarg});
-in2 = chunkerinterior(reverse(merge(cgrph.echnks(eneu))),{xtarg,ytarg}); 
-t1 = toc(start);
-in = in1 & ~in2;
-out = ~in;
-
 ids= chunkgraphinregion(cgrph,{xtarg,ytarg});
-nnz(in-(ids==2))
+in = ids == 2;
+in2 = ids == 3;
 
 fprintf('%5.2e s : time to find points in domain\n',t1)