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 7009 2024-11-26 14:03:17Z aadler $
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 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 % The stimulation matrix Vc is the voltage at each node (row) for a
0094 % stimulation (column)
0095 Vc = zeros(pp.n_node, pp.n_stim);
0096 Vc(nodeidx, :) = s_mat(nodeidx, nodeidx) \ pp.QQ(nodeidx,:);
0097 
0098 % Define the electrode resistance matrix Re
0099 % Calculate the resistance between electrodes (row) and all nodes (col)
0100 % N2E matrix maps each electrode to its node(s); we exclude GND
0101 Re = zeros(pp.n_elec, pp.n_node);
0102 Re(:, nodeidx) = pp.N2E(:, nodeidx) / s_mat(nodeidx, nodeidx);
0103 
0104 % FIXME: why do we calculate the negative??
0105 Re = -Re;
0106 
0107 
0108 
0109 % Calculate Meas jacobian from derivative on nodes
0110 % Input delVc
0111 % Ouput J
0112 function J= nodes_to_stim_jacobian( delV, fwd_model, pp )
0113 
0114 sz_out= size(delV,3);
0115 % Define the conductivity Jacobian Jc
0116 J = zeros(pp.n_meas, sz_out);
0117 % Calculate the Jacobian columns for each stimulation pattern
0118 rowidx = 0;
0119 for j = 1:pp.n_stim
0120     % Get the measurement pattern for the stimulation pattern j
0121     meas_pattern = fwd_model.stimulation(j).meas_pattern;
0122     n_measj = size(meas_pattern, 1);
0123     % Extract the voltage sensitivity for electrode j
0124     delVcj = reshape( delV(:,j,:), pp.n_elec, sz_out);
0125     % Calculate sensitivity block for measurements during stimulation j
0126     J(rowidx+(1:n_measj), :) = meas_pattern*delVcj;
0127     rowidx = rowidx+n_measj;
0128 end
0129 
0130 
0131 
0132 % CONDUCTIVITY JACOBIAN (Based on Andy Adler's 1996 algorithms)
0133 % Define the voltage sensitivity delVc on electrode I, for stimulation
0134 % pattern J, for a change in conductivity of element K as a 3D array
0135 function Jc = calc_conductivity_jacobian(pp, fwd_model, img);
0136 delVc = zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0137 % Calculate the values for the voltage sensitivity for each element
0138 for k = 1:pp.n_elem
0139     if ~mod(k,500)
0140         fprintf('   JC: element # %d\n',k);
0141     end
0142     % Extract the coordinates of the element's four nodes
0143     Ae = pp.NODE(:,pp.ELEM(:,k))';
0144     % Augment Ae by adding a column of ones to invert
0145     Ae = inv([ones(pp.n_dims+1,1), Ae]);
0146     % Define Be as the matrix Ae with row 1 deleted
0147     Be = Ae(2:pp.n_dims+1,:);
0148     % Calculate the system submatrix subSe for the element i
0149     %%%%%subSe = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 11/05/2012
0150     subSe = Be'*Be/pp.dfact/abs(det(Ae)); 
0151     % Select the same submatrix of Ce
0152     subidx = (pp.n_dims+1)*(k-1)+1 : (pp.n_dims+1)*k;
0153     % The system submatrix is given by the product
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 % MOVEMENT JACOBIAN
0184 function Jm = calc_movement_jacobian(pp, fwd_model, img)
0185 % The movement Jacobian is defined for each coordinate Jm = [Jmx Jmy Jmz]
0186 % Define the voltage sensitivity delVm on electrode I, for stimulation
0187 % pattern J, for a movement of electrode K as a 3D array
0188 delVm = zeros(pp.n_elec, pp.n_stim, pp.n_elec*pp.n_dims);
0189 
0190 % Precalculate
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 %    fprintf('   JM: direction # %d\n',colidx);
0195     % Calculate the values for the voltage sensitivity for each electrode
0196     for k = 1:pp.n_elec
0197         % Find which elements touch this electrode
0198         elec_nodes = fwd_model.electrode(k).nodes;
0199  
0200         %for compound electrodes, average jacobian for each node
0201 %         delVm_part = zeros(pp.n_elec,pp.n_stim);
0202         delVm_part = calc_delVm(elec_nodes(:)',pp,fwd_model,img,colidx,...
0203                 Re_Ce, cond_Ce_Vc);
0204 %         delVm_part = sparse(length(pp.kron_cond),length(pp.kron_cond));
0205 %         for each_elec_node= elec_nodes(:)';
0206 %            delVm_part =  delVm_part + ...
0207 %                 calc_delVm(each_elec_node,pp,fwd_model,img,colidx,...
0208 %                 Re_Ce, cond_Ce_Vc);
0209 %         end
0210 %         delVm_part = Re_Ce * delVm_part * cond_Ce_Vc;
0211         %%%%% delVm_part = delVm_part/length(elec_nodes); %MC 25/05/2012
0212 %         delVm_part = delVm_part;
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 % FIXME: AA+CG 30/1/12
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 ); % WHY NEGATIVE?
0230 
0231             %%%%% if nn > 5e-5 ; keyboard; end %MC 25/05/2012
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 % Define the system sensitivity matrix to movement delSm
0246 sz= (pp.n_dims+1)*pp.n_elem;
0247 %delSm = sparse(sz,sz);
0248 % For each touching element, calculate the perturbation
0249 jcount = 1;
0250 for j = elemidx'
0251     % Extract the coordinates of the element's four nodes
0252     Ae = pp.NODE(:,pp.ELEM(:, j))';
0253     % Define the invertible matrix P: augment Ae by adding a
0254     % column of ones to invert
0255     P = [ones(pp.n_dims+1,1), Ae];
0256     Ae = inv(P);
0257     absdetAe = abs(det(Ae));
0258     % Define Be as the matrix Ae with row 1 deleted
0259     Be = Ae(2:pp.n_dims+1,:);
0260     % For this coordinate, perturb P by [rowidx,colidx], which are
0261     % our paper's perturbation vectors [a,b]
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     % Calculate the system submatrix subSm for the element j by
0268     % asymmetric perturbation of the electrode node k
0269     deldetAe =   1/absdetAe*b'*Ae*a;
0270     delBe = -Ae*a*b'*Ae;
0271     delBe = delBe(2:pp.n_dims+1,:);
0272     %%%%%subSm = 2/pp.dfact*(...
0273     %%%%%    deldetAe*Be'*Be + ...
0274     %%%%%    delBe'*Be/absdetAe + ...
0275     %%%%%    Be'*delBe/absdetAe); %MC 11/05/2012
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         %%%%% subSe = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 25/05/2012
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         % Define Be as the matrix Ae with row 1 deleted
0292         Be = Ae(2:pp.n_dims+1,:);
0293         %%%%% subSe_delta = 2*Be'*Be/pp.dfact/abs(det(Ae)); %MC 25/05/2012
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     % Embed subSm into delSm such that subSm(1,1) is the
0309     % (4j+1,4j+1) element of delSm
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 %     delSm(se_idx, se_idx) = subSm;
0324 end
0325 end
0326 delSm = sparse(I,J,S,sz,sz);
0327 delVm = Re_Ce * delSm * cond_Ce_Vc;
0328 
0329 
0330 % The system submatrix is given by the product where delSm is
0331 % non-zero only in submatrices corresponding to touching elements
0332 %%%%% delVm = pp.Re * pp.Ce' * delSm * pp.Ce * pp.Vc; %MC 25/05/2012
0333 % delVm = Re_Ce * delSm * cond_Ce_Vc;
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     % delSe_pert shound be Ce'*delSe*Ce
0342     if norm(delSm -delSm_pert ,1) > 1e-5
0343         eidors_msg('delSm calc wrong',1);
0344         %%%%% delVm = pp.Re * pp.Ce' * delSm_pert * pp.Ce * pp.Vc; %MC 25/05/2012
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 % Calc system matrix for Andy Adler's EIT code
0355 % fwd_model = forward model
0356 % img       = image background for system matrix calc
0357 % s_mat = CC' * SS * conductivites * CC;
0358 % where:
0359 %   SS  = Unconnected system Matrix
0360 %   CC  = Connectivity Matrix
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); % Valid for d<=3
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 %for j=1:ELEMs
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 % TEST CODE FOR MATRIX DERIVATIVES
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 % TEST d/dt 1/abs(det(X + t*a*b')) = abs(1/det(X+t*a*b'))
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     % remove pixels outside the model
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     % calculate coarse2fine
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

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005