gaussian_HPF_prior

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

0001 function Reg= gaussian_HPF_prior( inv_model );
0002 % GAUSSIAN_HPF_PRIOR calculate image prior
0003 % Reg= gaussian_HPF_prior( inv_model )
0004 % Reg        => output regularization term
0005 % inv_model  => inverse model struct
0006 % Parameters:
0007 %   diam_frac= inv_model.fwd_model.gaussian_HPF_prior.diam_frac DEFAULT 0.1
0008 
0009 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0010 % $Id: gaussian_HPF_prior.html 2819 2011-09-07 16:43:11Z aadler $
0011 
0012 fwd_model= inv_model.fwd_model;
0013 try 
0014     diam_frac= fwd_model.gaussian_HPF_prior.diam_frac;
0015 catch
0016     diam_frac= 0.1;
0017 end
0018 
0019 cache_test_obj= {fwd_model.nodes, fwd_model.elems, diam_frac};
0020 Reg = eidors_obj('get-cache', cache_test_obj, 'gaussian_HPF_prior');
0021 if ~isempty(Reg)
0022    eidors_msg('gaussian_HPF_prior: using cached value', 3);
0023    return
0024 end
0025 
0026 Reg = calc_Gaussian_HPF( fwd_model, diam_frac );
0027 
0028 cache_test_obj= {fwd_model.nodes, fwd_model.elems, diam_frac};
0029 eidors_obj('set-cache', cache_test_obj, 'gaussian_HPF_prior', Reg);
0030 eidors_msg('gaussian_HPF_prior: setting cached value', 3);
0031 
0032 % Calculate Gaussian HP Filter as per Adler & Guardo 96
0033 % parameter is diam_frac (normally 0.1)
0034 function filt= calc_Gaussian_HPF( fmdl, diam_frac)
0035   ELEM= fmdl.elems';
0036   NODE= fmdl.nodes';
0037 
0038 
0039   e= size(ELEM, 2);
0040   np= 128;
0041   [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0042 
0043   v_yx= [-y,x];
0044   o= ones(np*np,1);
0045   filt= zeros(e);
0046   tourne= [0 -1 1;1 0 -1;-1 1 0];
0047 
0048   for j= 1:e
0049 %   if ~rem(j,20); fprintf('.'); end
0050     xy= NODE(:,ELEM(:,j))';
0051     a= xy([2;3;1],1).*xy([3;1;2],2) ...
0052          -xy([3;1;2],1).*xy([2;3;1],2);
0053     aire=abs(sum(a));
0054     endr=find(y<=max(xy(:,2)) & y>=min(xy(:,2)) ...
0055             & x<=max(xy(:,1)) & x>=min(xy(:,1)) )';
0056     aa= sum(abs(ones(length(endr),1)*a' ...
0057          +v_yx(endr,:)*xy'*tourne)');
0058     endr( find( abs(1 - aa / aire) > 1e-8 ) )=[];
0059     ll=length(endr); endr=endr-1;
0060 
0061     yp= rem(endr,np)/(np-1) - .5; % (rem(endr,np) corresponde a y
0062     ym= ones(e,1)*yp -yc*ones(1,ll);
0063     xp= floor(endr/np)/(np-1) - .5; % (floor(endr/np)) corresponde a x
0064     xm= ones(e,1)*xp -xc*ones(1,ll);
0065 
0066     beta=2.769/diam_frac.^2;
0067 %   filt(:,j)=-aire/2*beta/pi*mean(...
0068     filt(:,j)=-beta/pi*sum( exp(-beta*(ym.^2+xm.^2))')';
0069   end %for j=1:ELEM
0070 % filt=filt/taille(1)/taille(2)+eye(e);
0071   filt=filt/np^2+eye(e);
0072   filt= ( filt+filt' )/ 2;
0073   filt= sparse(filt.*(abs(filt)>.003)); 
0074 
0075 function [x,xc,y,yc] = interp_points(NODE,ELEM,np);
0076   taille=max(NODE')-min(NODE');
0077   e= size(ELEM, 2);
0078 
0079 % Triangles of each shape
0080   xt= reshape(NODE(1,ELEM(:)),3,e)';
0081   yt= reshape(NODE(2,ELEM(:)),3,e)';
0082 
0083 % We want center [1,1,1]/3 and edges [4,1,1]/6
0084   pts= [2,2,2;4,1,1;1,4,1;1,1,4]'/6;
0085   xp= xt*pts;
0086   yp= yt*pts;
0087   
0088   [x y]=meshgrid( ...
0089       linspace( min(NODE(1,:)), max(NODE(1,:)) ,np ), ...
0090       linspace( min(NODE(2,:)), max(NODE(2,:)) ,np )  ); 
0091 % Add the basic interpolation points to those based on the
0092 %  elements
0093   x= [x(:);xp(:)]; 
0094   y= [y(:);yp(:)]; 
0095 
0096   xc= mean(xt,2)/taille(1);
0097   yc= mean(yt,2)/taille(2);