prior_laplace

PURPOSE ^

PRIOR_LAPLACE calculate image prior

SYNOPSIS ^

function Reg= prior_laplace( inv_model );

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: prior_laplace.m 5112 2015-06-14 13:00:41Z aadler $
0024 
0025 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 switch inv_model.type
0028   case 'inv_model'; fwd_model = inv_model.fwd_model;
0029   case 'fwd_model'; fwd_model = inv_model;
0030   otherwise; error('PRIOR_LAPLACE requires input type of inv_model or fwd_model');
0031 end
0032 
0033 Reg = eidors_cache(@build_laplace, fwd_model, 'prior_laplace');
0034 
0035 function Reg = build_laplace(fwd_model)
0036 
0037    pp= fwd_model_parameters( fwd_model );
0038    Reg = speye( pp.n_elem );
0039 
0040    Iidx= [];
0041    Jidx= [];
0042    Vidx= [];
0043    for ii=1:pp.n_elem
0044       el_adj = find_adjoin( ii, pp.ELEM );
0045       for jj=el_adj(:)'
0046          Iidx= [Iidx, ii, ii, jj, jj];
0047          Jidx= [Jidx, ii, jj, ii, jj];
0048          Vidx= [Vidx,  1, -1, -1,  1];
0049       end
0050    end
0051 
0052    Reg = sparse(Iidx,Jidx, Vidx, pp.n_elem, pp.n_elem );
0053    
0054    
0055 % find elems which are connected to elems ee
0056 function elems= find_adjoin(ee, ELEM)
0057    nn= ELEM(:,ee);
0058    [d,e]= size(ELEM);
0059    ss = false(size(ELEM));
0060    for i=1:d
0061      ss= ss | ELEM==nn(i);
0062    end
0063    elems= find(sum(ss,1)==d-1);
0064 
0065 function do_unit_test
0066 
0067    imdl = mk_common_model('a2c2',16);
0068    RtR = prior_laplace( imdl );
0069    subplot(221); spy(RtR);
0070    unit_test_cmp('a2c2', nnz(RtR), 240);
0071 
0072    fmdl = mk_circ_tank(2,[],4);
0073    RtR = prior_laplace( fmdl );
0074    subplot(222); spy(RtR);
0075    unit_test_cmp('2-4',RtR(1:4,1:8), [ ...
0076       6, -2,  0, -2, -2,  0,  0,  0; -2,  6, -2,  0,  0, -2,  0,  0;
0077       0, -2,  6, -2,  0,  0, -2,  0; -2,  0, -2,  6,  0,  0,  0, -2]);

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005