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

   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 %   movefcn = FUCN NAME or FUNC HANDLE
0021 %      the function must accept the following parameters
0022 %      [xp,yp] = movefcn(f_frac, radius);
0023 %
0024 % OUTPUT:
0025 %   vh - homogeneous measurements M x 1
0026 %   vi - target simulations       M x n_points
0027 %   xyr_pt - x y and radius of each point 3 x n_points
0028 %
0029 % For small targets it is more accurate and much faster to
0030 %  use the function: simulate_movement.m
0031 
0032 % (C) 2005-2009 Andy Adler. Licensed under GPL v2 or v3
0033 % $Id: simulate_2d_movement.html 2819 2011-09-07 16:43:11Z aadler $
0034 
0035 if nargin <1
0036    n_sims = 200;
0037 end
0038 
0039 if nargin<2 || isempty(fmdl) % create our own fmdl
0040    n_circles = 36;
0041    n_elec= 16;
0042    fmdl= mk_fwd_model(n_circles, n_elec);
0043 end
0044 
0045 if nargin<3 || isempty(rad_pr)
0046    radius= 2/3;
0047    rp= .05;
0048 else
0049    radius= rad_pr(1);
0050    rp=     rad_pr(2);
0051 end
0052 
0053 if nargin<4
0054    movefcn = 1;
0055 end
0056 if isnumeric(movefcn)
0057    if     movefcn==1
0058       movefcn = @rotation_path;
0059    elseif movefcn==2
0060       movefcn = @straight_out;
0061    else
0062       error('value of movefcn not understood');
0063    end
0064 else
0065    % assume movefcn is a function
0066 end
0067 
0068     n_elems= size(fmdl.elems,1);
0069     img= eidors_obj('image','simulate_movement', ...
0070                     'fwd_model', fmdl, ...
0071                     'elem_data', ones(n_elems,1) );
0072     vh= fwd_solve(img);
0073 
0074 if 0
0075     np= 256;
0076     maxxy= max(fmdl.nodes);
0077     minxy= min(fmdl.nodes);
0078     [x,y]=meshgrid( linspace(minxy(1),maxxy(1),np), ...
0079                     linspace(minxy(2),maxxy(2),np) );
0080     [eptr,vol]= img_mapper2(fmdl.nodes', fmdl.elems', np, np);
0081 else
0082 
0083     mdl_pts = interp_mesh( fmdl, 8); % 45 per elem
0084     x= squeeze( mdl_pts(:,1,:) );
0085     y= squeeze( mdl_pts(:,2,:) );
0086 end
0087 
0088     % there is a faster way to do this with sparse, but it is confusing
0089 %   eptr_n= zeros(n_elems,1);
0090 %   for i=1:n_elems; eptr_n(i) = sum( eptr(:)==i ); end
0091 %   eptr_n= sparse( eptr(:)+1,1,1, n_elems+1, 1);
0092 %   eptr_n= full(eptr_n(2:end));
0093     
0094     target_conductivity= .1;
0095  
0096     for i=1:n_sims
0097        f_frac= (i-1)/n_sims;
0098        fprintf('simulating %d / %d \n',i,n_sims);
0099 
0100       [xp,yp]= feval(movefcn, f_frac, radius);
0101 
0102 if 0
0103        xyr_pt(:,i)= [xp;-yp;rp]; % -y because images and axes are reversed
0104        ff= find( (x(:)-xp).^2 + (y(:)-yp).^2 <= rp^2 )';
0105        obj_n= sparse( eptr(ff)+1,1,1, n_elems+1, 1);
0106        obj_n= full(obj_n(2:end));
0107 %      img.elem_data= 1 + target_conductivity * (obj_n./eptr_n);
0108        img.elem_data= 1 + target_conductivity * (obj_n./vol);
0109 else
0110        xyr_pt(:,i)= [xp;yp;rp];
0111        ff=  (x-xp).^2 + (y-yp).^2 <= rp^2;
0112        img.elem_data= 1 + target_conductivity * mean(ff,2);
0113 end
0114 
0115        vi(i)= fwd_solve( img );
0116 %show_fem(img);drawnow; keyboard
0117     end
0118 
0119 % convert to data matrix
0120 vi= [vi(:).meas]; 
0121 vh= [vh(:).meas];
0122 
0123 %   movefcn = 1 (Default)  rotational motion where the target starts
0124 %     at (rad_pr(1),0) and rotates clockwise
0125 % calculate x,y position of point, given f_frac of path
0126 function [xp,yp] = rotation_path(f_frac, radius);
0127    xp= radius * cos(f_frac*2*pi);
0128    yp= radius * sin(f_frac*2*pi);
0129 
0130 function [xp,yp] = straight_out(f_frac, radius);
0131    xp= radius*f_frac;
0132    yp= 0;
0133 
0134 % THis is the code copied from calc_slices
0135 % Search through each element and find the points which
0136 % are in that element
0137 function [EPTR,VOL]= img_mapper2(NODE, ELEM, npx, npy );
0138   xmin = min(NODE(1,:));    xmax = max(NODE(1,:));
0139   xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0140 
0141   ymin = min(NODE(2,:));    ymax = max(NODE(2,:));
0142   ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0143 
0144   range= max([xrange, yrange]);
0145   [x y]=meshgrid( ...
0146       linspace( xmean - range*0.50, xmean + range*0.50, npx ), ...
0147       linspace( ymean + range*0.50, ymean - range*0.50, npy ) );
0148   v_yx= [-y(:) x(:)];
0149   turn= [0 -1 1;1 0 -1;-1 1 0];
0150   EPTR=zeros(npy,npx);
0151   % for each element j, we get points on the simplex a,b,c
0152   %   area A = abc
0153   %   for each candidate point d,
0154   %      area AA = abd + acd + bcd
0155   %      d is in j if AA = A
0156   e= size(ELEM,2);
0157   VOL= zeros(e,1);
0158   for j= 1:e
0159     % calculate area of three subtrianges to each candidate point.
0160     xy= NODE(:,ELEM(:,j))';
0161     % come up with a limited set of candidate points which
0162     % may be within the simplex
0163     endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0164              & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0165     % a is determinant of matrix [i,j,k, xy]
0166     a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0167     VOL(j)= abs(sum(a));
0168 
0169     aa= sum(abs(ones(length(endr),1)*a'+ ...
0170                 v_yx(endr,:)*xy'*turn)');
0171     endr( abs( ( VOL(j)-aa ) ./ VOL(j) ) >1e-8)=[];
0172     EPTR(endr)= j;
0173   end %for j=1:ELEM
0174 
0175 
0176 function mdl_2d= mk_fwd_model(n_circles, n_elec)
0177     params= mk_circ_tank(n_circles, [], n_elec); 
0178     n_rings= 1;
0179     options= {'no_meas_current','no_rotate_meas','do_redundant'};
0180     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0181     params.stimulation= st;
0182     params.meas_select= els;
0183     params.solve=      'aa_fwd_solve';
0184     params.system_mat= 'aa_calc_system_mat';
0185     params.jacobian=   'aa_calc_jacobian';
0186     params.normalize_measurements= 1;
0187     params.np_fwd_solve.perm_sym= '{n}';
0188     mdl_2d   = eidors_obj('fwd_model', params);