aa_calc_jacobian

 AA_CALC_JACOBIAN: J= aa_calc_jacobian( fwd_model, img)
 Calculate Jacobian Matrix for current stimulation EIT
 J         = Jacobian matrix
 fwd_model = forward model

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

0001 function J= aa_calc_jacobian( fwd_model, img)
0002 % AA_CALC_JACOBIAN: J= aa_calc_jacobian( fwd_model, img)
0003 % Calculate Jacobian Matrix for current stimulation EIT
0004 % J         = Jacobian matrix
0005 % fwd_model = forward model
0006 %
0007 % fwd_model.normalize_measurements if param exists, calculate
0008 %                                  a Jacobian for normalized
0009 %                                  difference measurements
0010 % img = image background for jacobian calc
0011 
0012 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0013 % $Id: aa_calc_jacobian.html 2819 2011-09-07 16:43:11Z aadler $
0014 
0015 pp= aa_fwd_parameters( fwd_model );
0016 s_mat= calc_system_mat( fwd_model, img );
0017 
0018 d= pp.n_dims+1;
0019 e= pp.n_elem;
0020 n= pp.n_node;
0021 
0022 idx= 1:size(s_mat.E,1);
0023 idx( fwd_model.gnd_node ) = [];
0024 
0025 sv= zeros(n, pp.n_stim );
0026 sv( idx,:) = forward_solver(s_mat.E(idx,idx) , pp.QQ( idx,: ));
0027 
0028 zi2E= zeros(pp.n_elec, n);
0029 zi2E(:, idx)= pp.N2E(:,idx)/ s_mat.E(idx,idx) ;
0030 
0031 FC= aa_system_mat_fields( fwd_model );
0032 
0033 
0034 if isfield(fwd_model,'coarse2fine')
0035    DE = jacobian_calc(pp, zi2E, FC, sv, fwd_model.coarse2fine);
0036    nparam= size(fwd_model.coarse2fine,2);
0037 else
0038    DE = jacobian_calc(pp, zi2E, FC, sv);
0039    nparam= e;
0040 end
0041 
0042 J = zeros( pp.n_meas, nparam );
0043 idx=0;
0044 for j= 1:pp.n_stim
0045    meas_pat= fwd_model.stimulation(j).meas_pattern;
0046    n_meas  = size(meas_pat,1);
0047    DEj = reshape( DE(:,j,:), pp.n_elec, nparam );
0048    J( idx+(1:n_meas),: ) = meas_pat*DEj;
0049    idx= idx+ n_meas;
0050 end
0051 
0052 % calculate normalized Jacobian
0053 if pp.normalize
0054    data= fwd_solve( img );
0055    J= J ./ (data.meas(:)*ones(1,nparam));
0056    
0057 end
0058 
0059 % FIXME: The Jacobian calculated is inversed
0060 J= -J;
0061 
0062 % DE_{i,j,k} is dV_i,j / dS_k
0063 %  where V_i is change in voltage on electrode i for
0064 %        stimulation pattern j
0065 %        S_k is change in conductivity on element k
0066 function DE = jacobian_calc(pp, zi2E, FC, sv, c2f);
0067 d= pp.n_dims+1;
0068 dfact= (d-1)*(d-2); % Valid for d<=3
0069 
0070 do_c2f = ( nargin==5 );
0071 
0072 zi2E_FCt = zi2E * FC';
0073 FC_sv   = FC * sv;
0074 
0075 if ~do_c2f
0076    DE= zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0077    for k= 1:pp.n_elem
0078        idx= (d-1)*(k-1)+1 : (d-1)*k;
0079        dq= zi2E_FCt(:,idx) * FC_sv(idx,:);
0080        DE(:,:,k)= dq;
0081    end
0082 else
0083    DE= zeros(pp.n_elec, pp.n_stim, size(c2f,2) );
0084    if 0 % Code is slower
0085       de= pp.n_elem * (d-1);
0086       for k= 1:size(c2f,2);
0087           chg_col = kron( c2f(:,k), ones(d-1,1));
0088           dDD_dEj = spdiags(chg_col,0, de, de);
0089           dq= zi2E_FCt * dDD_dEj * FC_sv;
0090           DE(:,:,k)= dq;
0091       end
0092    else
0093       de= pp.n_elem * (d-1);
0094       for k= 1:size(c2f,2);
0095           ff = find( c2f(:,k) );
0096           lff= length(ff)*(d-1);
0097           ff1= ones(d-1,1) * ff(:)';
0098           ffd= (d-1)*ff1 + (-(d-2):0)'*ones(1,length(ff));
0099           dDD_dEj = spdiags(c2f(ff1,k), 0, lff, lff);
0100           dq= zi2E_FCt(:,ffd) * dDD_dEj * FC_sv(ffd,:);
0101           DE(:,:,k)= dq;
0102       end
0103    end
0104 end