Issue106/surf smooth

chunkie/+chnk/+smoother/expeval.m

@@ -0,0 +1,76 @@
+function [cvals] = expeval(czs,ifexps)
+
+    sz = size(czs);
+    czs = czs(:);
+    cvals= zeros(size(czs));
+
+    ifloc = find(abs(czs) <= 2);
+    iffar = find(abs(czs) >2);
+
+    if (numel(ifloc))
+    cloc = expeval_tayl(czs(ifloc));
+    cvals(ifloc) = cloc;
+    end
+    if (numel(iffar)>0)
+    cfar = expeval_lentz(czs(iffar));
+    cvals(iffar) = cfar;
+    end
+
+    if ((nargin>1) && ifexps)
+        cvals = cvals.*exp(-czs);
+    end
+
+    cvals = reshape(cvals,sz);
+
+end
+
+
+function [cvals] = expeval_tayl(czs)
+
+        gamma = 0.577215664901532860606512090082;
+
+        cvals = -gamma-log(czs);
+        cpow = czs;
+        rfac = 1;
+
+        eps = 1E-24;
+
+        for ii=1:31
+            cvals = cvals+cpow.*rfac/ii;
+            cpow = -czs.*cpow;
+            rfac = rfac/(ii+1);
+            % if (max(abs(cpow.*rfac)) < eps)
+            %     break;
+            % end
+        end
+
+         cvals = cvals.*exp(czs);
+
+end
+
+
+function [cvals] = expeval_lentz(czs)
+
+        nterms = 40;
+        cf = 1./czs;
+        cd = cf;
+        cv = 1/1E-30;
+
+        for ii=1:80
+        can = ii;
+        cd  = 1./(1+ii*cd);
+        cv  = 1 + ii./cv;
+        del = cd.*cv;
+        cf  = cf.*del;
+        cd  = 1./(czs+ii*cd);
+        cv  = czs + ii./cv;
+        del = cd.*cv;
+        cf  = cf.*del;
+        % if (max(abs(del-1)) < 1E-14)
+        %     break;
+        % end
+        end
+
+        cvals = cf;
+
+end

chunkie/+chnk/+smoother/get_mesh.m

@@ -0,0 +1,40 @@
+function qmesh = get_mesh(umesh, nchs, k)
+    [x, w] = lege.exps(k);
+
+    nchtot = sum(nchs(:));
+    n = nchtot*k;
+    qmesh.r = zeros(2,n);
+    qmesh.n = zeros(2,n);
+    qmesh.pseudo_normals = zeros(2,n);
+    qmesh.wts = zeros(n,1);
+
+    [~, ne] = size(umesh.verts);
+    verts_ext = [umesh.verts, umesh.verts(:,1)];
+    pseudo_normals_ext = [umesh.pseudo_normals, umesh.pseudo_normals(:,1)];
+    istart = 1;
+    for i=1:ne
+        [xext, wext] = get_nodes(x, w, nchs(i));
+        ll = length(xext);
+        iend = istart+ll-1;
+        qmesh.r(:,istart:iend) = verts_ext(:,i) + ...
+             (verts_ext(:,i+1) - verts_ext(:,i)).*xext.';
+        qmesh.n(:,istart:iend) = repmat(umesh.face_normals(:,i),[1,ll]);
+        qmesh.pseudo_normals(:,istart:iend) = ...
+             repmat(umesh.pseudo_normals(:,i), [1,ll]) + ...
+             (pseudo_normals_ext(:,i+1) - pseudo_normals_ext(:,i)).*xext.';
+        qmesh.wts(istart:iend) = wext.*umesh.lengths(i);
+        istart = istart + ll;
+    end
+
+end
+
+
+function [xext, wext] = get_nodes(x, w, nch)
+    tabs = 0:1/nch:1;
+    tas = tabs(1:nch);
+    tbs = tabs(2:end);
+    xext = (x(:)+1)/2.*(tbs - tas) + tas;
+    xext = xext(:);
+    wext = w(:)/2.*(tbs-tas);
+    wext = wext(:);    
+end

chunkie/+chnk/+smoother/get_sigs.m

