simulate_3d_movement

 SIMULATE_3D_MOVEMENT simulate rotational movement in 3D
 [vh,vi,xyzr_pt]= simulate_3d_movement( n_points, model, rad_pr, movefcn )

   rad_pr = [path_radius, target_radius, zmin, zmax]
      values are the fraction of the extent in each dimension
      DEFAULT: [2/3, .05, .1, .9 ]
 
   n_points = number of points to simulate (default = 200)

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

   movefcn = 1 (Default)  helical motion where the target starts
     at (rad_pr(1),0) and rotates clockwise moving bottom to top.
   movefcn = 2            radial movement in the vertical plane

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

 OUTPUT:
   vh - homogeneous measurements            M x 1
   vi - target simulations                  M x n_points
   xyzr_pt - x,y,z 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,xyzr_pt]= simulate_3d_movement( n_sims, fmdl, rad_pr,movefcn )
0002 % SIMULATE_3D_MOVEMENT simulate rotational movement in 3D
0003 % [vh,vi,xyzr_pt]= simulate_3d_movement( n_points, model, rad_pr, movefcn )
0004 %
0005 %   rad_pr = [path_radius, target_radius, zmin, zmax]
0006 %      values are the fraction of the extent in each dimension
0007 %      DEFAULT: [2/3, .05, .1, .9 ]
0008 %
0009 %   n_points = number of points to simulate (default = 200)
0010 %
0011 %   model = fwd_model to simulate
0012 %         (default use internal, or if model= []);
0013 %
0014 %   movefcn = 1 (Default)  helical motion where the target starts
0015 %     at (rad_pr(1),0) and rotates clockwise moving bottom to top.
0016 %   movefcn = 2            radial movement in the vertical plane
0017 %
0018 %   movefcn = FUCN NAME or FUNC HANDLE
0019 %      the function must accept the following parameters
0020 %      [xp,yp,zp] = movefcn(f_frac, radius, z_bottom,z_top);
0021 %
0022 % OUTPUT:
0023 %   vh - homogeneous measurements            M x 1
0024 %   vi - target simulations                  M x n_points
0025 %   xyzr_pt - x,y,z and radius of each point 3 x n_points
0026 %
0027 % For small targets it is more accurate and much faster to
0028 %  use the function: simulate_movement.m
0029 
0030 % (C) 2005-2009 Andy Adler. Licensed under GPL v2 or v3
0031 % $Id: simulate_3d_movement.m 5112 2015-06-14 13:00:41Z aadler $
0032 
0033 if nargin <1; n_sims = 200; end
0034 
0035 if nargin>=1 && ischar(n_sims) && strcmp(n_sims,'UNIT_TEST'); do_unit_test; return; end
0036 
0037 if nargin<2 || isempty(fmdl) % create our own fmdl
0038    fmdl= mk_library_model('cylinder_16x2el_fine');
0039    fmdl.normalize_measurements = 0;
0040    fmdl.electrode = fmdl.electrode(1:16);
0041    fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0042    fmdl.normalize_measurements = 0;
0043 end
0044 
0045 if nargin<3 || isempty(rad_pr);
0046    rad_pr= [2/3, 0.05, 0.1, 0.9];
0047 end
0048 
0049 if nargin<4; movefcn= 1; end
0050 
0051 copt.cache_obj = {n_sims,fmdl, rad_pr, movefcn};
0052 copt.fstr = 'simulate_3d_movement';
0053 
0054 [vh,vi,xyzr_pt]= eidors_cache(@do_simulate_3d_movement, ...
0055                               { n_sims, fmdl, rad_pr, movefcn}, copt );
0056                            
0057 
0058 function [vh,vi,xyzr_pt]= do_simulate_3d_movement( n_sims, mdl_3d, rad_pr, movefcn )
0059 
0060     mv_start = 0;
0061     mv_end   = 1;
0062 if isnumeric(movefcn)
0063    if length(movefcn)>=2; mv_start = movefcn(2); end
0064    if length(movefcn)>=3; mv_end   = movefcn(3); end
0065    if     movefcn(1)==1
0066       movefcn = @helical_path;
0067    elseif movefcn(1)==2
0068       movefcn = @radial_path;
0069    else
0070       error('value of movefcn not understood');
0071    end
0072 else
0073    % assume movefcn is a function
0074 end
0075 
0076 eidors_msg('simulate_3d_movement: step #1: homogeneous simulation',2);
0077 % create homogeneous image + simulate data
0078 sigma= ones( size(mdl_3d.elems,1) ,1);
0079 img= eidors_obj('image', 'homogeneous image', ...
0080     'elem_data', sigma, ...
0081     'fwd_model', mdl_3d );
0082 vh = fwd_solve( img);
0083 
0084 eidors_msg('simulate_3d_movement: step #2: find points',2);
0085 
0086     mdl_pts = interp_mesh( mdl_3d, 2); % 10 per elem
0087     x= mdl_pts(:,1,:);
0088     y= mdl_pts(:,2,:);
0089     z= mdl_pts(:,3,:);
0090    [radius,rp,z0,zt] = calc_point_grid(mdl_3d.nodes', rad_pr);
0091 
0092 target_conductivity= .1;
0093 
0094 for i=1:n_sims
0095    if rem(i, max( floor(i/10), 10))==1;
0096        eidors_msg( ...
0097        'simulate_3d_movement: step #3 (%d of %d): target simulations', ...
0098        i, n_sims, 2); 
0099    end
0100 
0101    f_frac= mv_start + ( (i-1)/n_sims ) * (mv_end - mv_start);
0102    [xp,yp,zp]= feval(movefcn, f_frac, radius, z0,zt);
0103    ff=  (x-xp).^2 + (y-yp).^2 + (z-zp).^2 <= rp^2;
0104    img.elem_data= 1 + target_conductivity * mean(ff,3);
0105 
0106    xyzr_pt(:,i)= [xp;-yp;zp;rp]; % -y because images and axes are reversed
0107    vi(i)= fwd_solve( img );% measurement
0108 end
0109 
0110 vi= [vi(:).meas];
0111 vh= [vh(:).meas];
0112 
0113 %   movefcn = 1 (Default)  helical motion where the target starts
0114 %     at (rad_pr(1),0) and rotates clockwise moving bottom to top.
0115 % calculate x,y,z position of point, given f_frac of path
0116 function [xp,yp,zp]= helical_path(f_frac, radius, z0,zt);
0117    xp= radius * cos(f_frac*2*pi);
0118    yp= radius * sin(f_frac*2*pi);
0119    % object moves from bottow to top
0120    zp = z0 + (zt - z0) * f_frac;
0121 
0122 %   movefcn = 2            radial movement in the vertical plane
0123 function [xp,yp,zp]= radial_path(f_frac, radius, z0,zt);
0124    rp= f_frac*radius; 
0125    cv= 2*pi*f_frac;
0126    xp= rp * cos(cv);
0127    yp= rp * sin(cv);
0128    zp = mean([zt,z0]);
0129 
0130 % modified img_mapper3 from calc_slices.m
0131 % this is like tsearch, but doesn't require
0132 % delaunay triangularization. I also wrote it, so I like it :-)
0133 function [EPTR, VOL] = img_mapper3a(NODE, ELEM, x,y,z );
0134    % calc and see if one point is in one element
0135    elems = size(ELEM,2);
0136 
0137    EPTR= zeros(prod(size(x)),1);
0138    VOL= zeros(elems,1);
0139 
0140    for j= 1: elems
0141        if rem(j,5000)==0; fprintf('mapping %d / %d \n',j,elems); end
0142        xyz= NODE(:,ELEM(:,j))';
0143        min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0144        min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0145        min_z= min(xyz(:,3)); max_z= max(xyz(:,3));
0146 
0147        % Simplex relative volume is det([v2-v1,v3-v1, ...])
0148        VOL(j)= abs(det(xyz'*[-1,1,0,0;-1,0,1,0;-1,0,0,1]'));
0149 
0150        xlmax= x<=max_x; if ~any(xlmax); continue; end
0151        xgmin= x>=min_x; if ~any(xgmin); continue; end
0152        ylmax= y<=max_y; if ~any(ylmax); continue; end
0153        ygmin= y>=min_y; if ~any(ygmin); continue; end
0154        zlmax= z<=max_z; if ~any(zlmax); continue; end
0155        zgmin= z>=min_z; if ~any(zgmin); continue; end
0156        % come up with a limited set of candidate points which
0157        % may be within the simplex
0158        endr=find( xlmax & xgmin & ylmax & ygmin & zlmax & zgmin);
0159        ll=  prod(size(endr));
0160        if isempty(ll);
0161           continue;
0162        end
0163 
0164        nn=  size(ELEM,1); %Simplex vertices
0165        vol=zeros(ll,nn);
0166        for i=1:nn
0167            i1= i; i2= rem(i,nn)+1; i3= rem(i+1,nn)+1;
0168            x1= xyz(i1,1)-x(endr); y1= xyz(i1,2)-y(endr); z1= xyz(i1,3)-z(endr);
0169            x2= xyz(i2,1)-x(endr); y2= xyz(i2,2)-y(endr); z2= xyz(i2,3)-z(endr);
0170            x3= xyz(i3,1)-x(endr); y3= xyz(i3,2)-y(endr); z3= xyz(i3,3)-z(endr);
0171            vol(:,i)= x1.*y2.*z3 - x1.*y3.*z2 - x2.*y1.*z3 + ...
0172                x3.*y1.*z2 + x2.*y3.*z1 - x3.*y2.*z1;
0173        end
0174 
0175        endr( sum(abs(vol),2) - VOL(j) >1e-8 )=[];
0176        EPTR(endr)= j;
0177    end %for j=1:ELEM
0178 
0179 function [radius, rp, zmin, zmax,x,y,z] = ...
0180          calc_point_grid(NODE, rad_pr, npx, npy, npz);
0181 
0182    xmin = min(NODE(1,:));    xmax = max(NODE(1,:));
0183    xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0184 
0185    ymin = min(NODE(2,:));    ymax = max(NODE(2,:));
0186    ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0187 
0188    zmin = min(NODE(3,:));    zmax = max(NODE(3,:));
0189    zmean= mean([zmin,zmax]); zrange= zmax-zmin;
0190 
0191    radius= rad_pr(1)*(xmax-xmin)/2;
0192    rp=     rad_pr(2)*(xmax-xmin)/2;
0193    zmin=   (rad_pr(3)-.5)*zrange + zmean;
0194    zmax=   (rad_pr(4)-.5)*zrange + zmean;
0195 
0196    if nargout<=4; return; end
0197    range= max([xrange, yrange,zrange]);
0198    [x y z]=ndgrid( ...
0199        linspace( xmean - range*0.5, xmean + range*0.5, npx ), ...
0200        linspace( ymean + range*0.5, ymean - range*0.5, npy ),...
0201        linspace( zmean - zrange*0.5, zmean + zrange*0.5, npz ));
0202 
0203 function do_unit_test
0204    N_TEST = 5;
0205    imdl = mk_common_model( 'c2c2', 16 );
0206 
0207    [vh,vi,xyzr_pt]=simulate_3d_movement(N_TEST);
0208    subplot(421);
0209    imgs = inv_solve(imdl, vh, vi);
0210    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0211 
0212    [vh,vi,xyzr_pt]=simulate_3d_movement(N_TEST, [], [0.3,0.01,.1,.9],1);
0213    subplot(422);
0214    imgs = inv_solve(imdl, vh, vi);
0215    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0216 
0217    [vh,vi,xyzr_pt]=simulate_3d_movement(N_TEST, [], [0.9,0.01,.1,.9],2);
0218    subplot(423);
0219    imgs = inv_solve(imdl, vh, vi);
0220    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0221    
0222    [vh,vi,xyzr_pt]=simulate_3d_movement(N_TEST, [], [],[1,0.5,0.4]);
0223    subplot(424);
0224    imgs = inv_solve(imdl, vh, vi);
0225    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0226 
0227    [vh,vi,xyzr_pt]=simulate_3d_movement(N_TEST, [], [],@test_movefcn);
0228    subplot(425);
0229    imgs = inv_solve(imdl, vh, vi);
0230    imgs.show_slices.img_cols = N_TEST; show_slices(imgs);
0231 
0232    fmdl = mk_common_model('n3r2',[16,2]); fmdl= fmdl.fwd_model;
0233    fmdl.electrode = fmdl.electrode(1:16);
0234    fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0235    fmdl.nodes = fmdl.nodes*1.5;
0236    [vh,vi,xyzr_pt]=simulate_3d_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,zp] = test_movefcn(f_frac, radius,z0,zt);
0242   ff =  radius/sqrt(2);
0243   xp= ff*f_frac; yp= ff*f_frac; zp = mean([z0,zt]);