deform_cylinder

 fwd_mdl = deform_cylinder( fwd_mdl, niv);
 Deform the boundary of the cylinder to make it like a torso
 niv= 1.. 5 => Torso shape from T5 - T12
 xyz_expand - rescale xyz - default should be [1];

 This function is called by the '?2t#' (eg c2t2) options
  in mk_common_model. It is kept for compatibility with 
  old tutorials, but we don't recommend it's use for new code.

0001 function fwd_mdl = deform_cylinder( fwd_mdl, geo);
0002 % fwd_mdl = deform_cylinder( fwd_mdl, niv);
0003 % Deform the boundary of the cylinder to make it like a torso
0004 % niv= 1.. 5 => Torso shape from T5 - T12
0005 % xyz_expand - rescale xyz - default should be [1];
0006 %
0007 % This function is called by the '?2t#' (eg c2t2) options
0008 %  in mk_common_model. It is kept for compatibility with
0009 %  old tutorials, but we don't recommend it's use for new code.
0010 
0011 % (C) 1995 Andy Adler. License: GPL version 2 or 3
0012 % $Id: deform_cylinder.m 5112 2015-06-14 13:00:41Z aadler $
0013 
0014 if ischar(fwd_mdl) && strcmp(fwd_mdl,'UNIT_TEST'); do_unit_test; return; end
0015 
0016     NODE= fwd_mdl.nodes';
0017     a_max= size(geo.xy,1);
0018     ab_geo=sqrt(sum(([ geo.xy; geo.xy(1,:) ]').^2)');
0019     nn= zeros(size(NODE));
0020     for i=1:size(NODE,2);
0021       angle = rem(a_max*atan2( NODE(2,i), ...
0022             NODE(1,i) )/2/pi+a_max,a_max)+1;
0023       fl_angl = floor(angle + eps );
0024       fac=(  (fl_angl + 1 - angle) * ab_geo(fl_angl) + ...
0025              (angle - fl_angl)     * ab_geo(fl_angl + 1)  );
0026       nn(1:2,i)= NODE(1:2,i)* fac;
0027     end  %for i=1:size
0028     if size(nn,1) == 3;
0029        nn(3,:) = NODE(3,:)*geo.z_mag;
0030     end
0031 
0032     xyz_expand = 1;
0033     fwd_mdl.nodes = nn'*eye(xyz_expand);
0034 
0035 function fmdl_out = deform_cylinder2( fmdl, deform_pararms)
0036 % fmdl_out = deform_cylinder( fmdl_in, deform_pararms)
0037 % DEFORM_CYLINDER: deform cylinder from circle to
0038 %   nearly circular shape. Need to be a bit careful,
0039 %   because too much deformation will make for poorly
0040 %   shaped FEM elements.
0041 %
0042 % fmdl_out = output (deformed) fmdl
0043 % fmdl_in  = input (non-deformed) fmdl
0044 %
0045 % deform_params values
0046 %    .axes (rotate axes before and after deform)
0047 %       deformation is in x,y axis
0048 %       either as vector [1,3,2] or 2x2 (2D) or 3x3 mat
0049 %       Default = [1,2,3]
0050 %    .mode    = 'add' = add deformation
0051 %               'mul' = multiply by deformation
0052 %    .xycentre  = [xctr, yctr] (DEFAULT = 0,0)
0053 %
0054 %    .xyspecpos = [xelec, yelec]
0055 %         specified new electrode positions
0056 %
0057 %    .zerorad = radius at which deform goes to zero
0058 %        NOT IMPLEMENTED YET
0059 
0060 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test;return;end
0061 
0062 pp = proc_params( deform_pararms, fmdl )
0063 
0064 function pp = proc_params( deform_pararms, fmdl )
0065    d = size(fmdl.nodes,2);
0066    if isfield(deform_pararms,'axes')
0067      dpax = deform_pararms.axes;
0068      [r,c] = size(dpax);
0069      if r==1;
0070         pp.axes = sparse(dpax,1:c,1,c,c); 
0071      else 
0072         pp.axes = sparse(dpax);
0073 
0074      end
0075    else 
0076      pp.axes = speye(d);
0077    end
0078 
0079    if isfield(deform_pararms,'mode')
0080      switch deform_pararms.mode
0081         case 'add'; pp.mode = 1;
0082         case 'mul'; pp.mode = 2;
0083         otherwise; error(['deform mode not understood']);
0084      end
0085    else
0086      pp.mode = 2;
0087    end
0088 
0089    if isfield(deform_params,'xycentre')
0090      pp.xctr = deform_pararms.xycentre(1);
0091      pp.yctr = deform_pararms.xycentre(2);
0092    else
0093      pp.xctr = 0;
0094      pp.yctr = 0;
0095    end
0096 
0097    % deform_pararms.xyspecpos must be given
0098    xy= deform_pararms.xyspecpos;
0099   [pp.th_s,pp.r_s] = cart2pol(xy(:,1)-pp.xctr, ...
0100                               xy(:,2)-pp.yctr);
0101    
0102    % Find rad, theta of electrodes
0103    for i=1:length(fmdl.electrode);
0104       enodes = fmdl.electrode(i).nodes;
0105       fepos(i,:) = mean(fmdl.nodes(enodes,:),1);
0106    end
0107 
0108   [pp.th_e,pp.r_e] = cart2pol(fepos(:,1)-pp.xctr, ...
0109                               fepos(:,2)-pp.yctr);
0110 
0111 function deform(fmdl, pp);
0112 th = mean([the,thm],2);
0113 dr = re-rm;
0114 dr = [dr;dr;dr];
0115 th = [th-2*pi;th;th+2*pi];
0116 [th,idx]=sort(th); dr= dr(idx);
0117 tpf = polyfit(th,idx,10);
0118 plot(th,dr,'*');
0119 
0120 [thn,rn] = cart2pol( fmdl.nodes(:,2)-yc, fmdl.nodes(:,3)-zc ); 
0121 drn = interp1(th,dr,thn);
0122 rn = rn + drn;
0123 [nodesy, nodesz] = pol2cart(thn,rn);
0124 fmdl.nodes(:,2:3) = [nodesy+yc, nodesz+zc];
0125 
0126 hh=plot(elec_pos(:,2),elec_pos(:,3),'*-',fepos(:,2),fepos(:,3),'o-');
0127 set(hh,'LineWidth',2);
0128 hold on;
0129 f1ml = fmdl; f1ml.nodes = f1ml.nodes(:,[2,3,1]);
0130 show_fem(f1ml);
0131 hold off;
0132 
0133 function do_unit_test
0134    fmdl= ng_mk_cyl_models(1,0,[]);
0135    th = linspace(0,2*pi,50)';
0136    geo.xy=[sin(th),2*cos(th)];
0137    geo.z_mag = 10;
0138    show_fem(deform_cylinder(fmdl,geo));
0139    view(2);