iso_f_smooth

function [Reg] = iso_f_smooth(simp,vtx,deg,w);

Calculates a first order discrete Gaussian smoothing operator.



simp = The simplices matrix.
vtx  = The vertices matrix.
deg  = 1 for nodes, 2 for edges and 3 for faces
w    = smoothing weight, w=1...k, default value = 1
Reg  = The first order smoothing regulariser.

0001 function [Reg] = iso_f_smooth(simp,vtx,deg,w);
0002 %function [Reg] = iso_f_smooth(simp,vtx,deg,w);
0003 %
0004 %Calculates a first order discrete Gaussian smoothing operator.
0005 %
0006 %
0007 %
0008 %simp = The simplices matrix.
0009 %vtx  = The vertices matrix.
0010 %deg  = 1 for nodes, 2 for edges and 3 for faces
0011 %w    = smoothing weight, w=1...k, default value = 1
0012 %Reg  = The first order smoothing regulariser.
0013 
0014 if nargin<2 
0015    w=1;
0016 end
0017 if w<0
0018    error('Weight must be possitive');
0019 end
0020 
0021 ndsrch = [simp,(1:size(simp,1))'];
0022 
0023 Reg = spalloc(size(simp,1),size(simp,1),20*size(simp,1));
0024 
0025 for i=1:size(ndsrch,1)
0026    
0027    t_id = ndsrch(i,1:size(simp,2));
0028    
0029    Xsimp = []; %The vector of simp indices that share deg nodes with simp(i)
0030    XX = [];
0031    
0032  for j = 1:length(t_id)
0033       
0034       t_nd = t_id(j);
0035       
0036       X1 = find(simp(:,1)==t_nd);
0037       X2 = find(simp(:,2)==t_nd);
0038       X3 = find(simp(:,3)==t_nd);
0039       X4 = find(simp(:,4)==t_nd);
0040       
0041       %The set of all indices containing the node t_nd
0042       XX = [XX;X1;X2;X3;X4];
0043       
0044   end %for j
0045       
0046       if deg == 1 %Neigbouring nodes
0047       Xu = unique(XX);
0048       Xsimp = [Xsimp;Xu];
0049       end
0050       
0051       if deg == 2 %Neihgouring edges
0052       Xs = sort(XX);
0053       Xd = diff(Xs);
0054       Xf = Xs(setdiff(1:length(Xd),find(Xd)));
0055       Xu = unique(Xf);
0056       Xsimp = [Xsimp;Xu];
0057       end
0058       
0059       if deg == 3 %Neihbouring faces
0060       Xs = sort(XX);
0061       Xd = diff(Xs);
0062       Xf = Xs(setdiff(1:length(Xd),find(Xd)));
0063       Xq = diff(Xf);
0064       Xz = Xf(setdiff(1:length(Xq),find(Xq)));
0065       Xu = unique(Xz);
0066       Xsimp = [Xsimp;Xu];
0067       end
0068       
0069       if deg > 3 || deg < 1 
0070           error('deg parameter can only be 1, 2 or 3')
0071       end
0072             
0073    %Remove the i'th simplex from the list of its neighbours
0074    Xdif = setdiff(Xsimp,i);
0075    
0076    if deg == 1
0077    Reg(i,Xdif) = -1/w;
0078    end
0079 
0080    if deg == 2
0081    for p=1:length(Xdif)
0082         Intr = intersect(t_id,simp(Xdif(p),:)); % 2 or 3 long vector
0083         dd=0;
0084        for h=1:length(Intr)-1
0085           [da] = db23d(vtx(Intr(h),1),vtx(Intr(h),2),vtx(Intr(h),3),...
0086                        vtx(Intr(h+1),1),vtx(Intr(h+1),2),vtx(Intr(h+1),3));
0087           dd = dd+da;
0088       end
0089    Reg(i,Xdif(p)) = -w/dd; 
0090    end
0091    end
0092 
0093    if deg == 3
0094     for p=1:length(Xdif)
0095        Intr = intersect(t_id,simp(Xdif(p),:)); % 3 long vector
0096        [ta] = triarea3d(vtx(Intr,:));
0097     end
0098    Reg(i,Xdif) = -w/ta;
0099    end
0100         
0101    Reg(i,i) = abs(sum(Reg(i,:)));
0102    
0103 end %for i'th simplex
0104 
0105 
0106 
0107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0108 % This is part of the EIDORS suite.
0109 % Copyright (c) N. Polydorides 2003
0110 % Copying permitted under terms of GNU GPL
0111 % See enclosed file gpl.html for details.
0112 % EIDORS 3D version 2.0
0113 % MATLAB version 5.3 R11
0114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0115