Description of fwd_model

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 6473 2022-12-25 22:41:48Z 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  = size(meas_pat,1);
0112         if n_elec ~= size(meas_pat,2)
0113            error('meas_pattern %d ~= n_elec',i)
0114         end
0115         v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat.';
0116     end
0117     v2meas = v2meas'; % Conjugate transpose. For derivation
0118                       % see Chapter 5 of Adler&Holder 2021
0119 
0120 
0121 % calculate element volume and surface area
0122 function VOLUME = element_volume( NODE, ELEM, e, d)
0123    VOLUME=zeros(e,1);
0124    ones_d = ones(1,d);
0125    d1fac = prod( 1:d-1 );
0126    if d > size(NODE,1)
0127       for i=1:e
0128           this_elem = NODE(:,ELEM(:,i)); 
0129           VOLUME(i)= abs(det([ones_d;this_elem])) / d1fac;
0130       end
0131    elseif d == 3 % 3D nodes in 2D mesh
0132       for i=1:e
0133           this_elem = NODE(:,ELEM(:,i)); 
0134           d12= det([ones_d;this_elem([1,2],:)])^2;
0135           d13= det([ones_d;this_elem([1,3],:)])^2;
0136           d23= det([ones_d;this_elem([2,3],:)])^2;
0137           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0138       end
0139    elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0140       for i=1:e
0141           this_elem = NODE(:,ELEM(:,i)); 
0142           d12= det([ones_d;this_elem([1],:)])^2;
0143           d13= det([ones_d;this_elem([2],:)])^2;
0144           d23= det([ones_d;this_elem([3],:)])^2;
0145           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0146       end
0147    else
0148       warning('mesh size not understood when calculating volumes')
0149       VOLUME = NaN;
0150    end
0151 
0152 
0153 function [N2E,cem_electrodes] = calculate_N2E( fwd_model, bdy, n_elec, n);
0154    cem_electrodes= 0; % num electrodes part of Compl. Elec Model
0155    N2E = sparse(n_elec, n+n_elec);
0156    for i=1:n_elec
0157       eleci = fwd_model.electrode(i);
0158 % The faces field is used only if
0159       if isfield(eleci,'faces')  && ~isempty(eleci.faces)
0160         if ~isempty(eleci.nodes)
0161            eidors_msg('Warning: electrode %d has both faces and nodes',i);
0162         end
0163         % This is a CEM electrode
0164         cem_electrodes = cem_electrodes+1;
0165         N2Ei= 1;
0166         N2Ei_nodes = n+cem_electrodes;
0167 
0168       elseif isfield(eleci,'nodes') 
0169           elec_nodes = fwd_model.electrode(i).nodes;
0170          [N2Ei,N2Ei_nodes,cem_electrodes] =  ...
0171              N2Ei_from_nodes(fwd_model, ...
0172               bdy, elec_nodes, cem_electrodes,n);
0173       else
0174           eidors_msg('Warning: electrode %d has no nodes',i);
0175           break; %Not a real electrode so don't include
0176       end
0177       N2E(i, N2Ei_nodes) = N2Ei;
0178    end
0179    % Extra nodes are added if nodes for electronics components are added
0180    try
0181       extra_nodes = fwd_model.extra_nodes;
0182    catch
0183       extra_nodes = 0;
0184    end
0185    N2E = N2E(:, 1:(n+cem_electrodes+extra_nodes));
0186 
0187 % If N2E can be made a logical (0-1) matrix, do it.
0188    if all(N2E(find(N2E(:)))==1)
0189       N2E = logical(N2E);
0190    end
0191 
0192 
0193 function [N2Ei,N2Ei_nodes,cem_electrodes]= N2Ei_from_nodes( ...
0194       fwd_model, bdy, elec_nodes, cem_electrodes,n);
0195   % Instrument nodes are added here like a CEM node but with no connections,
0196   % Use system_mat_instrument to add
0197   if ischar(elec_nodes) && strcmp(elec_nodes,'instrument')
0198      cem_electrodes = cem_electrodes+1;
0199      N2Ei= 1;
0200      N2Ei_nodes = n+cem_electrodes;
0201      return
0202   end
0203 
0204   if length(elec_nodes) ==1 % point electrode (maybe inside body)
0205      N2Ei = 1;
0206      N2Ei_nodes = elec_nodes;
0207   elseif length(elec_nodes) ==0
0208     error('EIDORS:fwd_model_parameters:electrode','zero length electrode specified');
0209   else
0210      bdy_idx= find_electrode_bdy( bdy, [], elec_nodes);
0211 
0212      if ~isempty(bdy_idx) % CEM electrode
0213         cem_electrodes = cem_electrodes+1;
0214         N2Ei= 1;
0215         N2Ei_nodes = n+cem_electrodes;
0216      else % a set of point electrodes
0217           % FIXME: make current defs between point electrodes and CEMs compatible
0218         [bdy_idx,srf_area]= find_electrode_bdy( bdy, ...
0219                        fwd_model.nodes, elec_nodes);
0220         s_srf_area =  sum(srf_area);
0221         if s_srf_area == 0;
0222            error('Surface area for elec#%d is zero. Is boundary correct?',i);
0223         end
0224         N2Ei = srf_area/s_srf_area;
0225         N2Ei_nodes= elec_nodes;
0226      end
0227   end
0228 
0229 
0230 function [QQ, VV, n_meas] = calc_QQ_slow(N2E, stim, p)
0231    QQ = sparse(size(N2E,2),p);
0232    VV = sparse(size(N2E,2),p); N2E0 = N2E>0;
0233    n_meas= 0; % sum total number of measurements
0234    for i=1:p
0235        src= zeros(size(N2E,2),1);
0236        try;  src =       N2E' * stim(i).stim_pattern; end
0237        try;  src = src + stim(i).interior_sources;    end
0238        if all(size(src) == [1,1]) && src==0
0239           error('no stim_patterns or interior_sources provided for pattern #%d',i);
0240        end
0241        
0242        QQ(:,i) = src;
0243        n_meas = n_meas + size(stim(i).meas_pattern,1);
0244 
0245        vlt= zeros(size(N2E,2),1);
0246        try;  vlt =      N2E0' * stim(i).volt_pattern; end
0247        VV(:,i) = vlt;
0248    end
0249 
0250 function [QQ, VV, n_meas] = calc_QQ_fast(N2E, stim, p)
0251    try
0252       ncols = arrayfun(@(x) size(x.stim_pattern,2), stim);
0253    catch
0254       error('EIDORS:fwd_model_parameters stim_pattern not specified');
0255    end
0256    if any(ncols>1);
0257       str = 'multiple columns in stim_pattern for patterns: ';
0258       error('EIDORS:fwd_model_parameters:stim_pattern', ...
0259             [str, sprintf('#%d ',find(ncols>1))]);
0260    end
0261    idx = 1:p; idx(ncols==0)= [];
0262 
0263    QQ = sparse(size(N2E,2),p);
0264    try
0265    QQ(:,idx) = N2E' * horzcat( stim(:).stim_pattern );
0266    end
0267    VV = sparse(size(N2E,2),p);
0268    % For voltages, we just need to know which N2E, not the size
0269 
0270 
0271 
0272    try
0273    ncols = arrayfun(@(x) size(x.volt_pattern,2), stim);
0274    end
0275    if any(ncols>1);
0276       str = 'multiple columns in volt_pattern for patterns: ';
0277       error('EIDORS:fwd_model_parameters:volt_pattern', ...
0278             [str, sprintf('#%d ',find(ncols>1))]);
0279    end
0280    idx = 1:p; idx(ncols==0)= [];
0281 
0282    try
0283    VV(:,idx) = (N2E>0)' * horzcat( stim(:).volt_pattern );
0284    end
0285 
0286    n_meas = size(vertcat(stim(:).meas_pattern),1);
0287 
0288 
0289 
0290 function do_unit_test
0291    imdl = mk_common_model('a2c2',16); fmdl = imdl.fwd_model;
0292    pp = fwd_model_parameters(fmdl);
0293    [QQ1, VV1, n1m] = calc_QQ_slow(pp.N2E, fmdl.stimulation, pp.n_stim);
0294    [QQ2, VV2, n2m] = calc_QQ_fast(pp.N2E, fmdl.stimulation, pp.n_stim);
0295    unit_test_cmp('calc_QQ', norm(QQ1-QQ2,'fro') + norm(n1m-n2m), 0, 1e-15);
0296    unit_test_cmp('calc_VV1', norm(VV1,'fro'), 0, 1e-15);
0297    unit_test_cmp('calc_VV2', norm(VV2,'fro'), 0, 1e-15);
0298 
0299    for i=1:6;
0300       imdl = mk_common_model('a2C0',4); fmdl = imdl.fwd_model;
0301       switch i
0302          case 1; fmdl.stimulation(3).stim_pattern = fmdl.stimulation(3).stim_pattern*[1,2]; 
0303                  expected_err = 'EIDORS:fwd_model_parameters:stim_pattern';
0304          case 2; fmdl.stimulation(1).stim_pattern = [];
0305                  expected_err = ''; expected = zeros(45,4);
0306                  expected(42:45,2:4) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0307                  param = 'QQ';
0308          case 3; fmdl.electrode(1).nodes = [];
0309                  expected_err = 'EIDORS:fwd_model_parameters:electrode';
0310          case 4; fmdl.stimulation(1).volt_pattern = [zeros(3,1);6];
0311                  expected_err = ''; expected = zeros(45,4); expected(45,1) = 6;
0312                  param = 'VV';
0313          case 5; fmdl.stimulation(3).volt_pattern = [ones(4,2)];
0314                  expected_err = 'EIDORS:fwd_model_parameters:volt_pattern';
0315          case 6; fmdl.electrode(3).faces = [1,2;2,3;3,1];
0316                  fmdl.electrode(3).nodes = [];
0317                  expected_err = ''; expected = zeros(4,45);
0318                  expected(:,42:45) = eye(4);
0319                  param = 'N2E';
0320       end
0321       err= '';
0322       try;  pp = fwd_model_parameters(fmdl);
0323       catch e
0324          err= e.identifier;
0325       end
0326       if length(expected_err)>0;
0327          unit_test_cmp(['expected error:',num2str(i)], err, expected_err);
0328       else
0329          unit_test_cmp(['case:',num2str(i)], full(pp.(param)), expected);
0330       end
0331    end
fwd_model_parameters

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
