prior_laplace

 PRIOR_LAPLACE calculate image prior
 Reg= prior_laplace( inv_model )
 Reg        => output regularization term
 inv_model  => inverse model struct
  or
 Reg= prior_laplace( fwd_model )

 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. 

 The EIDORS implementation is equivalent to twice the Graph-Laplacian

0001 function Reg= prior_laplace( inv_model )
0002 % PRIOR_LAPLACE calculate image prior
0003 % Reg= prior_laplace( inv_model )
0004 % Reg        => output regularization term
0005 % inv_model  => inverse model struct
0006 %  or
0007 % Reg= prior_laplace( fwd_model )
0008 %
0009 % This image prior is intended to be used as
0010 %  R'*R, but may be used as R for as well.
0011 %
0012 % The Laplace prior is a 2nd order high pass filter.
0013 % On a rectangular mesh, it is a convolution with
0014 %   [-1,-1,-1;      [ 0;-1; 0
0015 %    -1, 8,-1    or  -1; 4;-1
0016 %    -1,-1,-1]        0;-1; 0]
0017 %
0018 % On a finite element mesh, we define the it as
0019 % -1 for each adjacent element, and 3 (in 2D) or 4 (in 3D)
0020 % for the element itself.
0021 %
0022 % The EIDORS implementation is equivalent to twice the Graph-Laplacian
0023 
0024 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: prior_laplace.m 5609 2017-06-23 18:05:53Z htregidgo $
0026 
0027 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 switch inv_model.type
0030     case 'inv_model'; fwd_model = inv_model.fwd_model;
0031     case 'fwd_model'; fwd_model = inv_model;
0032     otherwise; error('PRIOR_LAPLACE requires input type of inv_model or fwd_model');
0033 end
0034 
0035 copt.cache_obj = cache_obj(fwd_model);
0036 copt.fstr = 'prior_laplace';
0037 
0038 Reg = eidors_cache(@build_laplace, fwd_model, copt);
0039 
0040 % new implementation builds by running through the faces provided by
0041 % fix_model.
0042 function Reg = build_laplace(fwd_model)
0043 
0044 % obtain face2elem matrix
0045 fmopt.face2elem=true;
0046 fwd_model=fix_model(fwd_model,fmopt);
0047 
0048 n_elem=size(fwd_model.elems,1);
0049 
0050 % limit to interior faces
0051 interiorindx=fwd_model.face2elem(:,2)~=0;
0052 facei=fwd_model.face2elem(interiorindx,1);
0053 facej=fwd_model.face2elem(interiorindx,2);
0054 
0055 % build Idx as Iidx= [Iidx, ii, ii, jj, jj]; for each face
0056 Iidx=kron(facei,uint32([1;1;0;0]))+kron(facej,uint32([0;0;1;1]));
0057 
0058 % build Jidx as Jidx= [Jidx, ii, jj, ii, jj]; for each face
0059 Jidx=kron(facei,uint32([1;0;1;0]))+kron(facej,uint32([0;1;0;1]));
0060 
0061 % assign values to be consistent with previous implementation
0062 Vidx=kron(ones(size(facei)),[2;-2;-2;2]);
0063 
0064 % build
0065 Reg = sparse(Iidx,Jidx, Vidx, n_elem, n_elem );
0066 
0067 
0068 
0069 % Mapping depends only on elems - remove the other stuff
0070 function c_obj = cache_obj(f_mdl)
0071 c_obj = { f_mdl.elems, f_mdl.nodes};
0072 
0073 
0074 function do_unit_test
0075 
0076 imdl = mk_common_model('a2c2',16);
0077 RtR = prior_laplace( imdl );
0078 subplot(221); spy(RtR);
0079 unit_test_cmp('a2c2', nnz(RtR), 240);
0080 
0081 
0082 fmdl = mk_circ_tank(2,[],4);
0083 RtR = prior_laplace( fmdl );
0084 subplot(222); spy(RtR);
0085 unit_test_cmp('2-4',RtR(1:4,1:8), [ ...
0086     6, -2,  0, -2, -2,  0,  0,  0; -2,  6, -2,  0,  0, -2,  0,  0;
0087     0, -2,  6, -2,  0,  0, -2,  0; -2,  0, -2,  6,  0,  0,  0, -2]);
0088 
0089 compmdl = mk_common_model('n3r2',[16 2]);
0090 RtRnew = prior_laplace( compmdl);
0091 RtRold = prior_laplace_old(compmdl);
0092 subplot(223);spy(RtRnew);
0093 subplot(224);spy(RtRold);
0094 unit_test_cmp('consistent behaviour',RtRnew,RtRold);
0095