fwd_model_parameters

PURPOSE ^

FWD_MODEL_PARAMETERS: data= fwd_solve_1st_order( fwd_model, image)

SYNOPSIS ^

function param = fwd_model_parameters( fwd_model, opt )

DESCRIPTION ^

 FWD_MODEL_PARAMETERS: data= fwd_solve_1st_order( fwd_model, image)
   Internal function to extract parameters from a fwd_model
   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_elec     => number of electrodes
   param.n_dims     => dimentions (2= 2D, 3=3D)
   param.n_meas     => number of measurements (total)
   param.boundary   => FEM boundary
   param.NODE       => vertex matrix
   param.ELEM       => connection matrix
   param.QQ         => Current into each NODE (Neuman Boundary Conditions)
   param.VV         => Voltage driven into each NODE (Dirichlet BC - only where QQ is NaN)
   param.YY         => Output Admittance (1/Impedance) of each node current (for each value in QQ)
   param.VOLUME     => Volume (or area) of each element
   param.normalize  => difference measurements normalized?
   param.N2E        => Node to electrode converter
   param.v2meas     => Convert node voltages to measurements

 If the stimulation patterns has a 'interior_sources' field,
   the node current QQ, is set to this value for this stimulation.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function param = fwd_model_parameters( fwd_model, opt )
