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];

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