prior_gaussian_HPF

 PRIOR_GAUSSIAN_HPF calculate image prior
 Reg= prior_gaussian_HPF( fwd_model )
     or accepts inv_model
 Reg        => output regularization term
 fwd_model  => forward model struct
 Parameters:
   diam_frac= fwd_model.prior_gaussian_HPF.diam_frac DEFAULT 0.1
   zero_thresh=fwd_model.prior_gaussian_HPF.zero_thresh DEFAULT 0.0001


 prior_gaussian_HPF is designed to be used as an R_prior, rather than a RtR_prior

 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( fwd_model );
0002 % PRIOR_GAUSSIAN_HPF calculate image prior
0003 % Reg= prior_gaussian_HPF( fwd_model )
0004 %     or accepts inv_model
0005 % Reg        => output regularization term
0006 % fwd_model  => forward model struct
0007 % Parameters:
0008 %   diam_frac= fwd_model.prior_gaussian_HPF.diam_frac DEFAULT 0.1
0009 %   zero_thresh=fwd_model.prior_gaussian_HPF.zero_thresh DEFAULT 0.0001
0010 %
0011 %
0012 % prior_gaussian_HPF is designed to be used as an R_prior, rather than a RtR_prior
0013 %
0014 % CITATION_REQUEST:
0015 % AUTHOR: A Adler & R Guardo
0016 % YEAR: 1996
0017 % TITLE: Electrical impedance tomography: regularized imaging and contrast
0018 % detection
0019 % JOURNAL: IEEE transactions on medical imaging
0020 % VOL: 15
0021 % NUM: 2
0022 % PAGE: 170–9
0023 % LINK: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=491418
0024 
0025 % NOTES:
0026 % - the old version had numerous bugs. It has been fixed and simplified for version 3.11
0027 % - Old code should be removed in future
0028 % TODO:
0029 % - Accept and work with coarse2fine
0030 
0031 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0032 % $Id: prior_gaussian_HPF.m 6510 2022-12-30 16:34:22Z aadler $
0033 
0034 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0035 
0036 if ~strcmp(fwd_model.type, 'fwd_model')
0037     fwd_model= fwd_model.fwd_model;
0038 end
0039 try 
0040     diam_frac= fwd_model.prior_gaussian_HPF.diam_frac;
0041 catch
0042     diam_frac= 0.1;
0043 end
0044 try 
0045     diam_frac= fwd_model.prior_gaussian_HPF.zero_thresh;
0046 catch
0047     zero_thresh = .0001;
0048 end
0049 
0050 copt.cache_obj= {fwd_model.nodes, fwd_model.elems, diam_frac};
0051 copt.fstr = 'prior_gaussian_HPF';
0052 Reg = eidors_cache(@calc_Gaussian_HPF, {fwd_model, diam_frac, zero_thresh}, copt );
0053 
0054 function filt= calc_Gaussian_HPF( fmdl, diam_frac, zero_thresh)
0055   pts = interp_mesh(fmdl,3); % elem ctr
0056   fmdl.interp_mesh.n_interp = 0; % for ctr
0057   ptc = interp_mesh(fmdl,0);
0058   % we divide by four to closely match
0059   % previous implementation
0060   mdl_dim = max(fmdl.nodes) - min(fmdl.nodes);
0061   mdl_dim = mean(mdl_dim(1:2)); 
0062    
0063   beta=2.769/(diam_frac * mdl_dim)^2;
0064   dim = elem_dim(fmdl);
0065 
0066   % Gaussian is integral of exp(-beta*r^2)
0067   % Integral of exp(-1/2 * x'*Sigma*x)
0068   %  = sqrt((2*pi)^k*det(Sigma))
0069   % Here Sigma = 2*eye()*beta
0070   %  det(Sigma) = (2*beta)^k
0071   % so int[exp(-beta*r2)]= (pi*beta)^(d/2)
0072 
0073   % integral of exp(-r2) =
0074   Abeta_pi = get_elem_volume(fmdl,'no_c2f') ...
0075              *(beta/pi)^(dim/2);
0076   for j= 1:num_elems(fmdl)
0077     dif = bsxfun(@minus, pts, ptc(j,:));
0078     r2  = sum(dif.^2, 2);
0079     mean_elem = mean( exp(-beta*r2 ),3);
0080     flt(:,j)=Abeta_pi.*mean_elem;
0081   end %for j=1:ELEM
0082   % flt should sum to 1, but the integral is numeric and could vary. Instead we normalize
0083   flt = flt./sum(flt);
0084   % Filter is 1 - Gaussian
0085   flt=eye(num_elems(fmdl)) - flt;
0086   flt= ( flt+flt' )/ 2;
0087   filt= sparse(flt.*(abs(flt)>zero_thresh)); 
0088 
0089 % Calculate Gaussian HP Filter as per Adler & Guardo 96
0090 % parameter is diam_frac (normally 0.1)
0091 function filt= calc_Gaussian_HPF_old( fmdl, diam_frac)
0092   ELEM= fmdl.elems';
0093   NODE= fmdl.nodes';
0094 
0095 
0096   e= size(ELEM, 2);
0097   np= 128;
0098   [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0099 
0100   v_yx= [-y,x];
0101   o= ones(np*np,1);
0102   filt= zeros(e);
0103   tourne= [0 -1 1;1 0 -1;-1 1 0];
0104 
0105   for j= 1:e
0106 %   if ~rem(j,20); fprintf('.'); end
0107     xy= NODE(:,ELEM(:,j))';
0108     a= xy([2;3;1],1).*xy([3;1;2],2) ...
0109          -xy([3;1;2],1).*xy([2;3;1],2);
0110     aire=abs(sum(a));
0111     mx_xy = max(xy);
0112     mn_xy = min(xy);
0113     endr=find(y<=mx_xy(2) & y>=mn_xy(2) ...
0114             & x<=mx_xy(1) & x>=mn_xy(1) )';
0115     aa= sum(abs(ones(length(endr),1)*a' ...
0116          +v_yx(endr,:)*xy'*tourne)');
0117     endr( find( abs(1 - aa / aire) > 1e-8 ) )=[];
0118     ll=length(endr); endr=endr-1;
0119 
0120     yp= rem(endr,np)/(np-1) - .5; % (rem(endr,np) corresponde a y
0121     ym= ones(e,1)*yp -yc*ones(1,ll);
0122     xp= floor(endr/np)/(np-1) - .5; % (floor(endr/np)) corresponde a x
0123     xm= ones(e,1)*xp -xc*ones(1,ll);
0124 
0125     beta=2.769/diam_frac.^2;
0126 %   filt(:,j)=-aire/2*beta/pi*mean(...
0127     filt(:,j)=-beta/pi*sum( exp(-beta*(ym.^2+xm.^2))')';
0128   end %for j=1:ELEM
0129 % filt=filt/taille(1)/taille(2)+eye(e);
0130   filt=filt/np^2+eye(e);
0131   filt= ( filt+filt' )/ 2;
0132   filt= sparse(filt.*(abs(filt)>.003)); 
0133 
0134 function [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0135   taille=max(NODE')-min(NODE');
0136   e= size(ELEM, 2);
0137 
0138 % Triangles of each shape
0139   xt= reshape(NODE(1,ELEM(:)),3,e)';
0140   yt= reshape(NODE(2,ELEM(:)),3,e)';
0141 
0142 % We want center [1,1,1]/3 and edges [4,1,1]/6
0143   pts= [2,2,2;4,1,1;1,4,1;1,1,4]'/6;
0144   xp= xt*pts;
0145   yp= yt*pts;
0146   
0147   [x y]=meshgrid( ...
0148       linspace( min(NODE(1,:)), max(NODE(1,:)) ,np ), ...
0149       linspace( min(NODE(2,:)), max(NODE(2,:)) ,np )  ); 
0150 % Add the basic interpolation points to those based on the
0151 %  elements
0152   x= [x(:);xp(:)]; 
0153   y= [y(:);yp(:)]; 
0154 
0155 % BUG HERE: why is xc,yc scaled and not x,y
0156 
0157   xc= mean(xt,2)/taille(1);
0158   yc= mean(yt,2)/taille(2);
0159 
0160 function do_unit_test
0161   imdl = mk_common_model('a2c0',16);
0162   RtR = prior_gaussian_HPF(imdl);
0163   tt=[0.557367798877852, -0.124016055611277, -0.029231532288002, -0.124016055611277];
0164   unit_test_cmp('a2c2 :1', RtR(1,1:4),tt,1e-10);
0165 
0166   % Old values are somewhat similar
0167   RtR= calc_Gaussian_HPF_old( imdl.fwd_model, 0.1);
0168   tt=[0.562239752317943, -0.117068756722254, -0.025875127622824, -0.117068756722254];
0169   unit_test_cmp('a2c2 _old :1', RtR(1,1:4),tt,1e-10);
0170 
0171   imdl = mk_common_model('a3cr',16);
0172   RtR = prior_gaussian_HPF(imdl);  %NOTE: Fix required
0173   tt=[0.793660877689329,-0.096517734290331;
0174      -0.096517734290331, 0.793660877689329;
0175      -0.004236643618993,-0.011516904603663;
0176      -0.011516904603663,-0.004236643618993];
0177   unit_test_cmp('a3cr :1', RtR(1:4,1:2),tt,1e-10);
0178 
0179   ls = -1:.1:1;
0180   fmdl = mk_grid_model([], ls, ls); 
0181   fmdl = rmfield(fmdl,'coarse2fine');
0182   RtR = prior_gaussian_HPF(fmdl);
0183   
0184   % Test space-independence ... 2D
0185   ll= 2*(length(ls)-1);
0186   unit_test_cmp('2D space indep', ...
0187      RtR(10*ll+(1:4),10*ll+(1:4)), ...
0188      RtR(12*ll+(1:4),12*ll+(1:4)), 1e-14);
0189 
0190 
0191   ls = -1:.2:1;
0192   fmdl = mk_grid_model([],ls,ls,(-4:2:4)/10); 
0193   fmdl = rmfield(fmdl,'coarse2fine');
0194   RtR = prior_gaussian_HPF(fmdl);
0195   
0196   % Test space-independence ... 2D
0197   ll= 6*(length(ls)-1);
0198   lv= 6*(length(ls)-1)^2;
0199   unit_test_cmp('3D space indep', ...
0200      RtR(lv+3*ll+(1:4),lv+3*ll+(1:4)), ...
0201      RtR(lv+5*ll+(1:4),lv+5*ll+(1:4)), 1e-14);
0202