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
0077
0078
0079
0080
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
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
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
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
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;
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;
0143 end
0144 if length(elec_nodes) ==1
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)
0152 cem_electrodes = cem_electrodes+1;
0153 N2E(i, n+cem_electrodes) =1;
0154 else
0155
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
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;
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
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