simulate_3d_movement

 SIMULATE_3D_MOVEMENT simulate rotational movement in 2D
 [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, mdl_3d, rad_pr,movefcn )
0002 % SIMULATE_3D_MOVEMENT simulate rotational movement in 2D
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.html 2819 2011-09-07 16:43:11Z aadler $
0032 
0033 if nargin <1
0034    n_sims = 200;
0035 end
0036 
0037 if nargin<2 || isempty(mdl_3d) % create our own fmdl
0038    mdl_3d= mk_common_model('n3r2',[16,2]); % NP's demo model
0039    mdl_3d= mdl_3d.fwd_model;
0040 end
0041 
0042 if nargin<3 || isempty(rad_pr);
0043    rad_pr= [2/3, 0.05, 0.1, 0.9];
0044 end
0045 
0046 if nargin<4
0047    movefcn = 1;
0048 end
0049 if isnumeric(movefcn)
0050    if     movefcn==1
0051       movefcn = @helical_path;
0052    elseif movefcn==2
0053       movefcn = @radial_path;
0054    else
0055       error('value of movefcn not understood');
0056    end
0057 else
0058    % assume movefcn is a function
0059 end
0060 
0061 
0062 if 0 % OLD CODE - now use solver provided
0063    mdl_3d.solve=      'np_fwd_solve';
0064    mdl_3d.system_mat= 'np_calc_system_mat';
0065    mdl_3d.jacobian=   'np_calc_jacobian';
0066    mdl_3d.misc.perm_sym=   '{n}';
0067    n_elems= size(mdl_3d.elems,1);
0068 end
0069 
0070 
0071 eidors_msg('simulate_3d_movement: step #1: homogeneous simulation',2);
0072 % create homogeneous image + simulate data
0073 sigma= ones( size(mdl_3d.elems,1) ,1);
0074 img= eidors_obj('image', 'homogeneous image', ...
0075     'elem_data', sigma, ...
0076     'fwd_model', mdl_3d );
0077 vh = fwd_solve( img);
0078 
0079 eidors_msg('simulate_3d_movement: step #2: find points',2);
0080 
0081 if 0 % Old Style
0082    npx=128;
0083    npy=128;
0084    npz=64;
0085    [radius,rp,z0,zt,x,y,z] = calc_point_grid(mdl_3d.nodes', rad_pr, npx, npy, npz);
0086 
0087    clear pts;
0088    for i=1:n_sims
0089        f_frac= (i-1)/n_sims;
0090 
0091        % call function to simulate data
0092        [xp,yp,zp]= feval(movefcn, f_frac, radius, z0,zt);
0093 
0094        ff= find( (x(:)-xp).^2 + (y(:)-yp).^2 + (z(:)-zp).^2 <= rp^2 )';
0095        pts{i} = ff;
0096    end
0097    pts_all = unique( [pts{:}] );
0098    pts_all = pts_all(:);
0099    for i=1:n_sims
0100       [jnk,idx_i]= intersect( pts_all, pts{i});
0101       pts_idx{i}= idx_i;
0102    end
0103 
0104    [eptr,vol]= img_mapper3a(mdl_3d.nodes', mdl_3d.elems',  ...
0105             x(pts_all), y(pts_all), z(pts_all));
0106 else
0107     mdl_pts = interp_mesh( mdl_3d, 2); % 10 per elem
0108     x= mdl_pts(:,1,:);
0109     y= mdl_pts(:,2,:);
0110     z= mdl_pts(:,3,:);
0111    [radius,rp,z0,zt] = calc_point_grid(mdl_3d.nodes', rad_pr);
0112 end
0113 
0114 target_conductivity= .1;
0115 
0116 for i=1:n_sims
0117    if rem(i, max( floor(i/10), 10))==1; eidors_msg( ...
0118        'simulate_3d_movement: step #3 (%d of %d): target simulations', ...
0119        i, n_sims, 2); 
0120    end
0121    if 0
0122       obj_n= sparse( eptr(pts_idx{i}),1,1, n_elems, 1);
0123       img.elem_data= 1 + target_conductivity * obj_n./vol;
0124   %   show_fem(img); view([-2,84]);pause;
0125    else
0126       f_frac= (i-1)/n_sims;
0127       [xp,yp,zp]= feval(movefcn, f_frac, radius, z0,zt);
0128       ff=  (x-xp).^2 + (y-yp).^2 + (z-zp).^2 <= rp^2;
0129       img.elem_data= 1 + target_conductivity * mean(ff,3);
0130    end
0131 
0132    xyzr_pt(:,i)= [xp;-yp;zp;rp]; % -y because images and axes are reversed
0133    vi(i)= fwd_solve( img );% measurement
0134 end
0135 
0136 vi= [vi(:).meas];
0137 vh= [vh(:).meas];
0138 
0139 %   movefcn = 1 (Default)  helical motion where the target starts
0140 %     at (rad_pr(1),0) and rotates clockwise moving bottom to top.
0141 % calculate x,y,z position of point, given f_frac of path
0142 function [xp,yp,zp]= helical_path(f_frac, radius, z0,zt);
0143    xp= radius * cos(f_frac*2*pi);
0144    yp= radius * sin(f_frac*2*pi);
0145    % object moves from bottow to top
0146    zp = z0 + (zt - z0) * f_frac;
0147 
0148 %   movefcn = 2            radial movement in the vertical plane
0149 function [xp,yp,zp]= radial_path(f_frac, radius, z0,zt);
0150    rp= f_frac*radius; 
0151    cv= 2*pi*f_frac;
0152    xp= rp * cos(cv);
0153    yp= rp * sin(cv);
0154    zp = mean([zt,z0]);
0155 
0156 % modified img_mapper3 from calc_slices.m
0157 % this is like tsearch, but doesn't require
0158 % delaunay triangularization. I also wrote it, so I like it :-)
0159 function [EPTR, VOL] = img_mapper3a(NODE, ELEM, x,y,z );
0160    % calc and see if one point is in one element
0161    elems = size(ELEM,2);
0162 
0163    EPTR= zeros(prod(size(x)),1);
0164    VOL= zeros(elems,1);
0165 
0166    for j= 1: elems
0167        if rem(j,5000)==0; fprintf('mapping %d / %d \n',j,elems); end
0168        xyz= NODE(:,ELEM(:,j))';
0169        min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0170        min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0171        min_z= min(xyz(:,3)); max_z= max(xyz(:,3));
0172 
0173        % Simplex relative volume is det([v2-v1,v3-v1, ...])
0174        VOL(j)= abs(det(xyz'*[-1,1,0,0;-1,0,1,0;-1,0,0,1]'));
0175 
0176        xlmax= x<=max_x; if ~any(xlmax); continue; end
0177        xgmin= x>=min_x; if ~any(xgmin); continue; end
0178        ylmax= y<=max_y; if ~any(ylmax); continue; end
0179        ygmin= y>=min_y; if ~any(ygmin); continue; end
0180        zlmax= z<=max_z; if ~any(zlmax); continue; end
0181        zgmin= z>=min_z; if ~any(zgmin); continue; end
0182        % come up with a limited set of candidate points which
0183        % may be within the simplex
0184        endr=find( xlmax & xgmin & ylmax & ygmin & zlmax & zgmin);
0185        ll=  prod(size(endr));
0186        if isempty(ll);
0187           continue;
0188        end
0189 
0190        nn=  size(ELEM,1); %Simplex vertices
0191        vol=zeros(ll,nn);
0192        for i=1:nn
0193            i1= i; i2= rem(i,nn)+1; i3= rem(i+1,nn)+1;
0194            x1= xyz(i1,1)-x(endr); y1= xyz(i1,2)-y(endr); z1= xyz(i1,3)-z(endr);
0195            x2= xyz(i2,1)-x(endr); y2= xyz(i2,2)-y(endr); z2= xyz(i2,3)-z(endr);
0196            x3= xyz(i3,1)-x(endr); y3= xyz(i3,2)-y(endr); z3= xyz(i3,3)-z(endr);
0197            vol(:,i)= x1.*y2.*z3 - x1.*y3.*z2 - x2.*y1.*z3 + ...
0198                x3.*y1.*z2 + x2.*y3.*z1 - x3.*y2.*z1;
0199        end
0200 
0201        endr( sum(abs(vol),2) - VOL(j) >1e-8 )=[];
0202        EPTR(endr)= j;
0203    end %for j=1:ELEM
0204 
0205 function [radius, rp, zmin, zmax,x,y,z] = ...
0206          calc_point_grid(NODE, rad_pr, npx, npy, npz);
0207 
0208    xmin = min(NODE(1,:));    xmax = max(NODE(1,:));
0209    xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0210 
0211    ymin = min(NODE(2,:));    ymax = max(NODE(2,:));
0212    ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0213 
0214    zmin = min(NODE(3,:));    zmax = max(NODE(3,:));
0215    zmean= mean([zmin,zmax]); zrange= zmax-zmin;
0216 
0217    radius= rad_pr(1)*(xmax-xmin)/2;
0218    rp=     rad_pr(2)*(xmax-xmin)/2;
0219    zmin=   (rad_pr(3)-.5)*zrange + zmean;
0220    zmax=   (rad_pr(4)-.5)*zrange + zmean;
0221 
0222    if nargout<=4; return; end
0223    range= max([xrange, yrange,zrange]);
0224    [x y z]=ndgrid( ...
0225        linspace( xmean - range*0.5, xmean + range*0.5, npx ), ...
0226        linspace( ymean + range*0.5, ymean - range*0.5, npy ),...
0227        linspace( zmean - zrange*0.5, zmean + zrange*0.5, npz ));