@@ -0,0 +1,40 @@
+function [val,grad] = get_sigs(umesh,rt,sig0,dlam)
+
+        rc = umesh.centroids;
+        dlens = umesh.lengths;
+
+        [~,ns] = size(rc);
+        [~,nt] = size(rt);
+
+        xc = repmat(rc(1,:),nt,1);
+        yc = repmat(rc(2,:),nt,1);
+        dl = repmat(dlens,nt,1);
+
+        xt = repmat(rt(1,:).',1,ns);
+        yt = repmat(rt(2,:).',1,ns);
+
+        dx = xt-xc;
+        dy = yt-yc;
+        dr = sqrt(dx.^2+dy.^2);
+        dmin = min(dr,[],2);
+
+        dsigj = dl/dlam;
+        dexpn = exp(-(dr.^2-dmin.^2)./(2*sig0^2));
+        dsignum = sum(dsigj.*dexpn,2);
+        dsigden = sum(dexpn,2);
+        dsignumx = -sum(dsigj.*dexpn.*(dx./sig0^2),2);
+        dsignumy = -sum(dsigj.*dexpn.*(dy./sig0^2),2);
+        dsig = dsignum./dsigden;
+        dsigx = dsignumx./dsigden;
+        dsigy = dsignumy./dsigden;
+
+        val = dsig;
+
+        sz = size(dsig);
+        grad = zeros([sz,2]);
+        grad(:,:,1) = dsigx;
+        grad(:,:,2) = dsigy;
+
+        grad = squeeze(grad);
+
+end

chunkie/+chnk/+smoother/get_umesh.m

@@ -0,0 +1,23 @@
+function [umesh] = get_umesh(verts)
+    umesh = [];
+    umesh.verts = verts;
+    [~, nv] = size(verts);
+    verts_ext = [verts, verts(:,1)];
+    d = verts_ext(:,2:end) - verts_ext(:,1:nv);
+    umesh.centroids = 0.5*(verts_ext(:,2:end) + verts_ext(:,1:nv));
+    dd = sqrt(sum(d.^2, 1));
+    umesh.lengths = dd;
+    face_normals = zeros(2,nv);
+    face_normals(1,:) = d(2,:)./dd;
+    face_normals(2,:) = -d(1,:)./dd;
+
+    umesh.face_normals = face_normals;
+    face_normals_ext = [face_normals(:,end), face_normals];
+    pseudo_normals = 0.5*(face_normals_ext(:,2:end) + ...
+               face_normals_ext(:,1:nv));
+    dd = sqrt(sum(pseudo_normals.^2, 1));
+    pseudo_normals = pseudo_normals./dd;
+
+    umesh.pseudo_normals = pseudo_normals;
+
+end

chunkie/+chnk/+smoother/gpsi_all.m

@@ -0,0 +1,65 @@
+function [val,grad,hess,hess_sig] = gpsi_all(dx,dy,dr,dsigt)
+
+    [val,vald] = chnk.smoother.psi_eval(dr,dsigt);
+
+    da = 1/(4*pi);
+    z  = dr.^2./(2*dsigt.*dsigt);
+    val = val + log(z)*da;
+
+    if (nargout==1) 
+        return
+    end
+
+    grad = zeros(numel(dx),2);
+    grad(:,1) = vald.*dx./dr;
+    grad(:,2) = vald.*dy./dr;
+
+    grad(:,1) = grad(:,1)+2*dx./(dr.*dr)*da;
+    grad(:,2) = grad(:,2)+2*dy./(dr.*dr)*da;
+
+    if (nargout==2)
+        return
+    end
+
+    hess = zeros(numel(dx),2,2);
+
+    hess(:,1,1) = 0;
+    hess(:,1,2) = 0;
+    hess(:,2,1) = 0;
+    hess(:,2,2) = 0;
+
+    eterm = exp(-dr.*dr./(2*dsigt.^2));
+
+    hess(:,1,1) = eterm./(dr.*dr);
+    hess(:,2,2) = eterm./(dr.*dr);
+
+    fact = -(2*dsigt.^2+dr.^2)./(dsigt.^2.*dr.^2).*eterm;
+    hess(:,1,1) = hess(:,1,1) + (dx.*dx./dr.^2).*fact;
+    hess(:,1,2) = hess(:,1,2) + (dx.*dy./dr.^2).*fact;
+    hess(:,2,1) = hess(:,2,1) + (dx.*dy./dr.^2).*fact;
+    hess(:,2,2) = hess(:,2,2) + (dy.*dy./dr.^2).*fact;
+
+    hess(:,1,1) = hess(:,1,1)/(2*pi);
+    hess(:,2,1) = hess(:,2,1)/(2*pi);
+    hess(:,1,2) = hess(:,1,2)/(2*pi);
+    hess(:,2,2) = hess(:,2,2)/(2*pi);
+
+    hess(:,1,1) = hess(:,1,1) - 2./(dr.*dr)*da;
+    hess(:,2,2) = hess(:,2,2) - 2./(dr.*dr)*da;
+
+    fact = da*4./dr.^4;
+    hess(:,1,1) = hess(:,1,1) + dx.*dx.*fact;
+    hess(:,1,2) = hess(:,1,2) + dx.*dy.*fact;
+    hess(:,2,1) = hess(:,2,1) + dx.*dy.*fact;
+    hess(:,2,2) = hess(:,2,2) + dy.*dy.*fact;
+
+    if (nargout==3)
+        return
+    end
+
+
+    fact = eterm./(dsigt.^3)/(2*pi);
+    hess_sig(:,1) = -dx.*fact;
+    hess_sig(:,2) = -dy.*fact;
+
+end

chunkie/+chnk/+smoother/gpsi_loc.m

@@ -0,0 +1,48 @@
+function [val,grad,hess,hess_sig] = gpsi_loc(dx,dy,dr,dsigt)
+
+    degam = 0.57721566490153286060651209008240243104;
+    da = sqrt(2)*dsigt;
+    z = dr./da;
+    z2 = z.*z;
+    val = (z2-z2.^2/4+z2.^3/18-z2.^4/96+z2.^5/600-z2.^6/4320);
+    val = val./(4*pi);
+    val = val - degam/(4*pi); 
+
+    if (nargout==1) 
+        return
+    end
+
+    grad = zeros(numel(dx),2);
+    vvd = (1-z2/2+z2.^2/6-z2.^3/24+z2.^4/120-z2.^5/720);
+    grad(:,1) = dx.*vvd./(2*pi*da.^2);
+    grad(:,2) = dy.*vvd./(2*pi*da.^2);
+
+    if (nargout==2) 
+        return
+    end
+
+    hess = zeros(numel(dx),2,2);
+    hess(:,1,1) = -vvd./(2*pi*da.^2);
+    hess(:,2,2) = -vvd./(2*pi*da.^2);
+    hess(:,1,2) = 0;
+    hess(:,2,1) = 0;
+
+    vvvd = (-1+4*z2/6-6*z2.^2/24+8*z2.^3/120-10*z2.^4/720);
+    vvvd = -vvvd./(2*pi*da.^2)./da./da;
+    hess(:,1,1) = hess(:,1,1) + dx.*dx.*vvvd;
+    hess(:,2,2) = hess(:,2,2) + dy.*dy.*vvvd;
+    hess(:,1,2) = hess(:,1,2) + dx.*dy.*vvvd;
+    hess(:,2,1) = hess(:,2,1) + dx.*dy.*vvvd;
+
+    if (nargout==3)
+        return
+    end
+
+    hess_sig = zeros(numel(dx),2);
+
+    eterm = exp(-dr.*dr./(2*dsigt.^2));
+    fact = eterm./(dsigt.^3)/(2*pi);
+    hess_sig(:,1) = -dx.*fact;
+    hess_sig(:,2) = -dy.*fact;
+
+end

chunkie/+chnk/+smoother/green.m

@@ -0,0 +1,62 @@
+function [val,grad,hess,hess_sig] = green(rs,rt,dsigt)
+
+    [~,ns] = size(rs);
+    [~,nt] = size(rt);
+
+    xs = repmat(rs(1,:),nt,1);
+    ys = repmat(rs(2,:),nt,1);
+
+    xt = repmat(rt(1,:).',1,ns);
+    yt = repmat(rt(2,:).',1,ns);
+    dst = repmat(dsigt.',1,ns);
+
+    dx = xt-xs;
+    dy = yt-ys;
+    dr = sqrt(dx.^2+dy.^2);
+
+    sz = size(dx);
+    dx = dx(:);
+    dy = dy(:);
+    dr = dr(:);
+    dst = dst(:);
+
+    da = sqrt(2.0d0)*dst;
+    z = dr./da;
+
+    iloc = find(z < 0.1);
+    ifar = find(z >= 0.1);
+
+    vs = zeros(ns*nt,1);
+    gs = zeros(ns*nt,2);
+    hs = zeros(ns*nt,2,2);
+    hsigs = zeros(ns*nt,2);
+
+    if (numel(iloc) ~=0)
+        [vl,gl,hl,hsl] = chnk.smoother.gpsi_loc(dx(iloc), dy(iloc), ...
+                dr(iloc), dst(iloc));
+        vs(iloc) = vl;
+        gs(iloc,:) = gl;
+        hs(iloc,:,:) = hl;
+        hsigs(iloc,:) = hsl;
+    end
+    if (numel(ifar) ~=0)
+        [vf,gf,hf,hsf] = chnk.smoother.gpsi_all(dx(ifar), dy(ifar), ...
+              dr(ifar), dst(ifar));
+        vs(ifar) = vf;
+        gs(ifar,:) = gf;
+        hs(ifar,:,:) = hf;
+        hsigs(ifar,:) = hsf;
+    end
+
+    val = reshape(vs,sz);
+    if (nargout > 1)
+        grad = squeeze(reshape(gs,[sz,2]));
+    end
+    if (nargout > 2)
+        hess = squeeze(reshape(hs,[sz,2,2]));
+    end
+    if (nargout > 3)
+        hess_sig = squeeze(reshape(hsigs,[sz,2]));
+    end
+
+end