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

 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 %
0021 % If the stimulation patterns has a 'interior_sources' field,
0022 %   the node current QQ, is set to this value for this stimulation.
0023 
0024 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: fwd_model_parameters.m 5457 2017-05-03 00:34:26Z aadler $
0026 
0027 if ischar(fwd_model) && strcmp(fwd_model, 'UNIT_TEST'); do_unit_test; return; end
0028 
0029 if nargin < 2
0030    opt.skip_VOLUME = 0;
0031 else
0032    assert(ischar(opt),'opt must be a string');
0033    assert(strcmp(opt,'skip_VOLUME'),'opt can only be ''skip_VOLUME''');
0034    opt = struct;
0035    opt.skip_VOLUME = 1;
0036 end
0037 
0038 copt.fstr = 'fwd_model_parameters';
0039 copt.log_level = 4;
0040 
0041 param = eidors_cache(@calc_param,{fwd_model, opt},copt);
0042 
0043 
0044 % perform actual parameter calculation
0045 function pp= calc_param( fwd_model, opt )
0046 
0047 pp.NODE= fwd_model.nodes';
0048 pp.ELEM= fwd_model.elems';
0049 
0050 n= size(pp.NODE,2);        %NODEs
0051 d= size(pp.ELEM,1);        %dimentions+1
0052 e= size(pp.ELEM,2);        %ELEMents
0053 try
0054    p = length(fwd_model.stimulation );
0055 catch 
0056    p = 0;
0057 end
0058 try
0059    n_elec= length( fwd_model.electrode );
0060 catch
0061    n_elec= 0;
0062    fwd_model.electrode = [];
0063 end
0064 
0065 if ~opt.skip_VOLUME
0066    copt.fstr = 'element_volume';
0067    copt.log_level = 4;
0068    pp.VOLUME= eidors_cache(@element_volume, {pp.NODE, pp.ELEM, e, d}, copt );
0069 end
0070 
0071 if isfield(fwd_model,'boundary')
0072     bdy = double( fwd_model.boundary ); % double because of stupid matlab bugs
0073 else
0074     bdy = find_boundary(fwd_model.elems);
0075 end
0076 
0077 % Matrix to convert Nodes to Electrodes
0078 % Complete electrode model for all electrodes
0079 %  N2E = sparse(1:n_elec, n+ (1:n_elec), 1, n_elec, n+n_elec);
0080 %  pp.QQ= sparse(n+n_elec,p);
0081 copt.cache_obj = {fwd_model.nodes,fwd_model.elems,fwd_model.electrode};
0082 copt.fstr = 'calculate_N2E';
0083 [N2E,cem_electrodes] = eidors_cache(@calculate_N2E,{fwd_model, bdy, n_elec, n}, copt);
0084 
0085 if p>0
0086   stim = fwd_model.stimulation;
0087   [pp.QQ, pp.VV, pp.n_meas] = calc_QQ_fast(N2E, stim, p);
0088 end
0089 
0090 % pack into a parameter return list
0091 pp.n_elem   = e;
0092 pp.n_elec   = n_elec;
0093 pp.n_node   = n;
0094 pp.n_stim   = p;
0095 pp.n_dims   = d-1;
0096 pp.N2E      = N2E;
0097 pp.boundary = bdy;
0098 pp.normalize = mdl_normalize(fwd_model);
0099 
0100 
0101 % calculate element volume and surface area
0102 function VOLUME = element_volume( NODE, ELEM, e, d)
0103    VOLUME=zeros(e,1);
0104    ones_d = ones(1,d);
0105    d1fac = prod( 1:d-1 );
0106    if d > size(NODE,1)
0107       for i=1:e
0108           this_elem = NODE(:,ELEM(:,i)); 
0109           VOLUME(i)= abs(det([ones_d;this_elem])) / d1fac;
0110       end
0111    elseif d == 3 % 3D nodes in 2D mesh
0112       for i=1:e
0113           this_elem = NODE(:,ELEM(:,i)); 
0114           d12= det([ones_d;this_elem([1,2],:)])^2;
0115           d13= det([ones_d;this_elem([1,3],:)])^2;
0116           d23= det([ones_d;this_elem([2,3],:)])^2;
0117           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0118       end
0119    elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0120       for i=1:e
0121           this_elem = NODE(:,ELEM(:,i)); 
0122           d12= det([ones_d;this_elem([1],:)])^2;
0123           d13= det([ones_d;this_elem([2],:)])^2;
0124           d23= det([ones_d;this_elem([3],:)])^2;
0125           VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0126       end
0127    else
0128       warning('mesh size not understood when calculating volumes')
0129       VOLUME = NaN;
0130    end
0131 
0132 
0133 
0134 function [N2E,cem_electrodes] = calculate_N2E( fwd_model, bdy, n_elec, n);
0135    cem_electrodes= 0; % num electrodes part of Compl. Elec Model
0136    N2E = sparse(n_elec, n+n_elec);
0137    for i=1:n_elec
0138        try
0139            elec_nodes = fwd_model.electrode(i).nodes;
0140        catch
0141            eidors_msg('Warning: electrode %d has no nodes',i);
0142            break; %Not a real electrode so don't include
0143        end
0144        if length(elec_nodes) ==1 % point electrode (maybe inside body)
0145           N2E(i, elec_nodes) = 1;
0146        elseif length(elec_nodes) ==0
0147          error('EIDORS:fwd_model_parameters:electrode','zero length electrode specified');
0148        else
0149           bdy_idx= find_electrode_bdy( bdy, [], elec_nodes);
0150 
0151           if ~isempty(bdy_idx) % CEM electrode
0152              cem_electrodes = cem_electrodes+1;
0153              N2E(i, n+cem_electrodes) =1;
0154           else % point electrodes
0155                % FIXME: make current defs between point electrodes and CEMs compatible
0156              [bdy_idx,srf_area]= find_electrode_bdy( bdy, ...
0157                             fwd_model.nodes, elec_nodes);
0158              N2E(i, elec_nodes) = srf_area/sum(srf_area);
0159           end
0160        end
0161    end
0162    N2E = N2E(:, 1:(n+cem_electrodes));
0163 
0164 % If N2E can be made a logical (0-1) matrix, do it.
0165    if all(N2E(find(N2E(:)))==1)
0166       N2E = logical(N2E);
0167    end
0168 
0169 
0170 function [QQ, VV, n_meas] = calc_QQ_slow(N2E, stim, p)
0171    QQ = sparse(size(N2E,2),p);
0172    VV = sparse(size(N2E,2),p); N2E0 = N2E>0;
0173    n_meas= 0; % sum total number of measurements
0174    for i=1:p
0175        src= zeros(size(N2E,2),1);
0176        try;  src =       N2E' * stim(i).stim_pattern; end
0177        try;  src = src + stim(i).interior_sources;    end
0178        if all(size(src) == [1,1]) && src==0
0179           error('no stim_patterns or interior_sources provided for pattern #%d',i);
0180        end
0181        
0182        QQ(:,i) = src;
0183        n_meas = n_meas + size(stim(i).meas_pattern,1);
0184 
0185        vlt= zeros(size(N2E,2),1);
0186        try;  vlt =      N2E0' * stim(i).volt_pattern; end
0187        VV(:,i) = vlt;
0188    end
0189 
0190 function [QQ, VV, n_meas] = calc_QQ_fast(N2E, stim, p)
0191    try
0192    ncols = arrayfun(@(x) size(x.stim_pattern,2), stim);
0193    end
0194    if any(ncols>1);
0195       str = 'multiple columns in stim_pattern for patterns: ';
0196       error('EIDORS:fwd_model_parameters:stim_pattern', ...
0197             [str, sprintf('#%d ',find(ncols>1))]);
0198    end
0199    idx = 1:p; idx(ncols==0)= [];
0200 
0201    QQ = sparse(size(N2E,2),p);
0202    try
0203    QQ(:,idx) = N2E' * horzcat( stim(:).stim_pattern );
0204    end
0205    VV = sparse(size(N2E,2),p);
0206    % For voltages, we just need to know which N2E, not the size
0207 
0208 
0209 
0210    try
0211    ncols = arrayfun(@(x) size(x.volt_pattern,2), stim);
0212    end
0213    if any(ncols>1);
0214       str = 'multiple columns in volt_pattern for patterns: ';
0215       error('EIDORS:fwd_model_parameters:volt_pattern', ...
0216             [str, sprintf('#%d ',find(ncols>1))]);
0217    end
0218    idx = 1:p; idx(ncols==0)= [];
0219 
0220    try
0221    VV(:,idx) = (N2E>0)' * horzcat( stim(:).volt_pattern );
0222    end
0223 
0224    n_meas = size(vertcat(stim(:).meas_pattern),1);
0225 
0226 
0227 
0228 function do_unit_test
0229    imdl = mk_common_model('a2c2',16); fmdl = imdl.fwd_model;
0230    pp = fwd_model_parameters(fmdl);
0231    [QQ1, VV1, n1m] = calc_QQ_slow(pp.N2E, fmdl.stimulation, pp.n_stim);
0232    [QQ2, VV2, n2m] = calc_QQ_fast(pp.N2E, fmdl.stimulation, pp.n_stim);
0233    unit_test_cmp('calc_QQ', norm(QQ1-QQ2,'fro') + norm(n1m-n2m), 0, 1e-15);
0234    unit_test_cmp('calc_VV1', norm(VV1,'fro'), 0, 1e-15);
0235    unit_test_cmp('calc_VV2', norm(VV2,'fro'), 0, 1e-15);
0236 
0237    for i=1:5;
0238       imdl = mk_common_model('a2C0',4); fmdl = imdl.fwd_model;
0239       switch i
0240          case 1; fmdl.stimulation(3).stim_pattern = fmdl.stimulation(3).stim_pattern*[1,2]; 
0241                  expected_err = 'EIDORS:fwd_model_parameters:stim_pattern';
0242          case 2; fmdl.stimulation(1).stim_pattern = [];
0243                  expected_err = ''; expected = zeros(45,4);
0244                  expected(42:45,2:4) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0245                  param = 'QQ';
0246          case 3; fmdl.electrode(1).nodes = [];
0247                  expected_err = 'EIDORS:fwd_model_parameters:electrode';
0248          case 4; fmdl.stimulation(1).volt_pattern = [zeros(3,1);6];
0249                  expected_err = ''; expected = zeros(45,4); expected(45,1) = 6;
0250                  param = 'VV';
0251          case 5; fmdl.stimulation(3).volt_pattern = [ones(4,2)];
0252                  expected_err = 'EIDORS:fwd_model_parameters:volt_pattern';
0253       end
0254       err= '';
0255       try;  pp = fwd_model_parameters(fmdl);
0256       catch e
0257          err= e.identifier;
0258       end
0259       if length(expected_err)>0;
0260          unit_test_cmp(['expected error:',num2str(i)], err, expected_err);
0261       else
0262          unit_test_cmp(['case:',num2str(i)], full(pp.(param)), expected);
0263       end
0264    end

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005