prior_exponential_covar

 PRIOR_EXPONENTIAL_COVAR image prior with exponential
      inter-element covariance, multiplied by a NOSER
      type weighting, if desired.
 Reg= prior_exponential_covar( inv_model )
 Reg        => output regularization term
 inv_model  => inverse model struct
 Parameters: exponential rate
   gamma= inv_model.prior_exponential_covar.gamma
 Parameters: NOSER exponent: diag(J'*J)^exponent )
   gamma= inv_model.prior_exponential_covar.noser_exp
       DEFAULT is 0 (ie. no NOSER exponent)

 FIXME: this code does an inv or the reg matrix. This
        is very inefficient.

0001 function Reg= prior_exponential_covar( inv_model );
0002 % PRIOR_EXPONENTIAL_COVAR image prior with exponential
0003 %      inter-element covariance, multiplied by a NOSER
0004 %      type weighting, if desired.
0005 % Reg= prior_exponential_covar( inv_model )
0006 % Reg        => output regularization term
0007 % inv_model  => inverse model struct
0008 % Parameters: exponential rate
0009 %   gamma= inv_model.prior_exponential_covar.gamma
0010 % Parameters: NOSER exponent: diag(J'*J)^exponent )
0011 %   gamma= inv_model.prior_exponential_covar.noser_exp
0012 %       DEFAULT is 0 (ie. no NOSER exponent)
0013 %
0014 % FIXME: this code does an inv or the reg matrix. This
0015 %        is very inefficient.
0016 
0017 % (C) 2007 Andy Adler. License: GPL version 2 or version 3
0018 % $Id: prior_exponential_covar.m 3109 2012-06-08 15:04:12Z bgrychtol $
0019 
0020 fwd_model= inv_model.fwd_model;
0021 
0022 try 
0023     gamma= inv_model.prior_exponential_covar.gamma;
0024 catch
0025     xy_diam = max( max(fwd_model.nodes(:,1:2)) -  ...
0026                    min(fwd_model.nodes(:,1:2)));
0027     gamma= 0.05*xy_diam;
0028 end
0029 
0030 Reg = calc_prior_exponential_covar( fwd_model, gamma);
0031 Reg= 0.5*(Reg+Reg'); % calculation should be symmetric, but is slightly off.
0032 
0033 Reg= inv(Reg); % FIXME - not part of inverse
0034 
0035 try 
0036     noser_exp= inv_model.prior_exponential_covar.noser_exp;
0037 catch
0038     noser_exp= 0;
0039 end
0040 
0041 if noser_exp>0
0042     
0043     J = calc_jacobian( mk_image( inv_model ));
0044 % If J is too big, then the jacobian includes other terms (like movement) remove them
0045 % TODO: cruft code
0046     n_elem = size(inv_model.fwd_model.elems,1);
0047     if size(J,2) > n_elem; J = J(:,1:n_elem); end
0048     l_prior= size(J,2);
0049 
0050     % Reg is spdiags(diag(J'*J),0, l_prior, l_prior);
0051     diag_col= sum(J.^2,1)';
0052     RegN = spdiags( diag_col.^noser_exp/2, 0, l_prior, l_prior);
0053 
0054     Reg= RegN * Reg * RegN;
0055 end
0056 
0057 % Calculate exponential LPF Filter
0058 % parameter is gamma (normally 0.1)
0059 function Reg= calc_prior_exponential_covar( fwd_model, gamma)
0060    [rad,ctr]= get_elem_rad_ctr( fwd_model );
0061 
0062    n_elem= size(ctr,1);
0063    oo= ones(n_elem,1);
0064    radh = rad/2;
0065    Reg=zeros(n_elem,n_elem);
0066    for i=1:size(ctr,1)
0067       ctr_i = sqrt( sum( (ctr - oo*ctr(i,:)).^2 , 2));
0068       Reg_i= integ_fn(-radh,radh,ctr_i-radh(i),ctr_i+radh(i), gamma);
0069       Reg(:,i)= Reg_i .*(Reg_i>1e-4);
0070    end
0071    Reg=sparse(Reg);
0072 
0073 function [rad,elem_ctr]= get_elem_rad_ctr( fwd_model );
0074    pp= fwd_model_parameters( fwd_model);
0075    if     pp.n_dims==2 % in 2d A=pi*r^2
0076       rad= sqrt(pp.VOLUME/pi);
0077    elseif pp.n_dims ==3 % in 3D V=4/3*pi*r^3
0078       rad= (pp.VOLUME*3/4/pi).^(1/3);
0079    elseif pp.n_dims ==1 % in 1D V=2*r
0080       rad= pp.VOLUME/2;
0081    else 
0082       error('dont know what to do with n_dims ==%d',pp.n_dims);
0083    end
0084 %  rad =   rad* ones(1,pp.n_elem); % copy to matrix
0085 
0086    node_map = fwd_model.nodes(pp.ELEM,:);
0087    elem_ctr = mean( reshape( node_map, pp.n_dims+1, [], pp.n_dims), 1);
0088    elem_ctr = squeeze( elem_ctr);
0089 % show_fem(fwd_model); hold on; plot(elem_ctr(:,1),elem_ctr(:,2),'*'); hold off
0090 
0091 %  ctr= zeros(pp.n_elem);
0092 %  for i=1:pp.n_dims
0093 %     v= elem_ctr(:,i)*ones(1,pp.n_elem); % make square
0094 %     ctr = ctr + (v - v').^2;
0095 %  end
0096 %  ctr = sqrt(ctr);
0097 
0098 % calculate the exponential integral over a space x1,x2
0099 %  given gamma, x1, x2, y1, y2
0100 function fi= integ_fn(x1i,x2i,y1i,y2i, gamma)
0101    i_gam = 1/gamma;
0102 
0103    x1= min(x1i,x2i); x2= max(x1i,x2i);
0104    y1= min(y1i,y2i); y2= max(y1i,y2i);
0105 
0106    Dx= abs(x1-x2); Dy= abs(y1-y2);
0107 
0108    xa1= x1; 
0109    xa2= max(x1,min(x2,y1)); xb1= xa2;
0110    xb2= max(x1,min(x2,y2)); xc1= xb2;
0111    xc2= x2;
0112 
0113    RA= exp(-i_gam*y1).*(1-exp(-i_gam*Dy));
0114    RC= exp( i_gam*y2).*(1-exp(-i_gam*Dy));
0115    fi = 1./Dx./Dy/i_gam^2 .* ( ... 
0116             RA.*(exp( i_gam*xa2) - exp( i_gam*xa1)) - ...
0117             RC.*(exp(-i_gam*xc2) - exp(-i_gam*xc1)) + ...
0118             2*i_gam*(xb2-xb1) - ...
0119             exp(-i_gam*y2).*(exp( i_gam*xb2) - exp( i_gam*xb1)) + ...
0120             exp( i_gam*y1).*(exp(-i_gam*xb2) - exp(-i_gam*xb1))...
0121            );
0122    
0123