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 %simp = The simplices matrix.
0007 %vtx  = The vertices matrix.
0008 %deg  = 1 for nodes, 2 for edges and 3 for faces
0009 %w    = smoothing weight, w=1...k, default value = 1
0010 %Reg  = The first order smoothing regulariser.
0011 
0012 if ischar(simp) && strcmp(simp,'UNIT_TEST'); do_unit_test; return; end
0013 
0014 warning('EIDORS:deprecated','ISO_F_SMOOTH is deprecated as of 06-Jun-2012. ');
0015 % It is recommended to use the newer prior functions in EIDORS. These
0016 % produce well documented prior matrices.
0017 
0018 
0019 
0020 if nargin<2
0021    w=1;
0022 end
0023 if w<0
0024    error('Weight must be possitive');
0025 end
0026 
0027 ndsrch = [simp,(1:size(simp,1))'];
0028 
0029 Reg = spalloc(size(simp,1),size(simp,1),20*size(simp,1));
0030 
0031 for i=1:size(ndsrch,1)
0032 
0033    t_id = ndsrch(i,1:size(simp,2));
0034 
0035    Xsimp = []; %The vector of simp indices that share deg nodes with simp(i)
0036    XX = [];
0037 
0038    for j = 1:length(t_id)
0039 
0040       t_nd = t_id(j);
0041 
0042       X1 = find(simp(:,1)==t_nd);
0043       X2 = find(simp(:,2)==t_nd);
0044       X3 = find(simp(:,3)==t_nd);
0045       X4 = find(simp(:,4)==t_nd);
0046 
0047       %The set of all indices containing the node t_nd
0048       XX = [XX;X1;X2;X3;X4];
0049 
0050    end %for j
0051 
0052    if deg == 1 %Neigbouring nodes
0053       Xu = unique(XX);
0054       Xsimp = [Xsimp;Xu];
0055    end
0056 
0057    if deg == 2 %Neihgouring edges
0058       Xs = sort(XX);
0059       Xd = diff(Xs);
0060       Xf = Xs(setdiff(1:length(Xd),find(Xd)));
0061       Xu = unique(Xf);
0062       Xsimp = [Xsimp;Xu];
0063    end
0064 
0065    if deg == 3 %Neihbouring faces
0066       Xs = sort(XX);
0067       Xd = diff(Xs);
0068       Xf = Xs(setdiff(1:length(Xd),find(Xd)));
0069       Xq = diff(Xf);
0070       Xz = Xf(setdiff(1:length(Xq),find(Xq)));
0071       Xu = unique(Xz);
0072       Xsimp = [Xsimp;Xu];
0073    end
0074 
0075    if deg > 3 || deg < 1
0076       error('deg parameter can only be 1, 2 or 3')
0077    end
0078 
0079    %Remove the i'th simplex from the list of its neighbours
0080    Xdif = setdiff(Xsimp,i);
0081 
0082    if deg == 1
0083       Reg(i,Xdif) = -1/w;
0084    end
0085 
0086    if deg == 2
0087       for p=1:length(Xdif)
0088          Intr = intersect(t_id,simp(Xdif(p),:)); % 2 or 3 long vector
0089          dd=0;
0090          for h=1:length(Intr)-1
0091             [da] = db23d(vtx(Intr(h),1),vtx(Intr(h),2),vtx(Intr(h),3),...
0092             vtx(Intr(h+1),1),vtx(Intr(h+1),2),vtx(Intr(h+1),3));
0093             dd = dd+da;
0094          end
0095          Reg(i,Xdif(p)) = -w/dd;
0096       end
0097    end
0098 
0099    if deg == 3
0100       for p=1:length(Xdif)
0101          Intr = intersect(t_id,simp(Xdif(p),:)); % 3 long vector
0102          [ta] = triarea3d(vtx(Intr,:));
0103       end
0104       Reg(i,Xdif) = -w/ta;
0105    end
0106 
0107    Reg(i,i) = abs(sum(Reg(i,:)));
0108 
0109 end %for i'th simplex
0110 
0111 
0112 function [dd] = db23d(x1,y1,z1,x2,y2,z2);
0113 %Auxiliary function that caclulates the distance between
0114 %two points or two sets of points in 3D
0115 %
0116 %(x1,y1,z1) = The coordinates of the first point(s) in 3D
0117 %(x2,y2,z2) = The coordinates of the second point(s) in 3D
0118 %dd         = Their distance(s)
0119 
0120 dd = sqrt((x2 - x1).^2 + (y2 - y1).^2 + (z2 - z1).^2);
0121 
0122 
0123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0124 % This is part of the EIDORS suite.
0125 % Copyright (c) N. Polydorides 2003
0126 % Copying permitted under terms of GNU GPL
0127 % See enclosed file gpl.html for details.
0128 % EIDORS 3D version 2.0
0129 % MATLAB version 5.3 R11
0130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0131 
0132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0133 % This is part of the EIDORS suite.
0134 % Copyright (c) N. Polydorides 2003
0135 % Copying permitted under terms of GNU GPL
0136 % See enclosed file gpl.html for details.
0137 % EIDORS 3D version 2.0
0138 % MATLAB version 5.3 R11
0139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0140 
0141 function do_unit_test
0142    fmdl = mk_circ_tank(1,[0,.5],4);
0143    simp = fmdl.elems;
0144    vtx  = fmdl.nodes;
0145 subplot(221)
0146    show_fem(fmdl,[1,1,1]);
0147 
0148    Ref = iso_f_smooth(simp,vtx,1,1);
0149    Ref = iso_f_smooth(simp,vtx,2,1);
0150    Ref = iso_f_smooth(simp,vtx,3,1);
0151 subplot(222); spy(Ref)
0152    Rel = prior_laplace(fmdl);
0153 subplot(224); spy(Rel);
0154 
0155 %% The shapes should be the same but the
0156 %% Matrices are not symmetric from iso_f_smooth
0157