aa_e_move_jacobian

 AA_E_MOVE_JACOBIAN: J= aa_e_move_jacobian( fwd_model, img)
 Calculate Jacobian Matrix for EIT, based on conductivity
   change and movement of electrodes
 J         = Jacobian matrix
 fwd_model = forward model

 fwd_model.conductivity_jacobian = fcn

 fwd_model.normalize_measurements if param exists, calculate
                                  a Jacobian for normalized
                                  difference measurements
 img = image background for jacobian calc

0001 function J= aa_e_move_jacobian( fwd_model, img)
0002 % AA_E_MOVE_JACOBIAN: J= aa_e_move_jacobian( fwd_model, img)
0003 % Calculate Jacobian Matrix for EIT, based on conductivity
0004 %   change and movement of electrodes
0005 % J         = Jacobian matrix
0006 % fwd_model = forward model
0007 %
0008 % fwd_model.conductivity_jacobian = fcn
0009 %
0010 % fwd_model.normalize_measurements if param exists, calculate
0011 %                                  a Jacobian for normalized
0012 %                                  difference measurements
0013 % img = image background for jacobian calc
0014 
0015 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0016 % $Id: aa_e_move_jacobian.html 2819 2011-09-07 16:43:11Z aadler $
0017 
0018 pp= aa_fwd_parameters( fwd_model );
0019 delta= 1e-6; % tests indicate this is a good value
0020              % too high and J is not linear, too low and numerical error
0021 
0022 if isfield(fwd_model,'conductivity_jacobian')
0023    Jc= feval(fwd_model.conductivity_jacobian, fwd_model, img );
0024 else
0025    Jc= conductivity_jacobian_perturb( pp, delta, img );
0026 end
0027 
0028 Jm= movement_jacobian( pp, delta, img );
0029 J= [Jc,Jm];
0030 
0031 % calculate normalized Jacobian
0032 if pp.normalize
0033    data= fwd_solve( img );
0034    J= J ./ (data.meas(:)*ones(1,size(J,2)));
0035 end
0036 
0037 
0038 % Calculate the conductivity jacobian based on a perturbation
0039 % of each element by a delta
0040 % Relative error mean(mean(abs(J-Jx)))/ mean(mean(abs(J)))
0041 %   10^-2   0.00308129369015
0042 %   10^-3   0.00030910899216
0043 %   10^-4   0.00003092078190
0044 %   10^-5   0.00000309281790
0045 %   10^-6   0.00000035468582
0046 %   10^-7   0.00000098672198
0047 %   10^-8   0.00000938262464
0048 %   10^-9   0.00009144743903
0049 
0050 function J= conductivity_jacobian_perturb( pp, delta, img );
0051 
0052 J = zeros( pp.n_meas, pp.n_elem );
0053 
0054 elem_data = img.elem_data;
0055 d0= fwd_solve( img );
0056 for i=1:pp.n_elem
0057    img.elem_data   = elem_data;
0058    img.elem_data(i)= elem_data(i) + delta;
0059    di= fwd_solve( img );
0060    J(:,i) = (1/delta) * (di.meas - d0.meas);
0061 end
0062 
0063 % xy-Movement Jacobian
0064 function J= movement_jacobian( pp, delta, img )
0065 
0066 J = zeros( pp.n_meas, pp.n_elec*pp.n_dims );
0067 
0068 node0= img.fwd_model.nodes;
0069 d0= fwd_solve( img );
0070 for d= 1:pp.n_dims
0071    for i= 1:pp.n_elec
0072       idx= img.fwd_model.electrode(i).nodes;
0073 
0074       img.fwd_model.nodes( idx, d)= node0(idx,d) + delta;
0075       di= fwd_solve( img );
0076       img.fwd_model.nodes( idx, d)= node0(idx,d);
0077 
0078       J_idx = pp.n_elec*(d-1) + i;
0079       J(:,J_idx) = (1/delta) * (di.meas - d0.meas);
0080    end
0081 end