jacobian_movement

PURPOSE ^

JACOBIAN_MOVEMENT Computes the Jacobian matrix for conductivity and

SYNOPSIS ^

function J = jacobian_movement(fwd_model, img)

DESCRIPTION ^

 JACOBIAN_MOVEMENT   Computes the Jacobian matrix for conductivity and
 electrode movement variables in 3D EIT.
 Args:     fwd_model - the EIDORS object forward model
            img - the image background conductivity

 fwd_model.conductivity_jacobian - function to calculate conductivity
                                   Jacobian (defaults to jacobian_adjoint)

 Returns:          J - the Jacobian matrix [Jc, Jm]

 WARNING: THIS CODE IS EXPERIMENTAL AND GIVES PROBLEMS
 SEE: Camille Gomez-Laberge, Andy Adler
 Direct EIT Jacobian calculations for conductivity change
  and electrode movement,  Physiol. Meas., 29:S89-S99, 208

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function J = jacobian_movement(fwd_model, img)
0002 % JACOBIAN_MOVEMENT   Computes the Jacobian matrix for conductivity and
0003 % electrode movement variables in 3D EIT.
0004 % Args:     fwd_model - the EIDORS object forward model
0005 %            img - the image background conductivity
0006 %
0007 % fwd_model.conductivity_jacobian - function to calculate conductivity
0008 %                                   Jacobian (defaults to jacobian_adjoint)
0009 %
0010 % Returns:          J - the Jacobian matrix [Jc, Jm]
0011 %
0012 % WARNING: THIS CODE IS EXPERIMENTAL AND GIVES PROBLEMS
0013 % SEE: Camille Gomez-Laberge, Andy Adler
0014 % Direct EIT Jacobian calculations for conductivity change
0015 %  and electrode movement,  Physiol. Meas., 29:S89-S99, 208
0016 
0017 % (C) 2007, Camille Gomez-Laberge and Andy Adler.
0018 %  License: GPL version 2 or version 3
0019 % $Id: jacobian_movement.m 5196 2016-03-04 07:12:05Z alistair_boyle $
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 % System matrix and its parameters
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 %Tensor product for conductivity matrix %MC 25/05/2012
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); % we normalize on our own
0058    Jc = jacobian_adjoint(img);
0059 %    Jc = calc_conductivity_jacobian(pp, fwd_model, img);
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 % Define the element connectivity matrix Ce
0074 function Ce= connectivity_matrix( pp );
0075 lengthX = (pp.n_dims+1)*pp.n_elem;
0076 %%%%%lengthY = pp.n_node; %MC 11/05/2012
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 % Calculate fwd_solution and Impedance mapper matrices
0085 function [Vc, Re] = Vc_Re_matrices( pp, fwd_model, s_mat );
0086 % Define the stimulation matrix Vc for each node I, and pattern
0087 % Ground node is never excited; remove it from the index
0088 %%%%%%nodeidx = 1:pp.n_node; %MC 11/05/2012
0089 nodeidx = 1:size(s_mat);
0090 nodeidx( fwd_model.gnd_node ) = [];
0091 
0092 % The stimulation matrix Vc is the voltage at each node (row) for a
0093 % stimulation (column)
0094 Vc = zeros(pp.n_node, pp.n_stim);
0095 Vc(nodeidx, :) = s_mat(nodeidx, nodeidx) \ pp.QQ(nodeidx,:);
0096 
0097 % Define the electrode resistance matrix Re
0098 % Calculate the resistance between electrodes (row) and all nodes (col)
0099 % N2E matrix maps each electrode to its node(s); we exclude GND
0100 Re = zeros(pp.n_elec, pp.n_node);
0101 Re(:, nodeidx) = pp.N2E(:, nodeidx) / s_mat(nodeidx, nodeidx);
0102 
0103 % FIXME: why do we calculate the negative??
0104 Re = -Re;
0105 
0106 
0107 
0108 % Calculate Meas jacobian from derivative on nodes
0109 % Input delVc
0110 % Ouput J
0111 function J= nodes_to_stim_jacobian( delV, fwd_model, pp )
0112 
0113 sz_out= size(delV,3);
0114 % Define the conductivity Jacobian Jc
0115 J = zeros(pp.n_meas, sz_out);
0116 % Calculate the Jacobian columns for each stimulation pattern
0117 rowidx = 0;
0118 for j = 1:pp.n_stim
0119     % Get the measurement pattern for the stimulation pattern j
0120     meas_pattern = fwd_model.stimulation(j).meas_pattern;
0121     n_measj = size(meas_pattern, 1);
0122     % Extract the voltage sensitivity for electrode j
0123     delVcj = reshape( delV(:,j,:), pp.n_elec, sz_out);
0124     % Calculate sensitivity block for measurements during stimulation j
0125     J(rowidx+(1:n_measj), :) = meas_pattern*delVcj;
0126     rowidx = rowidx+n_measj;
0127 end
0128 
0129 
0130 
0131 % CONDUCTIVITY JACOBIAN (Based on Andy Adler's 1996 algorithms)
0132 % Define the voltage sensitivity delVc on electrode I, for stimulation
0133 % pattern J, for a change in conductivity of element K as a 3D array
0134 function Jc = calc_conductivity_jacobian(pp, fwd_model, img);
0135 delVc = zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0136 % Calculate the values for the voltage sensitivity for each element
0137 for k = 1:pp.n_elem
0138     if ~mod(k,500)
0139         fprintf('   JC: element # %d\n',k);
0140     end
0141     % Extract the coordinates of the element's four nodes
0142     Ae = pp.NODE(:,pp.ELEM(:,k))';
0143     % Augment Ae by adding a column of ones to invert
0144     Ae = inv([ones(pp.n_dims+1,1), Ae]);
0145     % Define Be as the matrix Ae with row 1 deleted
0146     Be = Ae(2:pp.n_dims+1,:);
0147     % Calculate the system submatrix subSe for the element i
0148     %%%%%subSe = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 11/05/2012
0149     subSe = Be'*Be/pp.dfact/abs(det(Ae)); 
0150     % Select the same submatrix of Ce
0151     subidx = (pp.n_dims+1)*(k-1)+1 : (pp.n_dims+1)*k;
0152     % The system submatrix is given by the product
0153     if ~pp.DEBUG
0154         delVc(:,:,k) = pp.Re * pp.Ce(subidx,:)' * subSe * pp.Ce(subidx,:)...
0155             * pp.Vc;
0156     else
0157         sz= (pp.n_dims+1)*pp.n_elem;
0158         delSe = sparse(sz,sz);
0159         se_idx= (pp.n_dims+1)*k+(-pp.n_dims : 0);
0160         delSe(se_idx, se_idx) = subSe;
0161 
0162         delVc(:,:,k) = pp.Re * pp.Ce' * delSe * pp.Ce * pp.Vc;
0163 
0164         if mod(k,50) == 0
0165             delta=1e-6;
0166             img_delta = img;
0167             img_delta.elem_data(k) = img_delta.elem_data(k) + delta;
0168             ss_mat_delta= calc_unconnected_system_mat( fwd_model, img_delta);
0169             delSe_pert = (ss_mat_delta - pp.ss_mat) / delta;
0170 
0171             if norm(delSe -delSe_pert ,1) > 1e-6
0172                 eidors_msg('delSe calc wrong',1);
0173             end
0174         end
0175     end
0176 
0177 end
0178 Jc= nodes_to_stim_jacobian( delVc, fwd_model, pp );
0179 
0180 
0181 
0182 % MOVEMENT JACOBIAN
0183 function Jm = calc_movement_jacobian(pp, fwd_model, img)
0184 % The movement Jacobian is defined for each coordinate Jm = [Jmx Jmy Jmz]
0185 % Define the voltage sensitivity delVm on electrode I, for stimulation
0186 % pattern J, for a movement of electrode K as a 3D array
0187 delVm = zeros(pp.n_elec, pp.n_stim, pp.n_elec*pp.n_dims);
0188 
0189 % Precalculate
0190 Re_Ce      = pp.Re * pp.Ce';
0191 cond_Ce_Vc = pp.kron_cond * pp.Ce * pp.Vc;
0192 for colidx = 1:pp.n_dims
0193 %    fprintf('   JM: direction # %d\n',colidx);
0194     % Calculate the values for the voltage sensitivity for each electrode
0195     for k = 1:pp.n_elec
0196         % Find which elements touch this electrode
0197         elec_nodes = fwd_model.electrode(k).nodes;
0198  
0199         %for compound electrodes, average jacobian for each node
0200 %         delVm_part = zeros(pp.n_elec,pp.n_stim);
0201         delVm_part = calc_delVm(elec_nodes(:)',pp,fwd_model,img,colidx,...
0202                 Re_Ce, cond_Ce_Vc);
0203 %         delVm_part = sparse(length(pp.kron_cond),length(pp.kron_cond));
0204 %         for each_elec_node= elec_nodes(:)';
0205 %            delVm_part =  delVm_part + ...
0206 %                 calc_delVm(each_elec_node,pp,fwd_model,img,colidx,...
0207 %                 Re_Ce, cond_Ce_Vc);
0208 %         end
0209 %         delVm_part = Re_Ce * delVm_part * cond_Ce_Vc;
0210         %%%%% delVm_part = delVm_part/length(elec_nodes); %MC 25/05/2012
0211 %         delVm_part = delVm_part;
0212 
0213         vm_idx= k + pp.n_elec*(colidx-1);
0214         delVm(:,:,vm_idx) = delVm_part;
0215 
0216         if pp.DEBUG
0217             delta=1e-8;
0218             mdl_delta = fwd_model;
0219             mdl_delta.nodes(elec_nodes, colidx) = ...
0220                 mdl_delta.nodes(elec_nodes, colidx) + delta;
0221             img_delta = img;
0222             img_delta.fwd_model = mdl_delta;
0223             S= calc_system_mat(img_delta); 
0224             S=S.E;
0225 % FIXME: AA+CG 30/1/12
0226             [Vc_delta] = Vc_Re_matrices( pp, mdl_delta, S);
0227             delVm_pert = pp.N2E*(Vc_delta - pp.Vc) / delta;
0228             nn = norm(delVm_part - delVm_pert,1 ); % WHY NEGATIVE?
0229 
0230             %%%%% if nn > 5e-5 ; keyboard; end %MC 25/05/2012
0231             if nn > 5e-3 ; keyboard; end
0232 
0233         end
0234     end
0235 end
0236 Jm= nodes_to_stim_jacobian( delVm, fwd_model, pp );
0237 
0238 
0239 
0240 function delVm=  calc_delVm( elec_nodes_array, pp, fwd_model, img, colidx, Re_Ce, cond_Ce_Vc)
0241 I = []; J=[]; S= [];
0242 for elec_nodes= elec_nodes_array(:)';
0243 [rowidx, elemidx] = find(pp.ELEM == elec_nodes);
0244 % Define the system sensitivity matrix to movement delSm
0245 sz= (pp.n_dims+1)*pp.n_elem;
0246 %delSm = sparse(sz,sz);
0247 % For each touching element, calculate the perturbation
0248 jcount = 1;
0249 for j = elemidx'
0250     % Extract the coordinates of the element's four nodes
0251     Ae = pp.NODE(:,pp.ELEM(:, j))';
0252     % Define the invertible matrix P: augment Ae by adding a
0253     % column of ones to invert
0254     P = [ones(pp.n_dims+1,1), Ae];
0255     Ae = inv(P);
0256     absdetAe = abs(det(Ae));
0257     % Define Be as the matrix Ae with row 1 deleted
0258     Be = Ae(2:pp.n_dims+1,:);
0259     % For this coordinate, perturb P by [rowidx,colidx], which are
0260     % our paper's perturbation vectors [a,b]
0261     a = zeros(pp.n_dims+1,1);
0262     b = a;
0263     a(rowidx(jcount)) = 1;
0264     jcount = jcount + 1;
0265     b(colidx+1) = 1;
0266     % Calculate the system submatrix subSm for the element j by
0267     % asymmetric perturbation of the electrode node k
0268     deldetAe =   1/absdetAe*b'*Ae*a;
0269     delBe = -Ae*a*b'*Ae;
0270     delBe = delBe(2:pp.n_dims+1,:);
0271     %%%%%subSm = 2/pp.dfact*(...
0272     %%%%%    deldetAe*Be'*Be + ...
0273     %%%%%    delBe'*Be/absdetAe + ...
0274     %%%%%    Be'*delBe/absdetAe); %MC 11/05/2012
0275     subSm = 1/pp.dfact*(...
0276         deldetAe*Be'*Be + ...
0277         delBe'*Be/absdetAe + ...
0278         Be'*delBe/absdetAe);
0279 
0280     if pp.DEBUG
0281         delta=1e-8;
0282         %%%%% subSe = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 25/05/2012
0283         subSe = Be'*Be/pp.dfact/abs(det(Ae));
0284         d_NODE= pp.NODE;
0285         d_NODE(colidx,elec_nodes) =  d_NODE(colidx,elec_nodes) + delta;
0286         Ae = d_NODE(:,pp.ELEM(:, j))';
0287         Ae = inv( [ones(pp.n_dims+1,1), Ae] );
0288         absdetAe_pert = abs(det(Ae));
0289         deldetAe_pert = (absdetAe_pert - absdetAe) / delta;
0290         % Define Be as the matrix Ae with row 1 deleted
0291         Be = Ae(2:pp.n_dims+1,:);
0292         %%%%% subSe_delta = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 25/05/2012
0293         subSe_delta = Be'*Be/pp.dfact/abs(det(Ae));
0294         subSm_pert= (subSe_delta - subSe ) / delta;
0295         if norm(subSm_pert - subSm,1) > 1e-5
0296             eidors_msg('subSm calc wrong',1);
0297             dd= (subSm_pert - 2/pp.dfact/absdetAe *...
0298                 (delBe'*Be + Be'*delBe) )./(Be'*Be);
0299 
0300             fprintf('colidx=%d, j=%d std=%6.4f >',...
0301                 colidx,j, std(dd(:)));
0302             keyboard
0303             subSm= subSm_pert;
0304         end
0305     end
0306 
0307     % Embed subSm into delSm such that subSm(1,1) is the
0308     % (4j+1,4j+1) element of delSm
0309     se_idx= (pp.n_dims+1)*j+(-pp.n_dims : 0);
0310     switch pp.n_dims
0311        case 2
0312           Iidx = vertcat(se_idx,se_idx,se_idx);
0313           I = [I Iidx(:)];
0314           J = [J, se_idx,se_idx,se_idx];
0315           
0316        case 3
0317           Iidx = vertcat(se_idx,se_idx,se_idx,se_idx);
0318           I = [I Iidx(:)];
0319           J = [J, se_idx,se_idx,se_idx,se_idx];
0320     end
0321     S = [S subSm(:)];
0322 %     delSm(se_idx, se_idx) = subSm;
0323 end
0324 end
0325 delSm = sparse(I,J,S,sz,sz);
0326 delVm = Re_Ce * delSm * cond_Ce_Vc;
0327 
0328 
0329 % The system submatrix is given by the product where delSm is
0330 % non-zero only in submatrices corresponding to touching elements
0331 %%%%% delVm = pp.Re * pp.Ce' * delSm * pp.Ce * pp.Vc; %MC 25/05/2012
0332 % delVm = Re_Ce * delSm * cond_Ce_Vc;
0333 if pp.DEBUG
0334     delta=1e-8;
0335     mdl_delta = fwd_model;
0336     mdl_delta.nodes(elec_nodes, colidx) = ...
0337         mdl_delta.nodes(elec_nodes, colidx) + delta;
0338     ss_mat_delta= calc_unconnected_system_mat( mdl_delta, img );
0339      delSm_pert = (ss_mat_delta - pp.ss_mat) / delta;
0340     % delSe_pert shound be Ce'*delSe*Ce
0341     if norm(delSm -delSm_pert ,1) > 1e-5
0342         eidors_msg('delSm calc wrong',1);
0343         %%%%% delVm = pp.Re * pp.Ce' * delSm_pert * pp.Ce * pp.Vc; %MC 25/05/2012
0344         delVm = pp.Re * pp.Ce' * delSm_pert * pp.kron_cond * pp.Ce * pp.Vc;
0345     
0346         keyboard
0347     end
0348 end
0349 
0350 
0351 
0352 function SS= calc_unconnected_system_mat( fwd_model, img)
0353 % Calc system matrix for Andy Adler's EIT code
0354 % fwd_model = forward model
0355 % img       = image background for system matrix calc
0356 % s_mat = CC' * SS * conductivites * CC;
0357 % where:
0358 %   SS  = Unconnected system Matrix
0359 %   CC  = Connectivity Matrix
0360 
0361 p= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0362 
0363 d= p.n_dims+1;
0364 e= p.n_elem;
0365 n= p.n_node;
0366 
0367 SSiidx= floor([0:d*e-1]'/d)*d*ones(1,d) + ones(d*e,1)*(1:d) ;
0368 SSjidx= [1:d*e]'*ones(1,d);
0369 SSdata= zeros(d*e,d);
0370 dfact= (d-1)*(d-2); % Valid for d<=3
0371 for j=1:e
0372     a=  inv([ ones(d,1), p.NODE( :, p.ELEM(:,j) )' ]);
0373     idx= d*(j-1)+1 : d*j;
0374     SSdata(idx,1:d)= a(2:d,:)'*a(2:d,:)/dfact/abs(det(a));
0375 
0376 end %for j=1:ELEMs
0377 idx= 1:e*d;
0378 SS= sparse(SSiidx,SSjidx,SSdata) * ...
0379     sparse(idx,idx, img.elem_data(ceil(idx/d)) );
0380 
0381 
0382 function do_unit_test;
0383    unit_test_compare_approaches
0384    unit_test_matrix_derivatives
0385    unit_test_diff_jacobian_b2C_const_cond
0386    unit_test_diff_jacobian_n3r2_const_cond
0387    unit_test_diff_jacobian_b2C_rand_cond
0388    unit_test_diff_jacobian_n3r2_rand_cond
0389   unit_test_3d_inv_solve1
0390   unit_test_3d_inv_solve2
0391    
0392 function unit_test_compare_approaches
0393    inv_model = mk_common_model('d2t2',16);
0394    img  = mk_image( inv_model);
0395    fwd_model = inv_model.fwd_model;
0396 
0397    pp = fwd_model_parameters( fwd_model );
0398    pp.DEBUG = 0;
0399    pp.dfact = factorial(pp.n_dims);
0400    s_mat= calc_system_mat(img );
0401    [pp.Vc, pp.Re] = Vc_Re_matrices( pp, fwd_model, s_mat.E );
0402    pp.Ce= connectivity_matrix( pp );
0403 
0404    Jc1= calc_conductivity_jacobian(pp, fwd_model, img);
0405    Jc2= jacobian_adjoint(fwd_model,img);
0406    unit_test_cmp('Compare J d2t2', Jc1, Jc2, 1e-13);
0407 
0408    fwd_model.normalize_measurements = 1;
0409    unit_test_cmp('Compare J norm', Jc1, Jc2, 1e-13);
0410 
0411 
0412 
0413 % TEST CODE FOR MATRIX DERIVATIVES
0414 function unit_test_matrix_derivatives
0415 
0416 TEST= 'd/dt det(X + t*a*b'')';
0417 d= 1e-8;
0418 X= rand(5);
0419 a=zeros(5,1); a(ceil(5*rand))=1;
0420 b=zeros(5,1); b(ceil(5*rand))=1;
0421 dX_p= (det(X + d*a*b') - det(X) )/d;
0422 dX  = b'*inv(X)*a*det(X);
0423 unit_test_cmp(TEST, dX_p,dX,1e-6);
0424 
0425 TEST= 'd/dt inv(X + t*a*b'')';
0426 dX_p= (inv(X + d*a*b') - inv(X) )/d;
0427 dX  = -inv(X)*a*b'*inv(X);
0428 unit_test_cmp(TEST, dX_p,dX,1e-5);
0429 
0430 % TEST d/dt 1/abs(det(X + t*a*b')) = abs(1/det(X+t*a*b'))
0431 for i=1:10
0432     TEST = sprintf('d/dt abs(1/det(X+t*a*b'')) [%02d]',i);
0433     X= rand(5);
0434     a=zeros(5,1); a(ceil(5*rand))=1;
0435     b=zeros(5,1); b(ceil(5*rand))=1;
0436     dX_p= (1/abs(det(X + d*a*b')) - 1/abs(det(X)) )/d;
0437     dX  = - 1/abs(det(X))*b'*inv(X)*a;
0438     unit_test_cmp(TEST, norm([dX_p-dX]),0, 1e-5*norm(dX));
0439 end
0440 
0441 
0442    
0443 function unit_test_diff_jacobian_b2C_const_cond
0444    TEST= 'J_perturb-J_direct - b2C model (const sigma)';
0445    mdl3dim = mk_common_model( 'b2C' );
0446    img = mk_image(mdl3dim);
0447    J_pert=jacobian_movement_perturb(img);
0448    J_direct =jacobian_movement(img);
0449    unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0450 
0451    
0452 function unit_test_diff_jacobian_n3r2_const_cond
0453    TEST= 'J_perturb-J_direct - n3r2 model (const sigma)';
0454    mdl3dim = mk_common_model( 'n3r2', [16,2] );
0455    img = mk_image(mdl3dim);
0456    J_pert=jacobian_movement_perturb(img);
0457    J_direct =jacobian_movement(img);
0458    unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0459    
0460 function unit_test_diff_jacobian_b2C_rand_cond
0461    TEST= 'J_perturb-J_direct - b2C model (rand sigma)';
0462    mdl3dim = mk_common_model( 'b2C' );
0463    cond=0.5+rand(size(mdl3dim.fwd_model.elems,1),1);
0464    img = mk_image(mdl3dim,cond);
0465    J_direct =jacobian_movement(img);
0466    J_pert=jacobian_movement_perturb(img);   
0467    unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0468    
0469 function unit_test_diff_jacobian_n3r2_rand_cond
0470    TEST= 'J_perturb-J_direct - n3r2 model (rand sigma)';
0471    mdl3dim = mk_common_model( 'n3r2', [16,2] );
0472    cond=0.5+rand(size(mdl3dim.fwd_model.elems,1),1);
0473    img = mk_image(mdl3dim,cond);
0474    J_pert=jacobian_movement_perturb(img);
0475    J_direct =jacobian_movement(img);
0476    unit_test_cmp(TEST, norm([J_pert-J_direct]),0, 1e-5*norm(J_direct));
0477   
0478 function unit_test_3d_inv_solve1
0479    mdl3dim = mk_common_model( 'n3r2', [16,2] );
0480    img = mk_image(mdl3dim);
0481    vh = fwd_solve( img );
0482    mdl3dim.fwd_model.jacobian = @jacobian_movement;
0483 
0484    mdl3dim.RtR_prior = @prior_movement;
0485 
0486    imgM = inv_solve(mdl3dim, vh, vh);
0487 
0488  function unit_test_3d_inv_solve2
0489     fmdl = mk_library_model('adult_male_16el_lungs');
0490     [fmdl.stimulation, fmdl.meas_select] = ...
0491        mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0492     
0493     outline = shape_library('get','adult_male','boundary');
0494     
0495     minnode = min(fmdl.nodes);
0496     maxnode = max(fmdl.nodes);
0497     imgsz = [32 32];
0498     
0499     xgrid = linspace(minnode(1),maxnode(1),imgsz(1));
0500     ygrid = linspace(minnode(2),maxnode(2),imgsz(2));
0501     rmdl = mk_grid_model([],xgrid,ygrid);
0502     
0503     % remove pixels outside the model
0504     x_avg = conv2(xgrid, [1,1]/2,'valid');
0505     y_avg = conv2(ygrid, [1,1]/2,'valid');
0506     [x,y] = ndgrid( x_avg, y_avg);
0507     inside = inpolygon(x(:),y(:),outline(:,1),outline(:,2));
0508     ff = find(~inside);
0509     
0510     rmdl.elems([2*ff, 2*ff-1],:)= [];
0511     rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0512     rmdl.coarse2fine(:,ff)= [];
0513     rmdl.mk_coarse_fine_mapping.f2c_offset = [0 0 0.5];
0514     rmdl.mk_coarse_fine_mapping.z_depth = 0.25;
0515     
0516     
0517     % calculate coarse2fine
0518     fmdl.coarse2fine = mk_coarse_fine_mapping(fmdl,rmdl);
0519     
0520     imdl = select_imdl( fmdl,{'Basic GN dif'});
0521     imdl.rec_model = rmdl;
0522     
0523     img = mk_image(imdl,1);
0524    vh = fwd_solve( img );
0525    imdl.prior_use_fwd_not_rec = 1;
0526    imdl.fwd_model.jacobian = @jacobian_movement;
0527    imdl.RtR_prior = @prior_movement;
0528    imgM = inv_solve(imdl, vh, vh);
0529
Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005