0002 % FWD_MODEL_PARAMETERS: data= fwd_solve_1st_order( fwd_model, image)
0003 %   Internal function to extract parameters from a fwd_model
0004 %   param.n_elem     => number of elements
0005 %   param.n_elec     => number of electrodes
0006 %   param.n_node     => number of nodes (vertices)
0007 %   param.n_stim     => number of current stimulation patterns
0008 %   param.n_elec     => number of electrodes
0009 %   param.n_dims     => dimentions (2= 2D, 3=3D)
0010 %   param.n_meas     => number of measurements (total)
0011 %   param.boundary   => FEM boundary
0012 %   param.NODE       => vertex matrix
0013 %   param.ELEM       => connection matrix
0014 %   param.QQ         => Current into each NODE (Neuman Boundary Conditions)
0015 %   param.VV         => Voltage driven into each NODE (Dirichlet BC - only where QQ is NaN)
0016 %   param.YY         => Output Admittance (1/Impedance) of each node current (for each value in QQ)
0017 %   param.VOLUME     => Volume (or area) of each element
0018 %   param.normalize  => difference measurements normalized?
0019 %   param.N2E        => Node to electrode converter
0020 %   param.v2meas     => Convert node voltages to measurements
0021 %
0022 % If the stimulation patterns has a 'interior_sources' field,
0023 %   the node current QQ, is set to this value for this stimulation.
0024 
0025 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0026 % $Id: fwd_model_parameters.m 6796 2024-04-20 19:21:29Z aadler $
0027 
0028 if ischar(fwd_model) && strcmp(fwd_model, 'UNIT_TEST'); do_unit_test; return; end
0029 
0030 if nargin < 2
0031    opt.skip_VOLUME = 0;
0032 else
0033    assert(ischar(opt),'opt must be a string');
0034    assert(strcmp(opt,'skip_VOLUME'),'opt can only be ''skip_VOLUME''');
0035    opt = struct;
0036    opt.skip_VOLUME = 1;
0037 end
0038 
0039 copt.fstr = 'fwd_model_parameters';
0040 copt.log_level = 4;
0041 
0042 param = eidors_cache(@calc_param,{fwd_model, opt},copt);
0043 
0044 
0045 % perform actual parameter calculation
0046 function pp= calc_param( fwd_model, opt )
0047 
0048 pp.NODE= fwd_model.nodes';
0049 pp.ELEM= fwd_model.elems';
0050 
0051 n= size(pp.NODE,2);        %NODEs
0052 d= size(pp.ELEM,1);        %dimentions+1
0053 e= size(pp.ELEM,2);        %ELEMents
0054 try
0055    p = length(fwd_model.stimulation );
0056 catch 
0057    p = 0;
0058 end
0059 try
0060    n_elec= length( fwd_model.electrode );
0061 catch
0062    n_elec= 0;
0063    fwd_model.electrode = [];
0064 end
0065 
0066 if ~opt.skip_VOLUME
0067    copt.fstr = 'element_volume';
0068    copt.log_level = 4;
0069    pp.VOLUME= eidors_cache(@element_volume, {pp.NODE, pp.ELEM, e, d}, copt );
0070 end
0071 
0072 if isfield(fwd_model,'boundary')
0073     bdy = double( fwd_model.boundary ); % double because of stupid matlab bugs
0074 else
0075     bdy = find_boundary(fwd_model.elems);
0076 end
0077 try %add system_mat_fields.CEM_boundary if it exists
0078    bdy = [bdy;fwd_model.system_mat_fields.CEM_boundary];
0079 end
0080 
0081 % Matrix to convert Nodes to Electrodes
0082 % Complete electrode model for all electrodes
0083 %  N2E = sparse(1:n_elec, n+ (1:n_elec), 1, n_elec, n+n_elec);
0084 %  pp.QQ= sparse(n+n_elec,p);
0085 copt.cache_obj = {fwd_model.nodes,fwd_model.elems,fwd_model.electrode};
0086 copt.fstr = 'calculate_N2E';
0087 [N2E,cem_electrodes] = eidors_cache(@calculate_N2E,{fwd_model, bdy, n_elec, n}, copt);
0088 
0089 % pack into a parameter return list
0090 pp.n_elem   = e;
0091 pp.n_elec   = n_elec;
0092 pp.n_node   = n;
0093 pp.n_stim   = p;
0094 pp.n_dims   = d-1;
0095 pp.N2E      = N2E;
0096 pp.boundary = bdy;
0097 pp.normalize = mdl_normalize(fwd_model);
0098 
0099 if p>0
0100   stim = fwd_model.stimulation;
0101   [pp.QQ, pp.VV, pp.n_meas] = calc_QQ_fast(N2E, stim, p);
0102 
0103   pp.v2meas   = get_v2meas(pp.n_elec, pp.n_stim, stim);
0104 end
0105 
0106 
0107 function v2meas = get_v2meas(n_elec,n_stim,stim)
0108     v2meas = sparse(n_elec*n_stim,0);
0109     for i=1:n_stim
0110         meas_pat= stim(i).meas_pattern;
0111         [n_meas,n_mpat]  = size(meas_pat);
0112         if n_elec ~= n_mpat
0113            error('EIDORS:fwd_model_parameters:meas_pattern', ...
0114                  'meas_pattern#%d uses %d electrodes, but model has %d', ...
0115                  i, n_mpat, n_elec);
0116         end
0117         n_spat = size(stim(i).stim_pattern,1);
0118         if n_elec ~= n_spat
0119            error('EIDORS:fwd_model_parameters:stim_pattern', ...
0120                  'stim_pattern#%d uses %d electrodes, but model has %d', ...
0121                  i, n_spat, n_elec);
0122         end
0123         v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat.';
0124     end
0125     v2meas = v2meas'; % Conjugate transpose. For derivation
0126                       % see Chapter 5 of Adler&Holder 2021
0127 
0128 
0129 % calculate element volume and surface area
0130 function VOLUME = element_volume( NODE, ELEM, e, d)
0131    VOLUME=zeros(e,1);
0132    ones_d = ones(1,d);
0133    d1fac = prod( 1:d-1 );
0134    if d > size(NODE,1)
0135       for i=1:e
0136           this_elem = NODE(:,ELEM(:,i)); 
0137           VOLUME(i)= abs(det([ones_d;this_elem])) / d1fac;
0138       end
0139    elseif d == 3 % 3D nodes in 2D mesh
0140       for i=1:e
0141           this_elem = NODE(:,ELEM(:,i)); 
0142           d12= det([ones_d;this_elem([1,2],:)])^2;
0143           d13= det([ones_d;this_elem([1,3],:)])^2;
0144           d23= det([ones_d;this_elem([2,3],:)])^2;
0145           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0146       end
0147    elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0148       for i=1:e
0149           this_elem = NODE(:,ELEM(:,i)); 
0150           d12= det([ones_d;this_elem([1],:)])^2;
0151           d13= det([ones_d;this_elem([2],:)])^2;
0152           d23= det([ones_d;this_elem([3],:)])^2;
0153           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0154       end
0155    else
0156       warning('mesh size not understood when calculating volumes')
0157       VOLUME = NaN;
0158    end
0159 
0160 
0161 function [N2E,cem_electrodes] = calculate_N2E( fwd_model, bdy, n_elec, n);
0162    cem_electrodes= 0; % num electrodes part of Compl. Elec Model
0163    N2E = sparse(n_elec, n+n_elec);
0164    % Extra nodes are added if nodes for electronics components are added
0165    % Better to use "instrument" electrodes
0166    try
0167       extra_nodes = fwd_model.extra_nodes;
0168    catch
0169       extra_nodes = 0;
0170    end
0171    for i=1:n_elec
0172       eleci = fwd_model.electrode(i);
0173 % The faces field is used only if
0174       if isfield(eleci,'faces')  && ~isempty(eleci.faces)
0175         if ~isempty(eleci.nodes)
0176            eidors_msg('Warning: electrode %d has both faces and nodes',i);
0177         end
0178         % This is a CEM electrode
0179         cem_electrodes = cem_electrodes+1;
0180         N2Ei= 1;
0181         N2Ei_nodes = n+cem_electrodes;
0182 
0183       elseif isfield(eleci,'nodes') 
0184          elec_nodes = fwd_model.electrode(i).nodes;
0185          if ~ischar(elec_nodes)
0186             [N2Ei,N2Ei_nodes,cem_electrodes] =  ...
0187                 N2Ei_from_nodes(fwd_model, ...
0188                  bdy, elec_nodes, cem_electrodes,n);
0189          elseif strcmp(elec_nodes,'instrument')
0190             extra_nodes = extra_nodes + 1; 
0191             N2Ei_nodes = n+cem_electrodes+extra_nodes;
0192          else
0193             error('String "nodes" must be instrument');
0194          end
0195       else
0196           eidors_msg('Warning: electrode %d has no nodes',i);
0197           break; %Not a real electrode so don't include
0198       end
0199       N2E(i, N2Ei_nodes) = N2Ei;
0200    end
0201    N2E = N2E(:, 1:(n+cem_electrodes+extra_nodes));
0202 
0203 % If N2E can be made a logical (0-1) matrix, do it.
0204    if all(N2E(find(N2E(:)))==1)
0205       N2E = logical(N2E);
0206    end
0207 
0208 
0209 function [N2Ei,N2Ei_nodes,cem_electrodes]= N2Ei_from_nodes( ...
0210       fwd_model, bdy, elec_nodes, cem_electrodes,n);
0211   % Instrument nodes are added here like a CEM node but with no connections,
0212   % Use system_mat_instrument to add
0213   if ischar(elec_nodes) && strcmp(elec_nodes,'instrument')
0214      cem_electrodes = cem_electrodes+1;
0215      N2Ei= 1;
0216      N2Ei_nodes = n+cem_electrodes;
0217      return
0218   end
0219 
0220   if length(elec_nodes) ==1 % point electrode (maybe inside body)
0221      N2Ei = 1;
0222      N2Ei_nodes = elec_nodes;
0223   elseif length(elec_nodes) ==0
0224     error('EIDORS:fwd_model_parameters:electrode','zero length electrode specified');
0225   else
0226      bdy_idx= find_electrode_bdy( bdy, [], elec_nodes);
0227 
0228      if ~isempty(bdy_idx) % CEM electrode
0229         cem_electrodes = cem_electrodes+1;
0230         N2Ei= 1;
0231         N2Ei_nodes = n+cem_electrodes;
0232      else % a set of point electrodes
0233           % FIXME: make current defs between point electrodes and CEMs compatible
0234         [bdy_idx,srf_area]= find_electrode_bdy( bdy, ...
0235                        fwd_model.nodes, elec_nodes);
0236         s_srf_area =  sum(srf_area);
0237         if s_srf_area == 0;
0238            error('Surface area for elec#%d is zero. Is boundary correct?',i);
0239         end
0240         N2Ei = srf_area/s_srf_area;
0241         N2Ei_nodes= elec_nodes;
0242      end
0243   end
0244 
0245 
0246 function [QQ, VV, n_meas] = calc_QQ_slow(N2E, stim, p)
0247    QQ = sparse(size(N2E,2),p);
0248    VV = sparse(size(N2E,2),p); N2E0 = N2E>0;
0249    n_meas= 0; % sum total number of measurements
0250    for i=1:p
0251       stimi = stim(i);
0252       src= zeros(size(N2E,2),1);
0253       if isfield(stimi,'stim_pattern')
0254          src =       N2E' * stimi.stim_pattern;
0255       end
0256       if isfield(stimi,'interior_sources')
0257          src = src + stimi.interior_sources;
0258       end
0259       if all(src(:)==0)
0260          error('EIDORS:fwd_model_parameters:stim_pattern', ...
0261                'no stim_patterns or interior_sources provided for pattern #%d',i);
0262       end
0263        
0264       QQ(:,i) = src;
0265       n_meas = n_meas + size(stimi.meas_pattern,1);
0266 
0267        vlt= zeros(size(N2E,2),1);
0268        if isfield(stimi,'volt_pattern');
0269           vlt =      N2E0' * stimi.volt_pattern;
0270        end
0271        VV(:,i) = vlt;
0272    end
0273 
0274 function [QQ, VV, n_meas] = calc_QQ_fast(N2E, stim, p)
0275    try
0276       ncols = arrayfun(@(x) size(x.stim_pattern,2), stim);
0277    catch
0278       error('EIDORS:fwd_model_parameters:stim_pattern', ...
0279             'EIDORS:fwd_model_parameters stim_pattern not specified');
0280    end
0281    if any(ncols>1);
0282       str = 'multiple columns in stim_pattern for patterns: ';
0283       error('EIDORS:fwd_model_parameters:stim_pattern', ...
0284             [str, sprintf('#%d ',find(ncols>1))]);
0285    end
0286    idx = 1:p; idx(ncols==0)= [];
0287 
0288    QQ = sparse(size(N2E,2),p);
0289    try
0290    QQ(:,idx) = N2E' * horzcat( stim(:).stim_pattern );
0291    end
0292    VV = sparse(size(N2E,2),p);
0293    % For voltages, we just need to know which N2E, not the size
0294 
0295 
0296 
0297    try
0298    ncols = arrayfun(@(x) size(x.volt_pattern,2), stim);
0299    end
0300    if any(ncols>1);
0301       str = 'multiple columns in volt_pattern for patterns: ';
0302       error('EIDORS:fwd_model_parameters:volt_pattern', ...
0303             [str, sprintf('#%d ',find(ncols>1))]);
0304    end
0305    idx = 1:p; idx(ncols==0)= [];
0306 
0307    try
0308    VV(:,idx) = (N2E>0)' * horzcat( stim(:).volt_pattern );
0309    end
0310 
0311    n_meas = size(vertcat(stim(:).meas_pattern),1);
0312 
0313 
0314 
0315 function do_unit_test
0316    imdl = mk_common_model('a2c2',16); fmdl = imdl.fwd_model;
0317    pp = fwd_model_parameters(fmdl);
0318    [QQ1, VV1, n1m] = calc_QQ_slow(pp.N2E, fmdl.stimulation, pp.n_stim);
0319    [QQ2, VV2, n2m] = calc_QQ_fast(pp.N2E, fmdl.stimulation, pp.n_stim);
0320    unit_test_cmp('calc_QQ', norm(QQ1-QQ2,'fro') + norm(n1m-n2m), 0, 1e-15);
0321    unit_test_cmp('calc_VV1', norm(VV1,'fro'), 0, 1e-15);
0322    unit_test_cmp('calc_VV2', norm(VV2,'fro'), 0, 1e-15);
0323 
0324    i=0; while true; i=i+1;
0325       imdl = mk_common_model('a2C0',4); fmdl = imdl.fwd_model;
0326       switch i
0327          case 1; fmdl.stimulation(3).stim_pattern = fmdl.stimulation(3).stim_pattern*[1,2]; 
0328                  expected_err = 'EIDORS:fwd_model_parameters:stim_pattern';
0329          case 2; fmdl.stimulation(1).stim_pattern(:) = 0;
0330                  expected_err = ''; expected = zeros(45,4);
0331                  expected(42:45,2:4) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0332                  param = 'QQ';
0333          case 2; fmdl.stimulation(1).stim_pattern(:) = 0;
0334                  expected_err = ''; expected = zeros(45,4);
0335                  expected(42:45,2:4) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0336                  param = 'QQ';
0337          case 3; fmdl.stimulation(1)= [];
0338                  expected_err = ''; expected = zeros(45,3);
0339                  expected(42:45,1:3) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0340                  param = 'QQ';
0341          case 4; fmdl.stimulation(1).stim_pattern = [];
0342                  expected_err = 'EIDORS:fwd_model_parameters:stim_pattern';
0343          case 5; fmdl.electrode(1).nodes = [];
0344                  expected_err = 'EIDORS:fwd_model_parameters:electrode';
0345          case 6; fmdl.stimulation(1).volt_pattern = [zeros(3,1);6];
0346                  expected_err = ''; expected = zeros(45,4); expected(45,1) = 6;
0347                  param = 'VV';
0348          case 7; fmdl.stimulation(3).volt_pattern = [ones(4,2)];
0349                  expected_err = 'EIDORS:fwd_model_parameters:volt_pattern';
0350          case 8; fmdl.electrode(3).faces = [1,2;2,3;3,1];
0351                  fmdl.electrode(3).nodes = [];
0352                  expected_err = ''; expected = zeros(4,45);
0353                  expected(:,42:45) = eye(4);
0354                  param = 'N2E';
0355          otherwise; break
0356       end
0357       err= '';
0358       try;  pp = fwd_model_parameters(fmdl);
0359       catch e
0360          err= e.identifier;
0361       end
0362       if length(expected_err)>0;
0363          unit_test_cmp(['expected error:',num2str(i)], err, expected_err);
0364       else
0365          unit_test_cmp(['case:',num2str(i)], full(pp.(param)), expected);
0366       end
0367    end

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005