np_fwd_parameters

 NP_FWD_PARAMETERS: data= np_fwd_solve( fwd_model )
 Extract parameters from a 'fwd_model' struct which are 
 appropriate for Nick Polydorides EIDORS3D code
   param.n_elem   => number of elements
   param.n_elec   => number of electrodes
   param.n_node   => number of nodes (vertices)
   param.n_stim   => number of current stimulation patterns
   param.n_meas   => number of measurements (total)
   param.vtx      => vertex matrix
   param.simp     => connection matrix
   param.srf      => boundary triangles
   param.df       => vector of measurements for each current pattern
   param.elec     => nodes attached to each electrode
   param.zc       => vector of contact impedances
   param.indH     => electrodes used for each measurement
   param.I        => RHS (current term) for FEM solution
   param.Ib       => Current for electrodes
   param.perm_sym => 'sym' parameter
   param.gnd_ind  => node attached to ground
   param.normalize  => difference measurements normalized?

0001 function param = np_fwd_parameters( fwd_model )
0002 % NP_FWD_PARAMETERS: data= np_fwd_solve( fwd_model )
0003 % Extract parameters from a 'fwd_model' struct which are
0004 % appropriate for Nick Polydorides EIDORS3D code
0005 %   param.n_elem   => number of elements
0006 %   param.n_elec   => number of electrodes
0007 %   param.n_node   => number of nodes (vertices)
0008 %   param.n_stim   => number of current stimulation patterns
0009 %   param.n_meas   => number of measurements (total)
0010 %   param.vtx      => vertex matrix
0011 %   param.simp     => connection matrix
0012 %   param.srf      => boundary triangles
0013 %   param.df       => vector of measurements for each current pattern
0014 %   param.elec     => nodes attached to each electrode
0015 %   param.zc       => vector of contact impedances
0016 %   param.indH     => electrodes used for each measurement
0017 %   param.I        => RHS (current term) for FEM solution
0018 %   param.Ib       => Current for electrodes
0019 %   param.perm_sym => 'sym' parameter
0020 %   param.gnd_ind  => node attached to ground
0021 %   param.normalize  => difference measurements normalized?
0022 
0023 
0024 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: np_fwd_parameters.html 2819 2011-09-07 16:43:11Z aadler $
0026 
0027 param = eidors_obj('get-cache', fwd_model, 'np_fwd_parameters');
0028 
0029 if ~isempty(param)
0030    eidors_msg('np_fwd_parameters: using cached value', 3);
0031    return
0032 end
0033 
0034 param = calc_param( fwd_model );
0035 
0036 eidors_obj('set-cache', fwd_model, 'np_fwd_parameters', param);
0037 eidors_msg('np_fwd_parameters: setting cached value', 3);
0038 
0039 % perform actual parameter calculation
0040 function param= calc_param( fwd_model );
0041 
0042 vtx= fwd_model.nodes;
0043 simp= fwd_model.elems;
0044 % calc num electrodes, nodes, stim_patterns
0045 n_elem= size(simp,1);
0046 n_elec=  length(fwd_model.electrode );
0047 n_node = size(fwd_model.nodes,1);
0048 n_stim = length(fwd_model.stimulation );
0049 n_meas = 0;
0050 
0051 % Recreate 'df' from fwd_model.stimulation
0052 df= zeros(n_stim,1);
0053 for i=1:n_stim;
0054     df(i) = size(fwd_model.stimulation(i).meas_pattern ,1);
0055     n_meas = n_meas + df(i);
0056 end
0057 
0058 elec= [];
0059 zc  = zeros(n_elec, 1);
0060 
0061 if isfield(fwd_model,'boundary')
0062     srf = fwd_model.boundary;
0063 else
0064     srf= find_boundary(simp);
0065 end
0066 
0067 max_elec_nodes=0;
0068 % get electrode parameters
0069 for i=1:n_elec
0070     elec_nodes= fwd_model.electrode(i).nodes;
0071     if length(elec_nodes)>1
0072        e_bdy  = bdy_with_nodes(srf,  elec_nodes );
0073        n_bdy  = srf(e_bdy,:)';
0074     else
0075        n_bdy= elec_nodes;
0076     end
0077 % elec is a series of nodes matching bdy faces
0078     en_list{i}= n_bdy(:)';
0079     if length(n_bdy) > max_elec_nodes
0080         max_elec_nodes = length(n_bdy);
0081     end
0082 
0083 % contact impedance
0084     zc(i)    = fwd_model.electrode(i).z_contact;
0085 end
0086 
0087 elec= zeros(n_elec, max_elec_nodes);
0088 for i=1:n_elec
0089     en= en_list{i};
0090     elec(i,1:length(en)) = en;
0091 end
0092 
0093 % Recreate 'indH' from fwd_model.stimulation
0094 indH= zeros(n_stim, 2);
0095 idx=0;
0096 for i=1:n_stim
0097    meas_pat= fwd_model.stimulation(i).meas_pattern';
0098 
0099    sourcepos= find(meas_pat(:)== 1);
0100    sourcepos= rem( sourcepos-1 , n_elec) + 1;
0101 
0102    sinkpos  = find(meas_pat(:)==-1);
0103    sinkpos  = rem( sinkpos  -1 , n_elec) + 1;
0104 
0105    indH( idx+(1:df(i)) , : ) = [sourcepos, sinkpos];
0106    idx= idx+ df(i);
0107 end
0108 
0109 % calculate FEM RHS matrix, i.e., the current patterns padded with zeroes
0110 I = zeros( n_elec + n_node, n_stim );
0111 idx=0;
0112 for i=1:n_stim
0113    I( n_node + (1:n_elec), i ) = ...
0114          fwd_model.stimulation(i).stim_pattern;
0115 end
0116 I(fwd_model.gnd_node,:) = 0;
0117 Ib= I( n_node + (1:n_elec), : );
0118 
0119 % pack into a parameter return list
0120 param.n_elem   = n_elem;
0121 param.n_elec   = n_elec;
0122 param.n_node   = n_node;
0123 param.n_stim   = n_stim;
0124 param.n_meas   = n_meas;
0125 param.vtx      = vtx;
0126 param.simp     = simp;
0127 param.srf      = srf;
0128 param.df       = df;
0129 param.elec     = elec;
0130 param.zc       = zc;
0131 param.indH     = indH;
0132 param.I        = I;
0133 param.Ib       = Ib;
0134 try
0135    param.perm_sym = fwd_model.np_fwd_solve.perm_sym;
0136 catch
0137    param.perm_sym = '{n}';
0138 end
0139 param.gnd_ind  = fwd_model.gnd_node;
0140 
0141 if isfield(fwd_model,'normalize_measurements')
0142    param.normalize = fwd_model.normalize_measurements;
0143 else
0144    param.normalize = 0;
0145 end
0146 
0147 
0148 
0149 % get boundary faces which match nodes
0150 function e_bdy  = bdy_with_nodes(bdy,  elec_nodes );
0151    mbdy= zeros(size(bdy));
0152    for n= elec_nodes(:)'
0153       mbdy= mbdy + (bdy == n); 
0154    end 
0155    e_bdy = find( all(mbdy') );
0156 
0157 % get boundary faces which match any node
0158 % Use this for point electrodes where there are no bdy faces
0159 % This is sort of an abuse of the model, but at least it can
0160 % produce a reasonable result for pt electrode mdls.
0161 
0162    if isempty(e_bdy)
0163       e_bdy = find( sum(mbdy')>=2 );
0164    end
0165    if isempty(e_bdy)
0166       e_bdy = find( any(mbdy') );
0167    end