simulate_2d_movement

 SIMULATE_2D_MOVEMENT simulate rotational movement in 2D
 [vh,vi,xyr_pt]= simulate_2d_movement( n_points, model, rad_pr, movefcn )

 the target starts at (rad_pr(1),0) and rotates around 
  clockwise
 
   rad_pr = [path_radius, target_radius] = [2/3, .05] (default)
 
   n_points = number of points to simulate (default = 200)

   model = fwd_model to simulate 
         (default use internal, or if model= []);

   movefcn = 1 (Default)  radial motion where the target starts
     at (rad_pr(1),0) and rotates clockwise
   movefcn = 2 movement from centre to outward on to
     at (rad_pr(1),0) on x-axis

   For these movefcn values, the start and end of the complete
      movement may be specified as [movefcn, mv_start, mv_end].
      default values are [movefcn, 0, 1] (ie. the complete movement)

   movefcn = FUCN NAME or FUNC HANDLE
      the function must accept the following parameters
      [xp,yp] = movefcn(f_frac, radius);

 OUTPUT:
   vh - homogeneous measurements M x 1
   vi - target simulations       M x n_points
   xyr_pt - x y and radius of each point 3 x n_points

 For small targets it is more accurate and much faster to
  use the function: simulate_movement.m

0001 function [vh,vi,xyr_pt]= simulate_2d_movement( n_sims, fmdl, rad_pr, movefcn )
0002 % SIMULATE_2D_MOVEMENT simulate rotational movement in 2D
0003 % [vh,vi,xyr_pt]= simulate_2d_movement( n_points, model, rad_pr, movefcn )
0004 %
0005 % the target starts at (rad_pr(1),0) and rotates around
0006 %  clockwise
0007 %
0008 %   rad_pr = [path_radius, target_radius] = [2/3, .05] (default)
0009 %
0010 %   n_points = number of points to simulate (default = 200)
0011 %
0012 %   model = fwd_model to simulate
0013 %         (default use internal, or if model= []);
0014 %
0015 %   movefcn = 1 (Default)  radial motion where the target starts
0016 %     at (rad_pr(1),0) and rotates clockwise
0017 %   movefcn = 2 movement from centre to outward on to
0018 %     at (rad_pr(1),0) on x-axis
0019 %
0020 %   For these movefcn values, the start and end of the complete
0021 %      movement may be specified as [movefcn, mv_start, mv_end].
0022 %      default values are [movefcn, 0, 1] (ie. the complete movement)
0023 %
0024 %   movefcn = FUCN NAME or FUNC HANDLE
0025 %      the function must accept the following parameters
0026 %      [xp,yp] = movefcn(f_frac, radius);
0027 %
0028 % OUTPUT:
0029 %   vh - homogeneous measurements M x 1
0030 %   vi - target simulations       M x n_points
0031 %   xyr_pt - x y and radius of each point 3 x n_points
0032 %
0033 % For small targets it is more accurate and much faster to
0034 %  use the function: simulate_movement.m
0035 
0036 % (C) 2005-2009 Andy Adler. Licensed under GPL v2 or v3
0037 % $Id: simulate_2d_movement.m 5112 2015-06-14 13:00:41Z aadler $
0038 
0039 if nargin>=1 && ischar(n_sims) && strcmp(n_sims,'UNIT_TEST'); do_unit_test; return; end
0040 
0041 if nargin <1
0042    n_sims = 200;
0043 end
0044 
0045 if nargin<2 || isempty(fmdl) % create our own fmdl
0046    n_circles = 36;
0047    n_elec= 16;
0048    fmdl= mk_fwd_model(n_circles, n_elec);
0049 end
0050 
0051 if nargin<3; rad_pr= []; end
0052 if nargin<4; movefcn= 1; end
0053 
0054 copt.cache_obj = {n_sims,fmdl, rad_pr, movefcn};
0055 copt.fstr = 'simulate_2d_movement';
0056 [vh,vi,xyr_pt]= eidors_cache(@do_simulate_2d_movement,{n_sims, fmdl, rad_pr, movefcn},copt );
0057 
0058 
0059 function [vh,vi,xyr_pt]= do_simulate_2d_movement( n_sims, fmdl, rad_pr, movefcn )
0060 if isempty(rad_pr)
0061    radius= 2/3;
0062    rp= .05;
0063 else
0064    radius= rad_pr(1);
0065    rp=     rad_pr(2);
0066 end
0067 
0068     mv_start = 0;
0069     mv_end   = 1;
0070 if isnumeric(movefcn)
0071    if length(movefcn)>=2; mv_start = movefcn(2); end
0072    if length(movefcn)>=3; mv_end   = movefcn(3); end
0073    if     movefcn(1)==1
0074       movefcn = @rotation_path;
0075    elseif movefcn(1)==2
0076       movefcn = @straight_out;
0077    else
0078       error('value of movefcn not understood');
0079    end
0080 else
0081    % assume movefcn is a function
0082 end
0083 
0084     n_elems= size(fmdl.elems,1);
0085     img= eidors_obj('image','simulate_movement', ...
0086                     'fwd_model', fmdl, ...
0087                     'elem_data', ones(n_elems,1) );
0088     vh= fwd_solve(img);
0089 
0090 if 0
0091     np= 256;
0092     maxxy= max(fmdl.nodes);
0093     minxy= min(fmdl.nodes);
0094     [x,y]=meshgrid( linspace(minxy(1),maxxy(1),np), ...
0095                     linspace(minxy(2),maxxy(2),np) );
0096     [eptr,vol]= img_mapper2(fmdl.nodes', fmdl.elems', np, np);
0097 else
0098 
0099     mdl_pts = interp_mesh( fmdl, 8); % 45 per elem
0100     x= squeeze( mdl_pts(:,1,:) );
0101     y= squeeze( mdl_pts(:,2,:) );
0102 end
0103 
0104     % there is a faster way to do this with sparse, but it is confusing
0105 %   eptr_n= zeros(n_elems,1);
0106 %   for i=1:n_elems; eptr_n(i) = sum( eptr(:)==i ); end
0107 %   eptr_n= sparse( eptr(:)+1,1,1, n_elems+1, 1);
0108 %   eptr_n= full(eptr_n(2:end));
0109     
0110     target_conductivity= .1;
0111  
0112     for i=1:n_sims
0113        f_frac= mv_start + ( (i-1)/n_sims ) * (mv_end - mv_start);
0114        fprintf('simulating %d / %d (f_frac=%0.2f) \n',i,n_sims, f_frac);
0115 
0116       [xp,yp]= feval(movefcn, f_frac, radius);
0117 
0118 if 0
0119        xyr_pt(:,i)= [xp;-yp;rp]; % -y because images and axes are reversed
0120        ff= find( (x(:)-xp).^2 + (y(:)-yp).^2 <= rp^2 )';
0121        obj_n= sparse( eptr(ff)+1,1,1, n_elems+1, 1);
0122        obj_n= full(obj_n(2:end));
0123 %      img.elem_data= 1 + target_conductivity * (obj_n./eptr_n);
0124        img.elem_data= 1 + target_conductivity * (obj_n./vol);
0125 else
0126        xyr_pt(:,i)= [xp;yp;rp];
0127        ff=  (x-xp).^2 + (y-yp).^2 <= rp^2;
0128        img.elem_data= 1 + target_conductivity * mean(ff,2);
0129 end
0130 
0131        vi(i)= fwd_solve( img );
0132 %show_fem(img);drawnow; keyboard
0133     end
0134 
0135 % convert to data matrix
0136 vi= [vi(:).meas]; 
0137 vh= [vh(:).meas];
0138 
0139 %   movefcn = 1 (Default)  rotational motion where the target starts
0140 %     at (rad_pr(1),0) and rotates clockwise
0141 % calculate x,y position of point, given f_frac of path
0142 function [xp,yp] = rotation_path(f_frac, radius);
0143    xp= radius * cos(f_frac*2*pi);
0144    yp= radius * sin(f_frac*2*pi);
0145 
0146 function [xp,yp] = straight_out(f_frac, radius);
0147    xp= radius*f_frac;
0148    yp= 0;
0149 
0150 % THis is the code copied from calc_slices
0151 % Search through each element and find the points which
0152 % are in that element
0153 function [EPTR,VOL]= img_mapper2(NODE, ELEM, npx, npy );
0154   xmin = min(NODE(1,:));    xmax = max(NODE(1,:));
0155   xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0156 
0157   ymin = min(NODE(2,:));    ymax = max(NODE(2,:));
0158   ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0159 
0160   range= max([xrange, yrange]);
0161   [x y]=meshgrid( ...
0162       linspace( xmean - range*0.50, xmean + range*0.50, npx ), ...
0163       linspace( ymean + range*0.50, ymean - range*0.50, npy ) );
0164   v_yx= [-y(:) x(:)];
0165   turn= [0 -1 1;1 0 -1;-1 1 0];
0166   EPTR=zeros(npy,npx);
0167   % for each element j, we get points on the simplex a,b,c
0168   %   area A = abc
0169   %   for each candidate point d,
0170   %      area AA = abd + acd + bcd
0171   %      d is in j if AA = A
0172   e= size(ELEM,2);
0173   VOL= zeros(e,1);
0174   for j= 1:e
0175     % calculate area of three subtrianges to each candidate point.
0176     xy= NODE(:,ELEM(:,j))';
0177     % come up with a limited set of candidate points which
0178     % may be within the simplex
0179     endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0180              & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0181     % a is determinant of matrix [i,j,k, xy]
0182     a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0183     VOL(j)= abs(sum(a));
0184 
0185     aa= sum(abs(ones(length(endr),1)*a'+ ...
0186                 v_yx(endr,:)*xy'*turn)');
0187     endr( abs( ( VOL(j)-aa ) ./ VOL(j) ) >1e-8)=[];
0188     EPTR(endr)= j;
0189   end %for j=1:ELEM
0190 
0191 
0192 function mdl_2d= mk_fwd_model(n_circles, n_elec)
0193     %TODO: is there a mk_common_model call that can replace all that?
0194     params= mk_circ_tank(n_circles, [], n_elec); 
0195     n_rings= 1;
0196     options= {'no_meas_current','no_rotate_meas','do_redundant'};
0197     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0198     params.stimulation= st;
0199     params.meas_select= els;
0200     params.solve=      'eidors_default';
0201     params.system_mat= 'eidors_default';
0202     params.jacobian=   'eidors_default';
0203     params.normalize_measurements= 0;
0204     mdl_2d   = eidors_obj('fwd_model', params);
0205 
0206 function do_unit_test
0207    N_TEST = 5;
0208    imdl = mk_common_model( 'c2c2', 16 );
0209    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST);
0210    subplot(421);
0211    imgs = inv_solve(imdl, vh, vi);
0212    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0213 
0214    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [0.3,0.01],1);
0215    subplot(422);
0216    imgs = inv_solve(imdl, vh, vi);
0217    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0218 
0219    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [0.9,0.01],2);
0220    subplot(423);
0221    imgs = inv_solve(imdl, vh, vi);
0222    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0223    
0224    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [],[1,0.5,0.4]);
0225    subplot(424);
0226    imgs = inv_solve(imdl, vh, vi);
0227    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0228 
0229    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, [], [],@test_movefcn);
0230    subplot(425);
0231    imgs = inv_solve(imdl, vh, vi);
0232    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0233 
0234    fmdl = mk_common_model('a2c2',16); fmdl= fmdl.fwd_model;
0235    fmdl.nodes = fmdl.nodes*1.5;
0236    [vh,vi,xyr_pt]=simulate_2d_movement(N_TEST, fmdl, [],@test_movefcn);
0237    subplot(426);
0238    imgs = inv_solve(imdl, vh, vi);
0239    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0240 
0241 function [xp,yp] = test_movefcn(f_frac, radius);
0242   ff = radius/sqrt(2);
0243   xp= ff*f_frac; yp= ff*f_frac;