jacobian_movement_halfspace

 JACOBIAN_MOVEMENT_PERTURB: J= jacobian_movement_perturb( 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= jacobian_movement_perturb( fwd_model, img)
0002 % JACOBIAN_MOVEMENT_PERTURB: J= jacobian_movement_perturb( 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: jacobian_movement_halfspace.m 5720 2018-03-28 11:57:45Z aadler $
0017 
0018 if nargin == 1
0019    img= fwd_model;
0020 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0021    warning('EIDORS:DeprecatedInterface', ...
0022       ['Calling JACOBIAN_MOVEMENT with two arguments is deprecated and will cause' ...
0023        ' an error in a future version. First argument ignored.']);
0024    warning off EIDORS:DeprecatedInterface
0025 
0026 end
0027 fwd_model= img.fwd_model;
0028 
0029 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0030 delta= 1e-6; % tests indicate this is a good value
0031              % too high and J is not linear, too low and numerical error
0032 
0033 if isfield(fwd_model,'conductivity_jacobian')
0034    % reroute the calculation through calc_jacobian to correctly process
0035    % eidors_default
0036    tmp = img;
0037    tmp.fwd_model = rmfield(tmp.fwd_model,'conductivity_jacobian');
0038    tmp.fwd_model.jacobian = fwd_model.conductivity_jacobian;
0039    Jc = calc_jacobian(tmp);
0040 %    Jc= feval(fwd_model.conductivity_jacobian, fwd_model, img );
0041 else
0042    fwd_model.jacobian_perturb.delta = delta;
0043    fwd_model = mdl_normalize(fwd_model,0); % we normalize on our own
0044    Jc = jacobian_perturb(fwd_model, img);
0045 %    Jc= conductivity_jacobian_perturb( pp, delta, img );
0046 end
0047 
0048 Jm= movement_jacobian( pp, delta, img );
0049 J= [Jc,Jm];
0050 
0051 % calculate normalized Jacobian
0052 if pp.normalize
0053    data= fwd_solve( img );
0054    J= J ./ (data.meas(:)*ones(1,size(J,2)));
0055 end
0056 
0057 
0058 % Calculate the conductivity jacobian based on a perturbation
0059 % of each element by a delta
0060 % Relative error mean(mean(abs(J-Jx)))/ mean(mean(abs(J)))
0061 %   10^-2   0.00308129369015
0062 %   10^-3   0.00030910899216
0063 %   10^-4   0.00003092078190
0064 %   10^-5   0.00000309281790
0065 %   10^-6   0.00000035468582
0066 %   10^-7   0.00000098672198
0067 %   10^-8   0.00000938262464
0068 %   10^-9   0.00009144743903
0069 
0070 function J= conductivity_jacobian_perturb( pp, delta, img );
0071 
0072 J = zeros( pp.n_meas, pp.n_elem );
0073 
0074 elem_data = img.elem_data;
0075 d0= fwd_solve( img );
0076 for i=1:pp.n_elem
0077    img.elem_data   = elem_data;
0078    img.elem_data(i)= elem_data(i) + delta;
0079    di= fwd_solve( img );
0080    J(:,i) = (1/delta) * (di.meas - d0.meas);
0081 end
0082 
0083 % xy-Movement Jacobian
0084 function J= movement_jacobian( pp, delta, img )
0085 
0086 J = zeros( pp.n_meas, pp.n_elec*pp.n_dims );
0087 
0088 node0= img.fwd_model.nodes;
0089 d0= fwd_solve_halfspace( img );
0090 for d= 1:pp.n_dims
0091    for i= 1:pp.n_elec
0092       idx= img.fwd_model.electrode(i).nodes;
0093 
0094       img.fwd_model.nodes( idx, d)= node0(idx,d) + delta;
0095       di= fwd_solve_halfspace( img );
0096       img.fwd_model.nodes( idx, d)= node0(idx,d);
0097 
0098       J_idx = pp.n_elec*(d-1) + i;
0099       J(:,J_idx) = (1/delta) * (di.meas - d0.meas);
0100    end
0101 end