prior_gaussian_HPF

 PRIOR_GAUSSIAN_HPF calculate image prior
 Reg= prior_gaussian_HPF( inv_model )
 Reg        => output regularization term
 inv_model  => inverse model struct
 Parameters:
   diam_frac= inv_model.fwd_model.prior_gaussian_HPF.diam_frac DEFAULT 0.1

 CITATION_REQUEST:
 AUTHOR: A Adler & R Guardo
 YEAR: 1996
 TITLE: Electrical impedance tomography: regularized imaging and contrast
 detection 
 JOURNAL: IEEE transactions on medical imaging
 VOL: 15
 NUM: 2
 PAGE: 170–9
 LINK: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=491418

0001 function Reg= prior_gaussian_HPF( inv_model );
0002 % PRIOR_GAUSSIAN_HPF calculate image prior
0003 % Reg= prior_gaussian_HPF( inv_model )
0004 % Reg        => output regularization term
0005 % inv_model  => inverse model struct
0006 % Parameters:
0007 %   diam_frac= inv_model.fwd_model.prior_gaussian_HPF.diam_frac DEFAULT 0.1
0008 %
0009 % CITATION_REQUEST:
0010 % AUTHOR: A Adler & R Guardo
0011 % YEAR: 1996
0012 % TITLE: Electrical impedance tomography: regularized imaging and contrast
0013 % detection
0014 % JOURNAL: IEEE transactions on medical imaging
0015 % VOL: 15
0016 % NUM: 2
0017 % PAGE: 170–9
0018 % LINK: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=491418
0019 
0020 
0021 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0022 % $Id: prior_gaussian_HPF.m 5791 2018-05-22 15:16:02Z aadler $
0023 
0024 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0025 
0026 fwd_model= inv_model.fwd_model;
0027 try 
0028     diam_frac= fwd_model.prior_gaussian_HPF.diam_frac;
0029 catch
0030     diam_frac= 0.1;
0031 end
0032 
0033 copt.cache_obj= {fwd_model.nodes, fwd_model.elems, diam_frac};
0034 copt.fstr = 'prior_gaussian_HPF';
0035 if elem_dim(fwd_model) == 2;
0036    Reg = eidors_cache(@calc_Gaussian_HPF, {fwd_model, diam_frac}, copt );
0037 else
0038    warning('prior_gaussian_HPF: not yet able to generate 3D models');
0039    Reg = prior_laplace( inv_model );
0040 end
0041 
0042 
0043 % Calculate Gaussian HP Filter as per Adler & Guardo 96
0044 % parameter is diam_frac (normally 0.1)
0045 function filt= calc_Gaussian_HPF( fmdl, diam_frac)
0046   ELEM= fmdl.elems';
0047   NODE= fmdl.nodes';
0048 
0049 
0050   e= size(ELEM, 2);
0051   np= 128;
0052   [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0053 
0054   v_yx= [-y,x];
0055   o= ones(np*np,1);
0056   filt= zeros(e);
0057   tourne= [0 -1 1;1 0 -1;-1 1 0];
0058 
0059   for j= 1:e
0060 %   if ~rem(j,20); fprintf('.'); end
0061     xy= NODE(:,ELEM(:,j))';
0062     a= xy([2;3;1],1).*xy([3;1;2],2) ...
0063          -xy([3;1;2],1).*xy([2;3;1],2);
0064     aire=abs(sum(a));
0065     endr=find(y<=max(xy(:,2)) & y>=min(xy(:,2)) ...
0066             & x<=max(xy(:,1)) & x>=min(xy(:,1)) )';
0067     aa= sum(abs(ones(length(endr),1)*a' ...
0068          +v_yx(endr,:)*xy'*tourne)');
0069     endr( find( abs(1 - aa / aire) > 1e-8 ) )=[];
0070     ll=length(endr); endr=endr-1;
0071 
0072     yp= rem(endr,np)/(np-1) - .5; % (rem(endr,np) corresponde a y
0073     ym= ones(e,1)*yp -yc*ones(1,ll);
0074     xp= floor(endr/np)/(np-1) - .5; % (floor(endr/np)) corresponde a x
0075     xm= ones(e,1)*xp -xc*ones(1,ll);
0076 
0077     beta=2.769/diam_frac.^2;
0078 %   filt(:,j)=-aire/2*beta/pi*mean(...
0079     filt(:,j)=-beta/pi*sum( exp(-beta*(ym.^2+xm.^2))')';
0080   end %for j=1:ELEM
0081 % filt=filt/taille(1)/taille(2)+eye(e);
0082   filt=filt/np^2+eye(e);
0083   filt= ( filt+filt' )/ 2;
0084   filt= sparse(filt.*(abs(filt)>.003)); 
0085 
0086 function [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0087   taille=max(NODE')-min(NODE');
0088   e= size(ELEM, 2);
0089 
0090 % Triangles of each shape
0091   xt= reshape(NODE(1,ELEM(:)),3,e)';
0092   yt= reshape(NODE(2,ELEM(:)),3,e)';
0093 
0094 % We want center [1,1,1]/3 and edges [4,1,1]/6
0095   pts= [2,2,2;4,1,1;1,4,1;1,1,4]'/6;
0096   xp= xt*pts;
0097   yp= yt*pts;
0098   
0099   [x y]=meshgrid( ...
0100       linspace( min(NODE(1,:)), max(NODE(1,:)) ,np ), ...
0101       linspace( min(NODE(2,:)), max(NODE(2,:)) ,np )  ); 
0102 % Add the basic interpolation points to those based on the
0103 %  elements
0104   x= [x(:);xp(:)]; 
0105   y= [y(:);yp(:)]; 
0106 
0107   xc= mean(xt,2)/taille(1);
0108   yc= mean(yt,2)/taille(2);
0109 
0110 function do_unit_test
0111   imdl = mk_common_model('a2c0',16);
0112   RtR = prior_gaussian_HPF(imdl);
0113   tt=[0.562239752317943, -0.117068756722254, -0.025875127622824, -0.117068756722254;
0114      -0.117068756722254,  0.562239752317943, -0.117068756722254, -0.025875127622824;
0115      -0.025875127622824, -0.117068756722254,  0.562239752317943, -0.117068756722254;
0116      -0.117068756722254, -0.025875127622824, -0.117068756722254,  0.562239752317943];
0117   unit_test_cmp('a2c2 :1', RtR(1:4,1:4),tt,1e-10);
0118 
0119   imdl = mk_common_model('a3cr',16);
0120   RtR = prior_gaussian_HPF(imdl);  %NOTE: Fix required
0121   tt = [6    -2     0     0; -2     6     0     0;
0122         0     0     6    -2;  0     0    -2     6];
0123   unit_test_cmp('a3cr :1', RtR(1:4,1:4),tt,1e-10);
0124