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