0001 function param = fwd_model_parameters( fwd_model, opt )
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
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
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);
0052 d= size(pp.ELEM,1);
0053 e= size(pp.ELEM,2);
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 );
0074 else
0075 bdy = find_boundary(fwd_model.elems);
0076 end
0077 try
0078 bdy = [bdy;fwd_model.system_mat_fields.CEM_boundary];
0079 end
0080
0081
0082
0083
0084
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
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';
0126
0127
0128
0129
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
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
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;
0163 N2E = sparse(n_elec, n+n_elec);
0164
0165
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
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
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;
0198 end
0199 N2E(i, N2Ei_nodes) = N2Ei;
0200 end
0201 N2E = N2E(:, 1:(n+cem_electrodes+extra_nodes));
0202
0203
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
0212
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
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)
0229 cem_electrodes = cem_electrodes+1;
0230 N2Ei= 1;
0231 N2Ei_nodes = n+cem_electrodes;
0232 else
0233
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;
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
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