updating smoother

Author: mrachh

chunkie/+chnk/+smoother/get_sigs.m

@@ -1,7 +1,7 @@
-function [val,grad] = get_sigs(src_input,rt,sig0,dlam)
+function [val,grad] = get_sigs(umesh,rt,sig0,dlam)
 
-        rc = src_input.centroids;
-        dlens = src_input.lengths;
+        rc = umesh.centroids;
+        dlens = umesh.lengths;
 
         [~,ns] = size(rc);
         [~,nt] = size(rt);

chunkie/+chnk/+smoother/gpsi_all.m

@@ -33,7 +33,7 @@
     hess(:,1,1) = eterm./(dr.*dr);
     hess(:,2,2) = eterm./(dr.*dr);
 
-    fact = -(2*dsigt.^2+dr.^2)./(dsigt.^2*dr.^2).*eterm;
+    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;
@@ -59,7 +59,7 @@
 
 
     fact = eterm./(dsigt.^3)/(2*pi);
-    hess_sig(:,1) = -dx*fact;
-    hess_sig(:,2) = -dy*fact;
+    hess_sig(:,1) = -dx.*fact;
+    hess_sig(:,2) = -dy.*fact;
 
 end

chunkie/+chnk/+smoother/green.m

@@ -1,61 +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)*dsigt;
-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,dy,dr,dst);
-    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,dy,dr,dsigt);
-    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
+    
+    [~,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

chunkie/+chnk/+smoother/smooth_curve.m

@@ -1,7 +1,7 @@
 nv = 3;
 z = exp(1j*2*pi*(1:nv)/nv);
 verts = [real(z); imag(z)];
-nchs = randi([1,10], 1, nv);
+nchs = 3*ones(nv,1);
 k = 16;
 
 kquad = 32;
@@ -10,7 +10,7 @@
 
 umesh = chnk.smoother.get_umesh(verts);
 dmesh = chnk.smoother.get_mesh(umesh, nchs, k);
-qmesh = get_mesh(umesh, nchs, kquad);
+qmesh = chnk.smoother.get_mesh(umesh, nchs, kquad);
 
 figure(1)
 clf
@@ -23,5 +23,62 @@
 clf
 plot(dmesh.r(1,:), dmesh.r(2,:),'k.'); hold on;
 axis equal;
-% quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.n(1,:), dmesh.n(2,:),'k');
-quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
+quiver(dmesh.r(1,:), dmesh.r(2,:), dmesh.pseudo_normals(1,:), dmesh.pseudo_normals(2,:),'r');
+
+
+[~, nd] = size(dmesh.r);
+h = zeros(nd,1);
+n_newton = 5;
+
+sig0 = sqrt(5)*max(umesh.lengths);
+lam = 10;
+
+
+ww = qmesh.wts(:).';
+dpx = dmesh.pseudo_normals(1,:).';
+dpy = dmesh.pseudo_normals(2,:).';
+
+rnx = qmesh.n(1,:).';
+rny = qmesh.n(2,:).';
+
+for i=1:n_newton
+    rt = dmesh.r;
+    rt(1,:) = rt(1,:) + (h.*dpx).';
+    rt(2,:) = rt(2,:) + (h.*dpy).';
+
+    [sig, sig_grad] = chnk.smoother.get_sigs(umesh, rt, sig0, lam);
+    [val, grad, hess, hess_sig] = chnk.smoother.green(qmesh.r, rt, sig);
+    gx = grad(:,:,1); 
+    gy = grad(:,:,2);
+    
+    phi = -(gx*(rnx.*qmesh.wts) + gy*(rny.*qmesh.wts)) - 0.5;
+    
+    
+    h11 = hess(:,:,1,1).*ww;
+    h12 = hess(:,:,1,2).*ww;
+    h21 = hess(:,:,2,1).*ww;
+    h22 = hess(:,:,2,2).*ww;
+    
+    h1sig = hess_sig(:,:,1).*ww;
+    h2sig = hess_sig(:,:,2).*ww;
+
+    dx1 = h11*rnx + h12*rny;
+    dy1 = h21*rnx + h22*rny;
+    
+    dx2 = h1sig*rnx + h2sig*rny;
+    dy2 = dx2;
+    
+    dsigx = sig_grad(:,1);
+    dsigy = sig_grad(:,2);
+
+    dx2 = dx2.*dsigx;
+    dy2 = dy2.*dsigy;
+    
+    
+    dphidh = (dx1 + dx2).*dpx + (dy1 + dy2).*dpy;
+    h = h - phi./dphidh;
+
+    figure(i)
+    plot(rt(1,:), rt(2,:), 'k.')
+ 
+end