simulate_movement

 SIMULATE_MOVEMENT simulate small conductivity perturbations
 [vh,vi,xyzr, c2f]= simulate_movement( img, xyzr );

 Simulate small conductivity perturbations (using the 
  Jacobian calculation to be fast).

 INPUT:
   2D Models: specify xyzr = matrix of 3xN.
      Each perturbation (i) is at (x,y) = xyzr(1:2,i)
      with radius xyzr(3,i). 
  
   3D Models: specify xyzr = matrix of 4xN.
      Each perturbation (i) is at (x,y,z) = xyzr(1:3,i)
      with radius xyzr(4,i). 

   xyzr = scalar =N - single spiral of N in medium centre
  
   img = eidors image object (with img.fwd_model FEM model).

 OUTPUT:
   vh - homogeneous measurements M x 1
   vi - target simulations       M x n_points
   xyzr - x y and radius of each point 3 x n_points
   c2f - image representation of the simulations

 Example: Simulate small 3D object in centre:
    img = mk_image( mk_common_model('b3cr',16) ); % 2D Image
    [vh,vi] = simulate_movement( img, [0;0;0;0.05]);
 Example: Simulate small 2D object in at x near side:
    img = mk_image( mk_common_model('d2d3c',16) ); % 2D Image
    [vh,vi] = simulate_movement( img, [0.9;0;0.05]);
    show_fem(inv_solve(mk_common_model('c2c2',16),vh,vi)) % Reconst (new mesh)

0001 function [vh,vi,xyzr,c2f]= simulate_movement( img, xyzr, UNUSED );
0002 % SIMULATE_MOVEMENT simulate small conductivity perturbations
0003 % [vh,vi,xyzr, c2f]= simulate_movement( img, xyzr );
0004 %
0005 % Simulate small conductivity perturbations (using the
0006 %  Jacobian calculation to be fast).
0007 %
0008 % INPUT:
0009 %   2D Models: specify xyzr = matrix of 3xN.
0010 %      Each perturbation (i) is at (x,y) = xyzr(1:2,i)
0011 %      with radius xyzr(3,i).
0012 %
0013 %   3D Models: specify xyzr = matrix of 4xN.
0014 %      Each perturbation (i) is at (x,y,z) = xyzr(1:3,i)
0015 %      with radius xyzr(4,i).
0016 %
0017 %   xyzr = scalar =N - single spiral of N in medium centre
0018 %
0019 %   img = eidors image object (with img.fwd_model FEM model).
0020 %
0021 % OUTPUT:
0022 %   vh - homogeneous measurements M x 1
0023 %   vi - target simulations       M x n_points
0024 %   xyzr - x y and radius of each point 3 x n_points
0025 %   c2f - image representation of the simulations
0026 %
0027 % Example: Simulate small 3D object in centre:
0028 %    img = mk_image( mk_common_model('b3cr',16) ); % 2D Image
0029 %    [vh,vi] = simulate_movement( img, [0;0;0;0.05]);
0030 % Example: Simulate small 2D object in at x near side:
0031 %    img = mk_image( mk_common_model('d2d3c',16) ); % 2D Image
0032 %    [vh,vi] = simulate_movement( img, [0.9;0;0.05]);
0033 %    show_fem(inv_solve(mk_common_model('c2c2',16),vh,vi)) % Reconst (new mesh)
0034 %
0035 
0036 % (C) 2009 Andy Adler. Licensed under GPL v2 or v3
0037 % $Id: simulate_movement.m 6249 2022-03-24 18:28:16Z aadler $
0038 
0039    if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0040 
0041    if nargin>=3; 
0042       warning('change parameter unused');
0043    end
0044 
0045    if size(xyzr) == [1,1]
0046       path = linspace(0,1,xyzr); phi = 2*pi*path;
0047       meanodes= mean(    img.fwd_model.nodes  );
0048       lennodes= size( img.fwd_model.nodes,1); 
0049       img.fwd_model.nodes = img.fwd_model.nodes - ones(lennodes,1)*meanodes;
0050       maxnodes= max(max(abs( img.fwd_model.nodes(:,1:2) )));
0051       img.fwd_model.nodes = img.fwd_model.nodes / maxnodes;
0052 
0053       xyzr = [0.9*path.*sin(phi);0.9*path.*cos(phi);0*path; 0.05 + 0*path];
0054    end
0055    
0056    [vh,vi,xyzr,c2f]= eidors_cache( @do_simulate_movement, {img, xyzr});
0057 
0058 function [vh,vi,xyzr,c2f]= do_simulate_movement( img, xyzr);
0059 
0060    Nt = size(xyzr,2);
0061    [c2f failed] = mk_c2f_circ_mapping( img.fwd_model, xyzr);
0062    xyzr(:,failed) = [];
0063    c2f(:,failed) = [];
0064    Nt = Nt - nnz(failed);
0065    img.fwd_model.coarse2fine = c2f;
0066 
0067    % We don't want a normalized jacobian here
0068    img.fwd_model = mdl_normalize(img.fwd_model,0);
0069 
0070    J= calc_jacobian( img );
0071    J= move_jacobian_postprocess( J, img, Nt);
0072 
0073    vh= fwd_solve(img);
0074    vh=vh.meas;
0075 
0076    vi= vh*ones(1,Nt) + J;
0077    
0078    % Would this be the slow approach?:
0079 %    vi = vh*zeros(1,Nt);
0080 %    for i = 1: Nt
0081 %        img.elem_data = 1 - c2f(:,i);
0082 %        jnk = fwd_solve(img);
0083 %        vi(:,i) = jnk.meas;
0084 %    end
0085 
0086 function J= move_jacobian_postprocess( J, img, Nt)
0087    if size(J,2) == Nt; % No problem
0088        return ;
0089    % Check if movement jacobian introduced electrode movements (elecs * dims)
0090    elseif size(J,2) == Nt + ...
0091         length(img.fwd_model.electrode) * size(img.fwd_model.nodes,2)
0092       J = J(:,1:Nt);
0093    else
0094       error('Jacobian calculator is not doing the coarse2fine mapping. This is a bug.');
0095    end
0096 
0097 
0098 
0099 
0100 
0101 % QUESTON: the accuracy of J will depend on how well we interpolate the
0102 % mesh. How important is this?
0103 
0104 function do_unit_test
0105    imdl = mk_common_model('a2c0',8);
0106    img = mk_image(imdl,1); vh=fwd_solve(img);
0107    img.elem_data(1)= 1.001;vi=fwd_solve(img);
0108    sig = vi.meas - vh.meas;
0109    [wh,wi] = simulate_movement(img,[0.05;0.05;0.01]);
0110    sigw= wi - wh;
0111    ratio = sigw\sig; 
0112    unit_test_cmp('Ratio ', ratio,.0995,1e-5);
0113    unit_test_cmp('Ratio ', sig*ratio, sigw,1e-3);