laplace_image_prior

 LAPLACE_IMAGE_PRIOR calculate image prior
 Reg= laplace_image_prior( inv_model )
 Reg        => output regularization term
 inv_model  => inverse model struct

 This image prior is intended to be used as
  R'*R, but may be used as R for as well.

 The Laplace prior is a 2nd order high pass filter.
 On a rectangular mesh, it is a convolution with
   [-1,-1,-1;      [ 0;-1; 0
    -1, 8,-1    or  -1; 4;-1
    -1,-1,-1]        0;-1; 0]

 On a finite element mesh, we define the it as 
 -1 for each adjacent element, and 3 (in 2D) or 4 (in 3D)
 for the element itself

0001 function Reg= laplace_image_prior( inv_model );
0002 % LAPLACE_IMAGE_PRIOR calculate image prior
0003 % Reg= laplace_image_prior( inv_model )
0004 % Reg        => output regularization term
0005 % inv_model  => inverse model struct
0006 %
0007 % This image prior is intended to be used as
0008 %  R'*R, but may be used as R for as well.
0009 %
0010 % The Laplace prior is a 2nd order high pass filter.
0011 % On a rectangular mesh, it is a convolution with
0012 %   [-1,-1,-1;      [ 0;-1; 0
0013 %    -1, 8,-1    or  -1; 4;-1
0014 %    -1,-1,-1]        0;-1; 0]
0015 %
0016 % On a finite element mesh, we define the it as
0017 % -1 for each adjacent element, and 3 (in 2D) or 4 (in 3D)
0018 % for the element itself
0019 
0020 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0021 % $Id: laplace_image_prior.html 2819 2011-09-07 16:43:11Z aadler $
0022 
0023 pp= aa_fwd_parameters( inv_model.fwd_model );
0024 
0025 Reg = speye( pp.n_elem );
0026 
0027    Iidx= [];
0028    Jidx= [];
0029    Vidx= [];
0030    for ii=1:pp.n_elem
0031      el_adj = find_adjoin( ii, pp.ELEM );
0032      for jj=el_adj(:)'
0033          Iidx= [Iidx, ii, ii, jj, jj];
0034          Jidx= [Jidx, ii, jj, ii, jj];
0035          Vidx= [Vidx,  1, -1, -1,  1];
0036      end
0037    end
0038 
0039    Reg = sparse(Iidx,Jidx, Vidx, pp.n_elem, pp.n_elem );
0040 
0041 % find elems which are connected to elems ee
0042 function elems= find_adjoin(ee, ELEM)
0043    nn= ELEM(:,ee);
0044    [d,e]= size(ELEM);
0045    ss= zeros(1,e);
0046    for i=1:d
0047      ss= ss+ any(ELEM==nn(i));
0048    end
0049    elems= find(ss==d-1);