ng_mk_closed_contour

PURPOSE ^

NG_MK_CLOSED_CONTOUR: fit elliptical model to a contour

SYNOPSIS ^

function [fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos,elec_shape, extra_ng_code);

DESCRIPTION ^

 NG_MK_CLOSED_CONTOUR: fit elliptical model to a contour
[fmdl,mat_idx] = ng_mk_closed_contour(shape, elec_pos, ...
                 elec_shape, extra_ng_code);

 This function creates netgen models of a closed contour in x,y.
   An elliptical FEM is created and then perturbed using Fourier 
   descriptors into the final model.
 It should work well for cases where the contour is roughtly elliptical

 INPUT:
 ellip_shape = cell{height, xy_points, [maxsz]]}
    if height = 0 -> calculate a 2D shape
    xy_points     -> Npoints x 2. points on the object boundary
    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)

 ELECTRODE POSITIONS:
  elec_pos = [n_elecs_per_plane,z_planes] 
     OR
  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
 Note that electrode positions are in the ellipse before fitting

 ELECTRODE SHAPES::
  elec_shape = [width,height, maxsz]  % Rectangular elecs
     OR
  elec_shape = [radius, 0, maxsz ]    % Circular elecs
     OR 
  elec_shape = [0, 0, maxsz ]         % Point electrodes
    (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask)

 Specify either a common electrode shape or for each electrode

 EXTRA_NG_CODE
   string of extra code to put into netgen geo file. Normally this
   would be to insert extra materials into the space

 OUTPUT:
  fmdl    - fwd_model object
  mat_idx - indices of materials (if extra_ng_code is used)
    Note mat_idx does not work in 2D. Netgen does not provide it.


 USAGE EXAMPLES:
 Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0002                   elec_shape, extra_ng_code);
0003 % NG_MK_CLOSED_CONTOUR: fit elliptical model to a contour
0004 %[fmdl,mat_idx] = ng_mk_closed_contour(shape, elec_pos, ...
0005 %                 elec_shape, extra_ng_code);
0006 %
0007 % This function creates netgen models of a closed contour in x,y.
0008 %   An elliptical FEM is created and then perturbed using Fourier
0009 %   descriptors into the final model.
0010 % It should work well for cases where the contour is roughtly elliptical
0011 %
0012 % INPUT:
0013 % ellip_shape = cell{height, xy_points, [maxsz]]}
0014 %    if height = 0 -> calculate a 2D shape
0015 %    xy_points     -> Npoints x 2. points on the object boundary
0016 %    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)
0017 %
0018 % ELECTRODE POSITIONS:
0019 %  elec_pos = [n_elecs_per_plane,z_planes]
0020 %     OR
0021 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0022 % Note that electrode positions are in the ellipse before fitting
0023 %
0024 % ELECTRODE SHAPES::
0025 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0026 %     OR
0027 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0028 %     OR
0029 %  elec_shape = [0, 0, maxsz ]         % Point electrodes
0030 %    (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask)
0031 %
0032 % Specify either a common electrode shape or for each electrode
0033 %
0034 % EXTRA_NG_CODE
0035 %   string of extra code to put into netgen geo file. Normally this
0036 %   would be to insert extra materials into the space
0037 %
0038 % OUTPUT:
0039 %  fmdl    - fwd_model object
0040 %  mat_idx - indices of materials (if extra_ng_code is used)
0041 %    Note mat_idx does not work in 2D. Netgen does not provide it.
0042 %
0043 %
0044 % USAGE EXAMPLES:
0045 % Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
0046 %     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
0047 %
0048 
0049 % (C) Andy Adler, 2010. Licenced under GPL v2 or v3
0050 % $Id: ng_mk_closed_contour.m 2267 2010-09-23 11:04:24Z aadler $
0051 
0052 if nargin < 4; extra_ng_code = {'',''}; end
0053 % Caching should be done in the ng_mk_ellip_models
0054 
0055 [ellip_shape,params] = analyse_shape(ellip_shape) 
0056 [fmdl, mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0057                   elec_shape, extra_ng_code);
0058 fmdl = fit_to_shape( fmdl, params );
0059 
0060 function [ellip_shape,params] = analyse_shape(ellip_shape);
0061   height = ellip_shape{1}; 
0062   xy_shape=ellip_shape{2};
0063   if length(ellip_shape)>=3
0064     maxh= ellip_shape{3};
0065   else
0066     maxh= [];
0067   end
0068 
0069   [u,d,v]= svd(cov(xy_shape));
0070   params.rot = [1,0;0,-1]*u;
0071   params.xy_mean = mean(xy_shape,1);
0072 
0073   ellip_ax = sqrt(2*[d(1,1),d(2,2)])  
0074   ellip_shape= [height, ellip_ax, maxh];
0075 
0076   N = min(15, ceil( size(xy_shape,1) / 3 )); % can't fit too many components
0077   params.f_fit_point = fourier_fit(xy_shape, N);
0078 
0079 % Now, fit the ellipse to the point
0080   th = linspace(0,2*pi,4*N)';
0081   ellip_pts = [cos(th),sin(th)]*sqrt(2*d)*params.rot;
0082   params.f_fit_ellip = fourier_fit(ellip_pts, N);
0083 %plot(ellip_pts(:,1),ellip_pts(:,2),'b-*',xy_shape(:,1),xy_shape(:,2),'r-*')
0084 
0085 function C = fourier_fit(xy_shape,N);
0086 % Fourier Fit
0087 % [1, cos(t1), cos(2*t1), ... sin(t1) ...] * [C1] = [R1]
0088 %            ...
0089 % [1, cos(tN), cos(2*tN), ... sin(tN) ...] * [CN] = [RM]
0090 
0091    xc = xy_shape(:,1); xm= mean(xc);   xc= xc-xm;
0092    yc = xy_shape(:,2); ym= mean(yc);   yc= yc-ym;
0093 
0094    ang = atan2(yc,xc);
0095    rad = sqrt(yc.^2 + xc.^2);
0096    A = ones(length(ang), 2*N+1);
0097    for i=1:N
0098      A(:,i+1  ) = cos(i*ang);
0099      A(:,i+1+N) = sin(i*ang);
0100    end
0101    C = A\rad;
0102 
0103 function fmdl = fit_to_shape( fmdl, params );
0104    % Rotate and move the objects
0105    xnodes = fmdl.nodes(:,1:2)*params.rot(:,1);
0106    ynodes = fmdl.nodes(:,1:2)*params.rot(:,2);
0107 
0108    N = (length(params.f_fit_point)-1)/2;
0109    [ang,rad] = cart2pol(xnodes,ynodes);
0110    A = ones(length(ang), 2*N+1);
0111    for i=1:N
0112      A(:,i+1  ) = cos(i*ang);
0113      A(:,i+1+N) = sin(i*ang);
0114    end
0115 
0116    r_ellip = A*params.f_fit_ellip;
0117    r_point = A*params.f_fit_point;
0118    rad = rad.* (r_point ./ r_ellip);
0119    xnodes = rad.*cos(ang) + params.xy_mean(1);
0120    ynodes = rad.*sin(ang) + params.xy_mean(2);
0121 
0122    fmdl.nodes(:,1) = xnodes;
0123    fmdl.nodes(:,2) = ynodes;
0124

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005