0001 function J = jacobian_movement(fwd_model, img)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return ; end
0022
0023 if nargin == 1
0024 img= fwd_model;
0025 elseif strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0026 warning('EIDORS:DeprecatedInterface', ...
0027 ['Calling JACOBIAN_MOVEMENT with two arguments is deprecated and will cause' ...
0028 ' an error in a future version. First argument ignored.']);
0029 warning off EIDORS:DeprecatedInterface
0030
0031 end
0032 fwd_model= img.fwd_model;
0033
0034
0035
0036 pp = fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0037 pp.dfact = factorial(pp.n_dims);
0038 pp.DEBUG = 0;
0039 if pp.DEBUG
0040 pp.ss_mat = calc_unconnected_system_mat( fwd_model, img);
0041 pp.fwd_meas =fwd_solve( fwd_model, img);
0042 end
0043
0044
0045 I_nd=speye(pp.n_dims+1);
0046 sigma_mat=spdiags(img.elem_data,0,pp.n_elem,pp.n_elem);
0047 pp.kron_cond=kron(sigma_mat,I_nd);
0048
0049 pp.Ce= connectivity_matrix( pp );
0050 s_mat= calc_system_mat( img );
0051 [pp.Vc, pp.Re] = Vc_Re_matrices( pp, fwd_model, s_mat.E );
0052
0053 if isfield(fwd_model,'conductivity_jacobian')
0054 img.fwd_model.jacobian = fwd_model.conductivity_jacobian;
0055 Jc= calc_jacobian( img );
0056 else
0057 img.fwd_model = mdl_normalize(fwd_model, 0);
0058 Jc = jacobian_adjoint(img);
0059
0060 end
0061
0062
0063 Jm = calc_movement_jacobian(pp, fwd_model, img);
0064 J=[Jc,Jm];
0065
0066 if pp.normalize
0067 data= fwd_solve( img );
0068 J= J ./ (data.meas(:)*ones(1,size(J,2)));
0069 end
0070
0071
0072
0073
0074 function Ce= connectivity_matrix( pp );
0075 lengthX = (pp.n_dims+1)*pp.n_elem;
0076
0077 lengthY=size(pp.N2E,2);
0078 Xidx = pp.ELEM(:);
0079 Yidx = ones(lengthX, 1);
0080 Ce = sparse(1:lengthX, Xidx, Yidx, lengthX, lengthY);
0081
0082
0083
0084
0085 function [Vc, Re] = Vc_Re_matrices( pp, fwd_model, s_mat );
0086
0087
0088
0089 assert(size(s_mat,1) == size(s_mat,2),'s_mat must be square');
0090 nodeidx = 1:size(s_mat,1);
0091 nodeidx( fwd_model.gnd_node ) = [];
0092
0093
0094
0095 Vc = zeros(pp.n_node, pp.n_stim);
0096 Vc(nodeidx, :) = s_mat(nodeidx, nodeidx) \ pp.QQ(nodeidx,:);
0097
0098
0099
0100
0101 Re = zeros(pp.n_elec, pp.n_node);
0102 Re(:, nodeidx) = pp.N2E(:, nodeidx) / s_mat(nodeidx, nodeidx);
0103
0104
0105 Re = -Re;
0106
0107
0108
0109
0110
0111
0112 function J= nodes_to_stim_jacobian( delV, fwd_model, pp )
0113
0114 sz_out= size(delV,3);
0115
0116 J = zeros(pp.n_meas, sz_out);
0117
0118 rowidx = 0;
0119 for j = 1:pp.n_stim
0120
0121 meas_pattern = fwd_model.stimulation(j).meas_pattern;
0122 n_measj = size(meas_pattern, 1);
0123
0124 delVcj = reshape( delV(:,j,:), pp.n_elec, sz_out);
0125
0126 J(rowidx+(1:n_measj), :) = meas_pattern*delVcj;
0127 rowidx = rowidx+n_measj;
0128 end
0129
0130
0131
0132
0133
0134
0135 function Jc = calc_conductivity_jacobian(pp, fwd_model, img);
0136 delVc = zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0137
0138 for k = 1:pp.n_elem
0139 if ~mod(k,500)
0140 fprintf(' JC: element # %d\n',k);
0141 end
0142
0143 Ae = pp.NODE(:,pp.ELEM(:,k))';
0144
0145 Ae = inv([ones(pp.n_dims+1,1), Ae]);
0146
0147 Be = Ae(2:pp.n_dims+1,:);
0148
0149
0150 subSe = Be'*Be/pp.dfact/abs(det(Ae));
0151
0152 subidx = (pp.n_dims+1)*(k-1)+1 : (pp.n_dims+1)*k;
0153
0154 if ~pp.DEBUG
0155 delVc(:,:,k) = pp.Re * pp.Ce(subidx,:)' * subSe * pp.Ce(subidx,:)...
0156 * pp.Vc;
0157 else
0158 sz= (pp.n_dims+1)*pp.n_elem;
0159 delSe = sparse(sz,sz);
0160 se_idx= (pp.n_dims+1)*k+(-pp.n_dims : 0);
0161 delSe(se_idx, se_idx) = subSe;
0162
0163 delVc(:,:,k) = pp.Re * pp.Ce' * delSe * pp.Ce * pp.Vc;
0164
0165 if mod(k,50) == 0
0166 delta=1e-6;
0167 img_delta = img;
0168 img_delta.elem_data(k) = img_delta.elem_data(k) + delta;
0169 ss_mat_delta= calc_unconnected_system_mat( fwd_model, img_delta);
0170 delSe_pert = (ss_mat_delta - pp.ss_mat) / delta;
0171
0172 if norm(delSe -delSe_pert ,1) > 1e-6
0173 eidors_msg('delSe calc wrong',1);
0174 end
0175 end
0176 end
0177
0178 end
0179 Jc= nodes_to_stim_jacobian( delVc, fwd_model, pp );
0180
0181
0182
0183
0184 function Jm = calc_movement_jacobian(pp, fwd_model, img)
0185
0186
0187
0188 delVm = zeros(pp.n_elec, pp.n_stim, pp.n_elec*pp.n_dims);
0189
0190
0191 Re_Ce = pp.Re * pp.Ce';
0192 cond_Ce_Vc = pp.kron_cond * pp.Ce * pp.Vc;
0193 for colidx = 1:pp.n_dims
0194
0195
0196 for k = 1:pp.n_elec
0197
0198 elec_nodes = fwd_model.electrode(k).nodes;
0199
0200
0201
0202 delVm_part = calc_delVm(elec_nodes(:)',pp,fwd_model,img,colidx,...
0203 Re_Ce, cond_Ce_Vc);
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214 vm_idx= k + pp.n_elec*(colidx-1);
0215 delVm(:,:,vm_idx) = delVm_part;
0216
0217 if pp.DEBUG
0218 delta=1e-8;
0219 mdl_delta = fwd_model;
0220 mdl_delta.nodes(elec_nodes, colidx) = ...
0221 mdl_delta.nodes(elec_nodes, colidx) + delta;
0222 img_delta = img;
0223 img_delta.fwd_model = mdl_delta;
0224 S= calc_system_mat(img_delta);
0225 S=S.E;
0226
0227 [Vc_delta] = Vc_Re_matrices( pp, mdl_delta, S);
0228 delVm_pert = pp.N2E*(Vc_delta - pp.Vc) / delta;
0229 nn = norm(delVm_part - delVm_pert,1 );
0230
0231
0232 if nn > 5e-3 ; keyboard; end
0233
0234 end
0235 end
0236 end
0237 Jm= nodes_to_stim_jacobian( delVm, fwd_model, pp );
0238
0239
0240
0241 function delVm= calc_delVm( elec_nodes_array, pp, fwd_model, img, colidx, Re_Ce, cond_Ce_Vc)
0242 I = []; J=[]; S= [];
0243 for elec_nodes= elec_nodes_array(:)';
0244 [rowidx, elemidx] = find(pp.ELEM == elec_nodes);
0245
0246 sz= (pp.n_dims+1)*pp.n_elem;
0247
0248
0249 jcount = 1;
0250 for j = elemidx'
0251
0252 Ae = pp.NODE(:,pp.ELEM(:, j))';
0253
0254
0255 P = [ones(pp.n_dims+1,1), Ae];
0256 Ae = inv(P);
0257 absdetAe = abs(det(Ae));
0258
0259 Be = Ae(2:pp.n_dims+1,:);
0260
0261
0262 a = zeros(pp.n_dims+1,1);
0263 b = a;
0264 a(rowidx(jcount)) = 1;
0265 jcount = jcount + 1;
0266 b(colidx+1) = 1;
0267
0268
0269 deldetAe = 1/absdetAe*b'*Ae*a;
0270 delBe = -Ae*a*b'*Ae;
0271 delBe = delBe(2:pp.n_dims+1,:);
0272
0273
0274
0275
0276 subSm = 1/pp.dfact*(...
0277 deldetAe*Be'*Be + ...
0278 delBe'*Be/absdetAe + ...
0279 Be'*delBe/absdetAe);
0280
0281 if pp.DEBUG
0282 delta=1e-8;
0283
0284 subSe = Be'*Be/pp.dfact/abs(det(Ae));
0285 d_NODE= pp.NODE;
0286 d_NODE(colidx,elec_nodes) = d_NODE(colidx,elec_nodes) + delta;
0287 Ae = d_NODE(:,pp.ELEM(:, j))';
0288 Ae = inv( [ones(pp.n_dims+1,1), Ae] );
0289 absdetAe_pert = abs(det(Ae));
0290 deldetAe_pert = (absdetAe_pert - absdetAe) / delta;
0291
0292 Be = Ae(2:pp.n_dims+1,:);
0293
0294 subSe_delta = Be'*Be/pp.dfact/abs(det(Ae));
0295 subSm_pert= (subSe_delta - subSe ) / delta;
0296 if norm(subSm_pert - subSm,1) > 1e-5
0297 eidors_msg('subSm calc wrong',1);
0298 dd= (subSm_pert - 2/pp.dfact/absdetAe *...
0299 (delBe'*Be + Be'*delBe) )./(Be'*Be);
0300
0301 fprintf('colidx=%d, j=%d std=%6.4f >',...
0302 colidx,j, std(dd(:)));
0303 keyboard
0304 subSm= subSm_pert;
0305 end
0306 end
0307
0308
0309
0310 se_idx= (pp.n_dims+1)*j+(-pp.n_dims : 0);
0311 switch pp.n_dims
0312 case 2
0313 Iidx = vertcat(se_idx,se_idx,se_idx);
0314 I = [I Iidx(:)];
0315 J = [J, se_idx,se_idx,se_idx];
0316
0317 case 3
0318 Iidx = vertcat(se_idx,se_idx,se_idx,se_idx);
0319 I = [I Iidx(:)];
0320 J = [J, se_idx,se_idx,se_idx,se_idx];
0321 end
0322 S = [S subSm(:)];
0323
0324 end
0325 end
0326 delSm = sparse(I,J,S,sz,sz);
0327 delVm = Re_Ce * delSm * cond_Ce_Vc;
0328
0329
0330
0331
0332
0333
0334 if pp.DEBUG
0335 delta=1e-8;
0336 mdl_delta = fwd_model;
0337 mdl_delta.nodes(elec_nodes, colidx) = ...
0338 mdl_delta.nodes(elec_nodes, colidx) + delta;
0339 ss_mat_delta= calc_unconnected_system_mat( mdl_delta, img );
0340 delSm_pert = (ss_mat_delta - pp.ss_mat) / delta;
0341
0342 if norm(delSm -delSm_pert ,1) > 1e-5
0343 eidors_msg('delSm calc wrong',1);
0344
0345 delVm = pp.Re * pp.Ce' * delSm_pert * pp.kron_cond * pp.Ce * pp.Vc;
0346
0347 keyboard
0348 end
0349 end
0350
0351
0352
0353 function SS= calc_unconnected_system_mat( fwd_model, img)
0354
0355
0356
0357
0358
0359
0360
0361
0362 p= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0363
0364 d= p.n_dims+1;
0365 e= p.n_elem;
0366 n= p.n_node;
0367
0368 SSiidx= floor([0:d*e-1]'/d)*d*ones(1,d) + ones(d*e,1)*(1:d) ;
0369 SSjidx= [1:d*e]'*ones(1,d);
0370 SSdata= zeros(d*e,d);
0371 dfact= (d-1)*(d-2);
0372 for j=1:e
0373 a= inv([ ones(d,1), p.NODE( :, p.ELEM(:,j) )' ]);
0374 idx= d*(j-1)+1 : d*j;
0375 SSdata(idx,1:d)= a(2:d,:)'*a(2:d,:)/dfact/abs(det(a));
0376
0377 end
0378 idx= 1:e*d;
0379 SS= sparse(SSiidx,SSjidx,SSdata) * ...
0380 sparse(idx,idx, img.elem_data(ceil(idx/d)) );
0381
0382
0383 function do_unit_test;
0384 unit_test_compare_approaches
0385 unit_test_matrix_derivatives
0386 unit_test_diff_jacobian_b2C_const_cond
0387 unit_test_diff_jacobian_n3r2_const_cond
0388 unit_test_diff_jacobian_b2C_rand_cond
0389 unit_test_diff_jacobian_n3r2_rand_cond
0390 unit_test_3d_inv_solve1
0391 unit_test_3d_inv_solve2
0392
0393 function unit_test_compare_approaches
0394 inv_model = mk_common_model('d2t2',16);
0395 img = mk_image( inv_model);
0396 fwd_model = inv_model.fwd_model;
0397
0398 pp = fwd_model_parameters( fwd_model );
0399 pp.DEBUG = 0;
0400 pp.dfact = factorial(pp.n_dims);
0401 s_mat= calc_system_mat(img );
0402 [pp.Vc, pp.Re] = Vc_Re_matrices( pp, fwd_model, s_mat.E );
0403 pp.Ce= connectivity_matrix( pp );
0404
0405 Jc1= calc_conductivity_jacobian(pp, fwd_model, img);
0406 Jc2= jacobian_adjoint(fwd_model,img);
0407 unit_test_cmp('Compare J d2t2', Jc1, Jc2, 1e-13);
0408
0409 fwd_model.normalize_measurements = 1;
0410 unit_test_cmp('Compare J norm', Jc1, Jc2, 1e-13);
0411
0412
0413
0414
0415 function unit_test_matrix_derivatives
0416
0417 TEST= 'd/dt det(X + t*a*b'')';
0418 d= 1e-8;
0419 X= rand(5);
0420 a=zeros(5,1); a(ceil(5*rand))=1;
0421 b=zeros(5,1); b(ceil(5*rand))=1;
0422 dX_p= (det(X + d*a*b') - det(X) )/d;
0423 dX = b'*inv(X)*a*det(X);
0424 unit_test_cmp(TEST, dX_p,dX,1e-6);
0425
0426 TEST= 'd/dt inv(X + t*a*b'')';
0427 dX_p= (inv(X + d*a*b') - inv(X) )/d;
0428 dX = -inv(X)*a*b'*inv(X);
0429 unit_test_cmp(TEST, dX_p,dX,1e-5);
0430
0431
0432 for i=1:10
0433 TEST = sprintf('d/dt abs(1/det(X+t*a*b'')) [%02d]',i);
0434 X= rand(5);
0435 a=zeros(5,1); a(ceil(5*rand))=1;
0436 b=zeros(5,1); b(ceil(5*rand))=1;
0437 dX_p= (1/abs(det(X + d*a*b')) - 1/abs(det(X)) )/d;
0438 dX = - 1/abs(det(X))*b'*inv(X)*a;
0439 unit_test_cmp(TEST, norm([dX_p-dX]),0, 1e-5*norm(dX));
0440 end
0441
0442
0443
0444 function unit_test_diff_jacobian_b2C_const_cond
0445 TEST= 'J_perturb-J_direct - b2C model (const sigma)';
0446 mdl3dim = mk_common_model( 'b2C' );
0447 img = mk_image(mdl3dim);
0448 J_pert=jacobian_movement_perturb(img);
0449 J_direct =jacobian_movement(img);
0450 unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0451
0452
0453 function unit_test_diff_jacobian_n3r2_const_cond
0454 TEST= 'J_perturb-J_direct - n3r2 model (const sigma)';
0455 mdl3dim = mk_common_model( 'n3r2', [16,2] );
0456 img = mk_image(mdl3dim);
0457 J_pert=jacobian_movement_perturb(img);
0458 J_direct =jacobian_movement(img);
0459 unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0460
0461 function unit_test_diff_jacobian_b2C_rand_cond
0462 TEST= 'J_perturb-J_direct - b2C model (rand sigma)';
0463 mdl3dim = mk_common_model( 'b2C' );
0464 cond=0.5+rand(size(mdl3dim.fwd_model.elems,1),1);
0465 img = mk_image(mdl3dim,cond);
0466 J_direct =jacobian_movement(img);
0467 J_pert=jacobian_movement_perturb(img);
0468 unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0469
0470 function unit_test_diff_jacobian_n3r2_rand_cond
0471 TEST= 'J_perturb-J_direct - n3r2 model (rand sigma)';
0472 mdl3dim = mk_common_model( 'n3r2', [16,2] );
0473 cond=0.5+rand(size(mdl3dim.fwd_model.elems,1),1);
0474 img = mk_image(mdl3dim,cond);
0475 J_pert=jacobian_movement_perturb(img);
0476 J_direct =jacobian_movement(img);
0477 unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0478
0479 function unit_test_3d_inv_solve1
0480 mdl3dim = mk_common_model( 'n3r2', [16,2] );
0481 img = mk_image(mdl3dim);
0482 vh = fwd_solve( img );
0483 mdl3dim.fwd_model.jacobian = @jacobian_movement;
0484
0485 mdl3dim.RtR_prior = @prior_movement;
0486
0487 imgM = inv_solve(mdl3dim, vh, vh);
0488
0489 function unit_test_3d_inv_solve2
0490 fmdl = mk_library_model('adult_male_16el_lungs');
0491 [fmdl.stimulation, fmdl.meas_select] = ...
0492 mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0493
0494 outline = shape_library('get','adult_male','boundary');
0495
0496 minnode = min(fmdl.nodes);
0497 maxnode = max(fmdl.nodes);
0498 imgsz = [32 32];
0499
0500 xgrid = linspace(minnode(1),maxnode(1),imgsz(1));
0501 ygrid = linspace(minnode(2),maxnode(2),imgsz(2));
0502 rmdl = mk_grid_model([],xgrid,ygrid);
0503
0504
0505 x_avg = conv2(xgrid, [1,1]/2,'valid');
0506 y_avg = conv2(ygrid, [1,1]/2,'valid');
0507 [x,y] = ndgrid( x_avg, y_avg);
0508 inside = inpolygon(x(:),y(:),outline(:,1),outline(:,2));
0509 ff = find(~inside);
0510
0511 rmdl.elems([2*ff, 2*ff-1],:)= [];
0512 rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0513 rmdl.coarse2fine(:,ff)= [];
0514 rmdl.mk_coarse_fine_mapping.f2c_offset = [0 0 0.5];
0515 rmdl.mk_coarse_fine_mapping.z_depth = 0.25;
0516
0517
0518
0519 fmdl.coarse2fine = mk_coarse_fine_mapping(fmdl,rmdl);
0520
0521 imdl = select_imdl( fmdl,{'Basic GN dif'});
0522 imdl.rec_model = rmdl;
0523
0524 img = mk_image(imdl,1);
0525 vh = fwd_solve( img );
0526 imdl.prior_use_fwd_not_rec = 1;
0527 imdl.fwd_model.jacobian = @jacobian_movement;
0528 imdl.RtR_prior = @prior_movement;
0529 imgM = inv_solve(imdl, vh, vh);
0530