exponential_covar_prior

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