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:
Thorax shape
xy_points = .01*[
963 534;1013 538;1056 535;1102 517;1130 492;
1150 454;1169 419;1182 375;1172 345;1149 318;
1112 301;1095 297;1039 285; 982 284; 928 283;
870 279; 816 299; 776 328; 761 374; 782 429;
805 482; 844 518; 900 533; 961 532];
% refine cylendar from ellip_shape
ng_write_opt('MSZCYLINDER',[0,0,1,0,0,1.2,2,.05]);
fmdl= ng_mk_closed_contour({3,xy_points},[8,1.1,1.9],[0.1]);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 % Thorax shape 0046 % xy_points = .01*[ 0047 % 963 534;1013 538;1056 535;1102 517;1130 492; 0048 % 1150 454;1169 419;1182 375;1172 345;1149 318; 0049 % 1112 301;1095 297;1039 285; 982 284; 928 283; 0050 % 870 279; 816 299; 776 328; 761 374; 782 429; 0051 % 805 482; 844 518; 900 533; 961 532]; 0052 % % refine cylendar from ellip_shape 0053 % ng_write_opt('MSZCYLINDER',[0,0,1,0,0,1.2,2,.05]); 0054 % fmdl= ng_mk_closed_contour({3,xy_points},[8,1.1,1.9],[0.1]); 0055 % 0056 0057 % (C) Andy Adler, 2010. Licenced under GPL v2 or v3 0058 % $Id: ng_mk_closed_contour.m 6934 2024-06-12 20:27:34Z aadler $ 0059 0060 if ischar(ellip_shape) && strcmp(ellip_shape,'UNIT_TEST'); do_unit_test; return; end 0061 0062 if nargin < 4; extra_ng_code = {'',''}; end 0063 % Caching should be done in the ng_mk_ellip_models 0064 0065 [ellip_shape,params] = analyse_shape(ellip_shape); 0066 [fmdl, mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ... 0067 elec_shape, extra_ng_code); 0068 fmdl = fit_to_shape( fmdl, params ); 0069 0070 function [ellip_shape,params] = analyse_shape(ellip_shape); 0071 height = ellip_shape{1}; 0072 xy_shape=ellip_shape{2}; 0073 if length(ellip_shape)>=3 0074 maxh= ellip_shape{3}; 0075 else 0076 maxh= []; 0077 end 0078 0079 [u,d,v]= svd(cov(xy_shape)); 0080 params.rot = [1,0;0,-1]*u; 0081 params.xy_mean = mean(xy_shape,1); 0082 0083 ellip_ax = sqrt(2*[d(1,1),d(2,2)]); 0084 ellip_shape= [height, ellip_ax, maxh]; 0085 0086 N = min(15, ceil( size(xy_shape,1) / 3 )); % can't fit too many components 0087 params.f_fit_point = fourier_fit(xy_shape, N); 0088 0089 % Now, fit the ellipse to the point 0090 th = linspace(0,2*pi,4*N)'; 0091 ellip_pts = [cos(th),sin(th)]*sqrt(2*d)*params.rot; 0092 params.f_fit_ellip = fourier_fit(ellip_pts, N); 0093 %plot(ellip_pts(:,1),ellip_pts(:,2),'b-*',xy_shape(:,1),xy_shape(:,2),'r-*') 0094 0095 function C = fourier_fit(xy_shape,N); 0096 % Fourier Fit 0097 % [1, cos(t1), cos(2*t1), ... sin(t1) ...] * [C1] = [R1] 0098 % ... 0099 % [1, cos(tN), cos(2*tN), ... sin(tN) ...] * [CN] = [RM] 0100 0101 xc = xy_shape(:,1); xm= mean(xc); xc= xc-xm; 0102 yc = xy_shape(:,2); ym= mean(yc); yc= yc-ym; 0103 0104 ang = atan2(yc,xc); 0105 rad = sqrt(yc.^2 + xc.^2); 0106 A = ones(length(ang), 2*N+1); 0107 for i=1:N 0108 A(:,i+1 ) = cos(i*ang); 0109 A(:,i+1+N) = sin(i*ang); 0110 end 0111 C = A\rad; 0112 0113 function fmdl = fit_to_shape( fmdl, params ); 0114 % Rotate and move the objects 0115 xnodes = fmdl.nodes(:,1:2)*params.rot(:,1); 0116 ynodes = fmdl.nodes(:,1:2)*params.rot(:,2); 0117 0118 N = (length(params.f_fit_point)-1)/2; 0119 [ang,rad] = cart2pol(xnodes,ynodes); 0120 A = ones(length(ang), 2*N+1); 0121 for i=1:N 0122 A(:,i+1 ) = cos(i*ang); 0123 A(:,i+1+N) = sin(i*ang); 0124 end 0125 0126 r_ellip = A*params.f_fit_ellip; 0127 r_point = A*params.f_fit_point; 0128 rad = rad.* (r_point ./ r_ellip); 0129 xnodes = rad.*cos(ang) + params.xy_mean(1); 0130 ynodes = rad.*sin(ang) + params.xy_mean(2); 0131 0132 fmdl.nodes(:,1) = xnodes; 0133 fmdl.nodes(:,2) = ynodes; 0134 0135 function do_unit_test 0136 xy_points = [ 0137 963 534; 987 535;1013 538;1036 540;1056 535; 0138 1077 527;1102 517;1115 506;1130 492;1141 474; 0139 1150 454;1157 435;1169 419;1179 397;1182 375; 0140 1178 358;1172 345;1163 331;1149 318;1132 307; 0141 1112 301;1097 293;1095 297;1077 289;1039 285; 0142 1011 284; 982 284; 955 285; 928 283; 900 280; 0143 870 279; 845 289; 816 299; 788 310; 776 328; 0144 767 349; 761 374; 768 395; 782 429; 795 456; 0145 805 482; 823 505; 844 518; 868 523; 900 533; 0146 929 536; 961 532]/100; 0147 fmdl= ng_mk_closed_contour({3,xy_points},[8,1.1,1.9],[0.1]); 0148 show_fem(fmdl) 0149 unit_test_cmp('contour A',num_elecs(fmdl),16); 0150 mnod = max(fmdl.nodes,[],1); 0151 unit_test_cmp('contour B',mnod(3),3) 0152 unit_test_cmp('contour C',mnod(1:2),max(xy_points),3e-2) 0153 mnod = min(fmdl.nodes,[],1); 0154 unit_test_cmp('contour D',mnod(3),0) 0155 unit_test_cmp('contour E',mnod(1:2),min(xy_points),3e-2) 0156