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 );
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.html 2819 2011-09-07 16:43:11Z aadler $
0038 
0039    if size(xyzr) == [1,1]
0040       path = linspace(0,1,xyzr); phi = 2*pi*path;
0041       meanodes= mean(    img.fwd_model.nodes  );
0042       lennodes= size( img.fwd_model.nodes,1); 
0043       img.fwd_model.nodes = img.fwd_model.nodes - ones(lennodes,1)*meanodes;
0044       maxnodes= max(max(abs( img.fwd_model.nodes(:,1:2) )));
0045       img.fwd_model.nodes = img.fwd_model.nodes / maxnodes;
0046 
0047       xyzr = [0.9*path.*sin(phi);0.9*path.*cos(phi);0*path; 0.05 + 0*path];
0048    end
0049 
0050    Nt = size(xyzr,2);
0051    [c2f failed] = mk_c2f_circ_mapping( img.fwd_model, xyzr);
0052    xyzr(:,failed) = [];
0053    c2f(:,failed) = [];
0054    Nt = Nt - nnz(failed);
0055    img.fwd_model.coarse2fine = c2f;
0056 
0057    % We don't want a normalized jacobian here
0058    img.fwd_model.normalize_measurements = 0;
0059 
0060    J= calc_jacobian( img );
0061    J= move_jacobian_postprocess( J, img, Nt);
0062 
0063    vh= fwd_solve(img);
0064    vh=vh.meas;
0065 
0066    vi= vh*ones(1,Nt) + J;
0067    
0068    % Would this be the slow approach?:
0069 %    vi = vh*zeros(1,Nt);
0070 %    for i = 1: Nt
0071 %        img.elem_data = 1 - c2f(:,i);
0072 %        jnk = fwd_solve(img);
0073 %        vi(:,i) = jnk.meas;
0074 %    end
0075 
0076 function J= move_jacobian_postprocess( J, img, Nt)
0077    if size(J,2) == Nt; % No problem
0078        return ;
0079    % Check if movement jacobian introduced electrode movements (elecs * dims)
0080    elseif size(J,2) == Nt + ...
0081         length(img.fwd_model.electrode) * size(img.fwd_model.nodes,2)
0082       J = J(:,1:Nt);
0083    else
0084       error('Jacobian calculator is not doing the coarse2fine mapping. This is a bug.');
0085    end
0086 
0087 
0088 
0089 
0090 
0091 % QUESTON: the accuracy of J will depend on how well we interpolate the
0092 % mesh. How important is this?