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 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
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);
0051 d= size(pp.ELEM,1);
0052 e= size(pp.ELEM,2);
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 );
0073 else
0074 bdy = find_boundary(fwd_model.elems);
0075 end
0076 try
0077 bdy = [bdy;fwd_model.system_mat_fields.CEM_boundary];
0078 end
0079
0080
0081
0082
0083
0084 copt.cache_obj = {fwd_model.nodes,fwd_model.elems,fwd_model.electrode};
0085 copt.fstr = 'calculate_N2E';
0086 [N2E,cem_electrodes] = eidors_cache(@calculate_N2E,{fwd_model, bdy, n_elec, n}, copt);
0087
0088 if p>0
0089 stim = fwd_model.stimulation;
0090 [pp.QQ, pp.VV, pp.n_meas] = calc_QQ_fast(N2E, stim, p);
0091 end
0092
0093
0094 pp.n_elem = e;
0095 pp.n_elec = n_elec;
0096 pp.n_node = n;
0097 pp.n_stim = p;
0098 pp.n_dims = d-1;
0099 pp.N2E = N2E;
0100 pp.boundary = bdy;
0101 pp.normalize = mdl_normalize(fwd_model);
0102
0103
0104
0105 function VOLUME = element_volume( NODE, ELEM, e, d)
0106 VOLUME=zeros(e,1);
0107 ones_d = ones(1,d);
0108 d1fac = prod( 1:d-1 );
0109 if d > size(NODE,1)
0110 for i=1:e
0111 this_elem = NODE(:,ELEM(:,i));
0112 VOLUME(i)= abs(det([ones_d;this_elem])) / d1fac;
0113 end
0114 elseif d == 3
0115 for i=1:e
0116 this_elem = NODE(:,ELEM(:,i));
0117 d12= det([ones_d;this_elem([1,2],:)])^2;
0118 d13= det([ones_d;this_elem([1,3],:)])^2;
0119 d23= det([ones_d;this_elem([2,3],:)])^2;
0120 VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0121 end
0122 elseif d == 2
0123 for i=1:e
0124 this_elem = NODE(:,ELEM(:,i));
0125 d12= det([ones_d;this_elem([1],:)])^2;
0126 d13= det([ones_d;this_elem([2],:)])^2;
0127 d23= det([ones_d;this_elem([3],:)])^2;
0128 VOLUME(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0129 end
0130 else
0131 warning('mesh size not understood when calculating volumes')
0132 VOLUME = NaN;
0133 end
0134
0135
0136 function [N2E,cem_electrodes] = calculate_N2E( fwd_model, bdy, n_elec, n);
0137 cem_electrodes= 0;
0138 N2E = sparse(n_elec, n+n_elec);
0139 for i=1:n_elec
0140 eleci = fwd_model.electrode(i);
0141
0142 if isfield(eleci,'faces') && ~isempty(eleci.faces)
0143 if ~isempty(eleci.nodes)
0144 eidors_msg('Warning: electrode %d has both faces and nodes',i);
0145 end
0146
0147 cem_electrodes = cem_electrodes+1;
0148 N2Ei= 1;
0149 N2Ei_nodes = n+cem_electrodes;
0150
0151 elseif isfield(eleci,'nodes')
0152 elec_nodes = fwd_model.electrode(i).nodes;
0153 [N2Ei,N2Ei_nodes,cem_electrodes] = ...
0154 N2Ei_from_nodes(fwd_model, ...
0155 bdy, elec_nodes, cem_electrodes,n);
0156 else
0157 eidors_msg('Warning: electrode %d has no nodes',i);
0158 break;
0159 end
0160 N2E(i, N2Ei_nodes) = N2Ei;
0161 end
0162 N2E = N2E(:, 1:(n+cem_electrodes));
0163
0164
0165 if all(N2E(find(N2E(:)))==1)
0166 N2E = logical(N2E);
0167 end
0168
0169
0170 function [N2Ei,N2Ei_nodes,cem_electrodes]= N2Ei_from_nodes( ...
0171 fwd_model, bdy, elec_nodes, cem_electrodes,n);
0172 if length(elec_nodes) ==1
0173 N2Ei = 1;
0174 N2Ei_nodes = elec_nodes;
0175 elseif length(elec_nodes) ==0
0176 error('EIDORS:fwd_model_parameters:electrode','zero length electrode specified');
0177 else
0178 bdy_idx= find_electrode_bdy( bdy, [], elec_nodes);
0179
0180 if ~isempty(bdy_idx)
0181 cem_electrodes = cem_electrodes+1;
0182 N2Ei= 1;
0183 N2Ei_nodes = n+cem_electrodes;
0184 else
0185
0186 [bdy_idx,srf_area]= find_electrode_bdy( bdy, ...
0187 fwd_model.nodes, elec_nodes);
0188 s_srf_area = sum(srf_area);
0189 if s_srf_area == 0;
0190 error('Surface area for elec#%d is zero. Is boundary correct?',i);
0191 end
0192 N2Ei = srf_area/s_srf_area;
0193 N2Ei_nodes= elec_nodes;
0194 end
0195 end
0196
0197
0198 function [QQ, VV, n_meas] = calc_QQ_slow(N2E, stim, p)
0199 QQ = sparse(size(N2E,2),p);
0200 VV = sparse(size(N2E,2),p); N2E0 = N2E>0;
0201 n_meas= 0;
0202 for i=1:p
0203 src= zeros(size(N2E,2),1);
0204 try; src = N2E' * stim(i).stim_pattern; end
0205 try; src = src + stim(i).interior_sources; end
0206 if all(size(src) == [1,1]) && src==0
0207 error('no stim_patterns or interior_sources provided for pattern #%d',i);
0208 end
0209
0210 QQ(:,i) = src;
0211 n_meas = n_meas + size(stim(i).meas_pattern,1);
0212
0213 vlt= zeros(size(N2E,2),1);
0214 try; vlt = N2E0' * stim(i).volt_pattern; end
0215 VV(:,i) = vlt;
0216 end
0217
0218 function [QQ, VV, n_meas] = calc_QQ_fast(N2E, stim, p)
0219 try
0220 ncols = arrayfun(@(x) size(x.stim_pattern,2), stim);
0221 end
0222 if any(ncols>1);
0223 str = 'multiple columns in stim_pattern for patterns: ';
0224 error('EIDORS:fwd_model_parameters:stim_pattern', ...
0225 [str, sprintf('#%d ',find(ncols>1))]);
0226 end
0227 idx = 1:p; idx(ncols==0)= [];
0228
0229 QQ = sparse(size(N2E,2),p);
0230 try
0231 QQ(:,idx) = N2E' * horzcat( stim(:).stim_pattern );
0232 end
0233 VV = sparse(size(N2E,2),p);
0234
0235
0236
0237
0238 try
0239 ncols = arrayfun(@(x) size(x.volt_pattern,2), stim);
0240 end
0241 if any(ncols>1);
0242 str = 'multiple columns in volt_pattern for patterns: ';
0243 error('EIDORS:fwd_model_parameters:volt_pattern', ...
0244 [str, sprintf('#%d ',find(ncols>1))]);
0245 end
0246 idx = 1:p; idx(ncols==0)= [];
0247
0248 try
0249 VV(:,idx) = (N2E>0)' * horzcat( stim(:).volt_pattern );
0250 end
0251
0252 n_meas = size(vertcat(stim(:).meas_pattern),1);
0253
0254
0255
0256 function do_unit_test
0257 imdl = mk_common_model('a2c2',16); fmdl = imdl.fwd_model;
0258 pp = fwd_model_parameters(fmdl);
0259 [QQ1, VV1, n1m] = calc_QQ_slow(pp.N2E, fmdl.stimulation, pp.n_stim);
0260 [QQ2, VV2, n2m] = calc_QQ_fast(pp.N2E, fmdl.stimulation, pp.n_stim);
0261 unit_test_cmp('calc_QQ', norm(QQ1-QQ2,'fro') + norm(n1m-n2m), 0, 1e-15);
0262 unit_test_cmp('calc_VV1', norm(VV1,'fro'), 0, 1e-15);
0263 unit_test_cmp('calc_VV2', norm(VV2,'fro'), 0, 1e-15);
0264
0265 for i=1:6;
0266 imdl = mk_common_model('a2C0',4); fmdl = imdl.fwd_model;
0267 switch i
0268 case 1; fmdl.stimulation(3).stim_pattern = fmdl.stimulation(3).stim_pattern*[1,2];
0269 expected_err = 'EIDORS:fwd_model_parameters:stim_pattern';
0270 case 2; fmdl.stimulation(1).stim_pattern = [];
0271 expected_err = ''; expected = zeros(45,4);
0272 expected(42:45,2:4) = [0,0,1;-1,0,0;1,-1,0;0,1,-1]*10;
0273 param = 'QQ';
0274 case 3; fmdl.electrode(1).nodes = [];
0275 expected_err = 'EIDORS:fwd_model_parameters:electrode';
0276 case 4; fmdl.stimulation(1).volt_pattern = [zeros(3,1);6];
0277 expected_err = ''; expected = zeros(45,4); expected(45,1) = 6;
0278 param = 'VV';
0279 case 5; fmdl.stimulation(3).volt_pattern = [ones(4,2)];
0280 expected_err = 'EIDORS:fwd_model_parameters:volt_pattern';
0281 case 6; fmdl.electrode(3).faces = [1,2;2,3;3,1];
0282 fmdl.electrode(3).nodes = [];
0283 expected_err = ''; expected = zeros(4,45);
0284 expected(:,42:45) = eye(4);
0285 param = 'N2E';
0286 end
0287 err= '';
0288 try; pp = fwd_model_parameters(fmdl);
0289 catch e
0290 err= e.identifier;
0291 end
0292 if length(expected_err)>0;
0293 unit_test_cmp(['expected error:',num2str(i)], err, expected_err);
0294 else
0295 unit_test_cmp(['case:',num2str(i)], full(pp.(param)), expected);
0296 end
0297 end