0001 function Reg= exponential_covar_prior( inv_model );
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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');
0034
0035 Reg= inv(Reg);
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
0047
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
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
0060
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
0078 rad= sqrt(pp.VOLUME/pi);
0079 elseif pp.n_dims ==3
0080 rad= (pp.VOLUME*3/4/pi).^(1/3);
0081 elseif pp.n_dims ==1
0082 rad= pp.VOLUME/2;
0083 else
0084 error('dont know what to do with n_dims ==%d',pp.n_dims);
0085 end
0086
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
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
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