point roundedpolygon back to chunkerpoly, fixes chnkr.h bug

chunkie/@kernel/kernel.m

@@ -85,17 +85,21 @@
           elseif ( isa(kern, 'function_handle') )
               obj.eval = kern;
           elseif ( isa(kern, 'kernel') )
-              if ( numel(kern) == 1 )
-                  obj = kern;
-              else
-                  % The input is a matrix of kernels.
-                  % Create a single kernel object by interleaving the
-                  % outputs of each sub-kernel's eval() and fmm() routines.
-                  % TODO: Check that opdims are consistent
-                  obj = interleave(kern);
+              if (numel(kern) == 1)
+		obj = kern;
+	      else
+		obj = interleave(kern);
               end
+          elseif isa(kern,'cell')
+            sz = size(kern); assert(length(sz)==2,'KERNEL: first input not of a supported type');
+            obj(sz(1),sz(2)) = kernel();
+            for j = 1:sz(2)
+                for i = 1:sz(1)
+                    obj(i,j) = kernel(kern{i,j});
+                end
+            end
           else
-              error('First input must be a string or function handle.');
+              error('KERNEL: first input not of a supported type');
           end
 
       end

chunkie/chunkerkernevalmat.m

@@ -1,4 +1,4 @@
-function mat = chunkerkernevalmat(chnkr,kern,targobj,opts)
+    function mat = chunkerkernevalmat(chnkr,kern,targobj,opts)
 %CHUNKERKERNEVALMAT compute the matrix which maps density values on 
 % the chunk geometry to the value of the convolution of the given
 % integral kernel with the density at the specified target points

chunkie/chunkermat.m

@@ -92,33 +92,28 @@
 % opts.sing provides a default value for singularities if not 
 % defined for kernels
 
-if isa(kern,'function_handle')
-    kern2 = kernel(kern);
-    kern = kern2;
-elseif isa(kern,'cell')
-    sz = size(kern);
-    kern2(sz(1),sz(2)) = kernel();
-    for j = 1:sz(2)
-        for i = 1:sz(1)
-            if isa(kern{i,j},'function_handle')
-                kern2(i,j) = kernel(kern{i,j});
-            elseif isa(kern{i,j},'kernel')
-                kern2(i,j) = kern{i,j};
-            else
-                msg = "Second input is not a kernel object, function handle, " ...
-                    + "or cell array";
-                error(msg);
-            end
-        end
+
+% Flag for determining whether input object is a chunkergraph
+icgrph = 0;
+
+if (class(chnkobj) == "chunker")
+    chnkrs = chnkobj;
+elseif(class(chnkobj) == "chunkgraph")
+    icgrph = 1;
+    chnkrs = chnkobj.echnks;
+else
+    msg = "CHUNKERMAT: first input is not a chunker or chunkgraph object";
+    error(msg)
+end
+
+if ~isa(kern,'kernel')
+    try 
+        kern = kernel(kern);
+    catch
+        error('CHUNKERMAT: second input kern not of supported type');
     end
-    kern = kern2;
-
-elseif ~isa(kern,'kernel')
-    msg = "Second input is not a kernel object, function handle, " ...
-                + "or cell array";
-    error(msg);
 end
-    
+
 if nargin < 3
     opts = [];
 end
@@ -167,18 +162,6 @@
     adaptive_correction = opts.adaptive_correction;
 end
 
-% Flag for determining whether input object is a chunkergraph
-icgrph = 0;
-
-if (class(chnkobj) == "chunker")
-    chnkrs = chnkobj;
-elseif(class(chnkobj) == "chunkgraph")
-    icgrph = 1;
-    chnkrs = chnkobj.echnks;
-else
-    msg = "First input is not a chunker or chunkgraph object";
-    error(msg)
-end
 
 nchunkers = length(chnkrs);
 

chunkie/demo/demo_barbell.m

@@ -21,12 +21,13 @@
 % parameters for curve rounding/chunking routine to oversample boundary
 cparams = [];
 cparams.widths = 0.1*ones(size(verts,2),1);% width to cut about each corner
-cparams.eps = 1e-6; % tolerance at which to resolve curve
+cparams.eps = 1e-12; % tolerance at which to resolve curve
 cparams.rounded = true;
 
 % call smoothed-polygon chunking routine
 % a smoothed version of edgevals is returned in 
 % chnkr.data
