diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m index 30803a6..5bbf4d5 100644 --- a/chunkie/@kernel/kernel.m +++ b/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 diff --git a/chunkie/chunkerkernevalmat.m b/chunkie/chunkerkernevalmat.m index 4169d53..08ccc15 100644 --- a/chunkie/chunkerkernevalmat.m +++ b/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 diff --git a/chunkie/chunkermat.m b/chunkie/chunkermat.m index f3fd7bd..68da513 100644 --- a/chunkie/chunkermat.m +++ b/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); diff --git a/chunkie/demo/demo_barbell.m b/chunkie/demo/demo_barbell.m index 96cc7ce..65965c8 100644 --- a/chunkie/demo/demo_barbell.m +++ b/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(); diff --git a/chunkie/roundedpolygon.m b/chunkie/roundedpolygon.m index 867e3ca..5ba0c08 100644 --- a/chunkie/roundedpolygon.m +++ b/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); \ No newline at end of file diff --git a/devtools/test/chunkerfuncTest.m b/devtools/test/chunkerfuncTest.m index 67de6ca..1820d76 100644 --- a/devtools/test/chunkerfuncTest.m +++ b/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)