calc_transferimpedance

 CALC_TRANSFERIMPEDANCE: Calculates transfer impedance matrix
 
   transfimp = calc_transferimpedance( img)

 fwd_model is a fwd_model structure
 img       is an image structure

 transfimp calculates electrode voltages based on electrode currents as
     V = transfimp*I
 For example, for 4 electrodes, the voltage across [1,2] when 3A is across [3,4]
    is given by [1,-1,0,0] * transfimp * [0;0;3;-3];

0001 function transfimp = calc_transferimpedance( img)
0002 % CALC_TRANSFERIMPEDANCE: Calculates transfer impedance matrix
0003 %
0004 %   transfimp = calc_transferimpedance( img)
0005 %
0006 % fwd_model is a fwd_model structure
0007 % img       is an image structure
0008 %
0009 % transfimp calculates electrode voltages based on electrode currents as
0010 %     V = transfimp*I
0011 % For example, for 4 electrodes, the voltage across [1,2] when 3A is across [3,4]
0012 %    is given by [1,-1,0,0] * transfimp * [0;0;3;-3];
0013 
0014 % (C) 2006 Bill Lionheart. License: GPL version 2 or version 3
0015 % $Id: calc_transferimpedance.m 6068 2020-03-16 14:21:57Z aadler $
0016 
0017 % create new stim patterns
0018 % stimulate with one ref electrode and then each in turn
0019 % make an equiv set of measurements
0020 
0021 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0022 
0023 copt.cache_obj = {img.fwd_model, img.elem_data};
0024 copt.fstr = 'calc_transferimpedance';
0025 transfimp = eidors_cache(@calc_T, img, copt);
0026 
0027 
0028 function transfimp = calc_T( img);
0029    n_elecs= length( img.fwd_model.electrode );
0030 
0031     [stim_pat, meas_pat]= trigonometric( n_elecs );
0032    %[stim_pat, meas_pat]= electrode_wise( n_elecs );
0033    %[stim_pat, meas_pat]= monopolar( n_elecs );
0034    %[stim_pat, meas_pat]= monopolar_even( n_elecs );
0035    img.fwd_model.stimulation = stim_pat;
0036 
0037    imeas_pat= pinv(meas_pat);
0038 
0039    data = fwd_solve(img);
0040 
0041    sz= length(img.fwd_model.stimulation);
0042    transfimp = reshape( data.meas, sz, sz);
0043    transfimp = imeas_pat * transfimp * imeas_pat';
0044 
0045 function [stim_pat, meas_pat] = trigonometric( n_elecs )
0046     stim_pat = struct;
0047     idx= linspace(0,2*pi,n_elecs+1)'; idx(end)= [];
0048     omega= idx*[1:n_elecs/2];
0049     meas_pat= [cos(omega), sin(omega) ]';
0050     for i=1:n_elecs
0051         stim_pat(i).stimulation='Amp';
0052         stim_pat(i).stim_pattern= meas_pat(i,:)';
0053         stim_pat(i).meas_pattern= meas_pat;
0054     end
0055 
0056 function [stim_pat, meas_pat] = electrode_wise( n_elecs)
0057     stim_pat = struct;
0058     meas_pat= [-ones(n_elecs-1,1), speye(n_elecs-1)];
0059     for i=2:n_elecs
0060         stim_pat(i-1).stimulation='Amp';
0061         stim_pat(i-1).stim_pattern= sparse([1,i],1,[-1,1],n_elecs,1);
0062         stim_pat(i-1).meas_pattern= meas_pat;
0063     end
0064 
0065 function [stim_pat, meas_pat] = monopolar( n_elecs)
0066     stim_pat = struct;
0067     meas_pat= speye(n_elecs);
0068     for i=1:n_elecs
0069         stim_pat(i).stimulation='Amp';
0070         stim_pat(i).stim_pattern= sparse(i,1,1,n_elecs,1);
0071         stim_pat(i).meas_pattern= meas_pat;
0072     end
0073 
0074 function [stim_pat, meas_pat] = monopolar_even( n_elecs)
0075     stim_pat = struct;
0076     meas_pat= eye(n_elecs) - ones(n_elecs)/n_elecs;
0077     for i=1:n_elecs
0078         stim_pat(i).stimulation='Amp';
0079         stim_pat(i).stim_pattern= meas_pat(i,:)';
0080         stim_pat(i).meas_pattern= meas_pat;
0081     end
0082 
0083 function do_unit_test
0084    current = 4; measure=1;
0085    [R,img] = test_2d_resistor(current,measure)
0086    T = calc_transferimpedance(img);
0087    unit_test_cmp('2D resistor',[1,-1]*T*[1;-1],sum(R), 1e-10);
0088    
0089 
0090 function [R,img] = test_2d_resistor(current,measure)
0091    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0092    z_contact= .1; wid = 3; len = 12; 
0093 
0094    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0095    fmdl.electrode(1).nodes = find(fmdl.nodes(:,2) ==   0);
0096    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == len);
0097    [fmdl.electrode(:).z_contact] = deal(z_contact);
0098    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0099    img= mk_image(fmdl,conduc);
0100 
0101    Block_R = len / wid / conduc;
0102    Contact_R = z_contact/wid;
0103    R = [Block_R, 2*Contact_R];