+
 chnkr = chunkerpoly(verts,cparams,[],edgevals);
 chnkr = chnkr.refine();
 

chunkie/roundedpolygon.m

@@ -53,210 +53,4 @@
 %
 % See also CHUNKERFUNC, CHUNKER, CHUNKERPREF, REFINE
 
-[dimv,nv] = size(verts);
-
-assert(dimv > 1)
-
-if nargin < 2
-    cparams = [];
-end
-if nargin < 3
-    p = [];
-    p.dim = dimv;
-    pref = chunkerpref(p);
-else
-    pref = chunkerpref(pref);
-    if pref.dim ~= dimv
-        warning('dimensions dont match, overwriting with vertex dim');
-        pref.dim = dimv;
-    end
-end
-if nargin < 4
-    edgevals = [];
-end
-
-autowidths = false;
-autowidthsfac = 0.1;
-ifclosed = true;
-eps = 1e-6;
-nover = 0;
-
-if isfield(cparams,'ifclosed')
-   ifclosed = cparams.ifclosed;
-end
-if isfield(cparams,'eps')
-   eps = cparams.eps;
-end
-if isfield(cparams,'nover')
-   nover = cparams.nover;
-end
-
-if (ifclosed)
-   verts2 = [verts(:,end), verts, verts(:,1)];
-   edges2 = sqrt(sum(diff(verts2,1,2).^2,1));
-else
-   edges2 = [0, sqrt(sum(diff(verts,1,2).^2,1)), 0];
-end
-
-widths_not_set = true;
-
-if isfield(cparams,'widths')
-   widths = cparams.widths;
-   widths_not_set = false;
-end
-if isfield(cparams,'autowidths')
-   autowidths = cparams.autowidths;
-end
-if isfield(cparams,'autowidthsfac')
-   autowidthsfac = cparams.autowidthsfac;
-end
-
-if (autowidths || widths_not_set)
-   widths = autowidthsfac*...
-   	  min(edges2(1:end-1),edges2(2:end));
-end
-
-if ifclosed
-    widths = [widths(:); widths(1)];
-    verts = [verts, verts(:,1)];
-    nv = size(verts,2);
-end
-
-nedge = nv - 1;
-nvals = 0;
-if numel(edgevals) ~= 0
-    assert(rem(numel(edgevals),nedge) == 0, ...
-        'number of edge values should be multiple of number of edges');
-    nvals = numel(edgevals)/nedge;
-    edgevals = reshape(edgevals,nvals,nedge);
-end
-
-chnkr = chunker(pref);
-chnkr = chnkr.makedatarows(nvals);
-k = chnkr.k; dim = chnkr.dim;
-[t,w] = lege.exps(k);
-onek = ones(k,1);
-
-
-if nvals > 0
-    aint = lege.intmat(chnkr.k);
-end
-
-for i = 1:nv-1
-    % grab vertices
-    r1 = verts(:,i); r2 = verts(:,i+1);
-    w1 = widths(i); w2 = widths(i+1);
-    l = sqrt(sum((r1-r2).^2));
-    assert(l > w1+w2+2*eps(1)*l,'widths too large for side');
-
-    % make chunk in middle
-    v = (r2-r1)/l;
-    ts = w1+(l-w2-w1)*(t+1)/2.0;
-    chnkr = chnkr.addchunk();
-    nch = chnkr.nch;
-    if nvals > 0
-        val1 = edgevals(:,i);
-        chnkr.data(:,:,nch) = bsxfun(@times,val1,onek(:).');
-    end
-    chnkr.r(:,:,nch) = r1 + bsxfun(@times,v(:),(ts(:)).');
-    chnkr.d(:,:,nch) = repmat(v(:),1,k);
-    chnkr.d2(:,:,nch) = zeros(dim,k);
-    chnkr.adj(1,nch) = 0; chnkr.adj(2,nch) = 0;
-    if nch > 1
-        chnkr.adj(1,nch) = nch-1;
-        chnkr.adj(2,nch-1) = nch;
-    end
-    chnkr.h(nch) = (l-w2-w1)/2.0;
-
-    if or(i < nv-1,ifclosed)
-			% chunk up smoothed corner made by three verts
-      if (i==nv-1)
-        if nvals > 0
-          val2 = edgevals(:,1);
-        end
-        r3 = verts(:,2);
-      else
-        if nvals > 0
-          val2 = edgevals(:,i+1);
-        end
-        r3 = verts(:,i+2);
-      end
-      l2 = sqrt(sum((r2-r3).^2));
-      v = -v;
-      v2 = (r3-r2)/l2;
-      cosphi = dot(v2,v); %cosine of angle between edges
-      sinphi = sqrt(1-cosphi*cosphi);
-      trange = w2*sqrt((1-cosphi)/2); %parameter space range of gaussian
-      hbell = abs(trange)/8.0; % width parameter of gaussian
-      m = sinphi/(1-cosphi); % slope of abs approx
-      cpt.ta = -trange;
-      cpt.tb = trange;
-      cpt.eps = eps; cpt.levrestr = 0; cpt.ifclosed = 0;
-      chnkrt = sort(chunkerfunc(@(t)fround(t,m,hbell,dim),cpt,pref));
-
-		  % do optimal procrustes match of left and right ends
-      rl = fround(-trange,m,hbell,dim); rr = fround(trange,m,hbell,dim);
-      [um,~,vm] = svd([w2*v w2*v2]*([rl rr].'));
-      rotmat = um*vm.';
-
-			       % copy in rotated and translated chunks
-      ncht = chnkrt.nch;
-      chnkr = chnkr.addchunk(ncht);
-      chnkr.r(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.r(:,:),dim,k,ncht)+r2;
-      chnkr.d(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d(:,:),dim,k,ncht);
-      chnkr.d2(:,:,nch+1:nch+ncht) = reshape(rotmat*chnkrt.d2(:,:),dim,k,ncht);
-      chnkr.adj(:,nch+1:nch+ncht) = chnkrt.adj+nch;
-      chnkr.adj(2,nch) = nch+1;
-      chnkr.adj(1,nch+1) = nch;
-      chnkr.h(nch+1:nch+ncht) = chnkrt.h;
-
-      if nvals > 0
-        ds = bsxfun(@times,reshape(sum((rotmat*chnkrt.d(:,:)).^2,1),k,ncht), ...
-                    (chnkrt.h(:)).');
-        dsw = bsxfun(@times,w(:),ds);
-        dssums = sum(dsw,1);
-        dssums2 = cumsum([0,dssums(1:end-1)]);
-        dsint = aint*ds;
-        dsint = bsxfun(@plus,dsint,dssums2);
-        lencorner = sum(dsw(:));
-        ss = -lencorner/2.0 + dsint;
-        ss = ss/lencorner*16;
-        erfss = erf(ss);
-        datass = reshape((val2(:)-val1(:))/2*((erfss(:)).'+1) ...
-			 +val1(:),nvals,k,ncht);
-        chnkr.data(:,:,nch+1:nch+ncht) = datass;
-      end
-    end
-
-end
-
-if ifclosed
-    nch = chnkr.nch;
-    chnkr.adj(1,1) = nch;
-    chnkr.adj(2,nch) = 1;
-end
-
-chnkr = chnkr.refine(cparams);
-
-% compute normals
-
-chnkr.n = normals(chnkr);
-
-end
-
-function [r,d,d2] = fround(t,m,h,dim)
-
-[y,dy,d2y] = chnk.spcl.absconvgauss(t,m,0.0,h);
-
-r = zeros(dim,length(t));
-d = zeros(dim,length(t));
-d2 = zeros(dim,length(t));
-
-r(1,:) = t;
-r(2,:) = y;
-d(1,:) = 1.0;
-d(2,:) = dy;
-d2(1,:) = 0.0;
-d2(2,:) = d2y;
-
-end
+chnkr = chunkerpoly(verts,cparams,pref,edgevals);

devtools/test/chunkerfuncTest.m

@@ -83,6 +83,9 @@
 
 a = area(chnkr);
 assert(abs(a - pi*r^2) < 1e-12,'area wrong for circle domain')
+chnkr = refine(chnkr,struct('nover',1));
+a = area(chnkr);
+assert(abs(a - pi*r^2) < 1e-12,'area wrong for circle domain')
 
 
 % subplot(1,3,3)