NG_MAKE_CYL_MODELS: create cylindrical models using netgen
[fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ...
elec_shape, extra_ng_code);
INPUT:
cyl_shape = {height, [radius, [maxsz]]}
if height = 0 -> calculate a 2D shape
radius (OPT) -> (default = 1)
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)
ELECTRODE SHAPES::
elec_shape = [width,height, maxsz] % Rectangular elecs
OR
elec_shape = [radius, 0, maxsz ] % Circular elecs
OR
elec_shape = [0, 0, maxsz ] % Point elecs
(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 cylinder. Radius = 1. No electrodes
fmdl= ng_mk_cyl_models(3,[0],[]);
Simple 2D cylinder. Radius = 2. Set minsize to refine
fmdl= ng_mk_cyl_models([0,2,.2],[0],[]);
3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1
fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]);
3D cylinder. Radius = 1. 6 circ elecs with elec refinement
fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]);
3D cylinder. Radius = 1. 7 rect elecs with no refinement
fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]);
2D cylinder. Radius = 1. 11 rect elecs with refinement
fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]);
2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement
fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]);
2D cylinder. elecs at 0, 90 and 120 degrees
fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]);
2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees
fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]);
3D cylinder. elecs at 0, 30, 60, 90 in planes
fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]);
3D cylinder. Various elecs at 0, 30, 60, 90 in planes
el_pos = [0,0.5;30,1;60,1.5;90,2.0];
el_sz = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
fmdl= ng_mk_cyl_models(3,el_pos,el_sz);
Simple 3D cylinder with a ball
extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'}
[fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra);
img= eidors_obj('image','ball'); img.fwd_model= fmdl;
img.elem_data(mat_idx{1}) = 1; img.elem_data(mat_idx{2}) = 2;
3D cylinder with 8 electrodes and cube
extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'}
[fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra);
3D cylinder with inner cylinder
extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'}
[fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra);
2D cylinder with 8 electrodes and hole
extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'}
fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra);
2D cylinder with 9 electrodes and inner cylinder
extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'}
fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra);
img= eidors_obj('image','ball'); img.fwd_model= fmdl;
ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2;
img.elem_data = 1 + 0.1*(ctr<0.2^2);0001 function [fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ... 0002 elec_shape, extra_ng_code); 0003 % NG_MAKE_CYL_MODELS: create cylindrical models using netgen 0004 %[fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ... 0005 % elec_shape, extra_ng_code); 0006 % INPUT: 0007 % cyl_shape = {height, [radius, [maxsz]]} 0008 % if height = 0 -> calculate a 2D shape 0009 % radius (OPT) -> (default = 1) 0010 % maxsz (OPT) -> max size of mesh elems (default = course mesh) 0011 % 0012 % ELECTRODE POSITIONS: 0013 % elec_pos = [n_elecs_per_plane,z_planes] 0014 % OR 0015 % elec_pos = [degrees,z] centres of each electrode (N_elecs x 2) 0016 % 0017 % ELECTRODE SHAPES:: 0018 % elec_shape = [width,height, maxsz] % Rectangular elecs 0019 % OR 0020 % elec_shape = [radius, 0, maxsz ] % Circular elecs 0021 % OR 0022 % elec_shape = [0, 0, maxsz ] % Point elecs 0023 % (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask) 0024 % 0025 % Specify either a common electrode shape or for each electrode 0026 % 0027 % EXTRA_NG_CODE 0028 % string of extra code to put into netgen geo file. Normally this 0029 % would be to insert extra materials into the space 0030 % 0031 % OUTPUT: 0032 % fmdl - fwd_model object 0033 % mat_idx - indices of materials (if extra_ng_code is used) 0034 % Note mat_idx does not work in 2D. Netgen does not provide it. 0035 % 0036 % 0037 % USAGE EXAMPLES: 0038 % Simple 3D cylinder. Radius = 1. No electrodes 0039 % fmdl= ng_mk_cyl_models(3,[0],[]); 0040 % Simple 2D cylinder. Radius = 2. Set minsize to refine 0041 % fmdl= ng_mk_cyl_models([0,2,.2],[0],[]); 0042 % 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1 0043 % fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]); 0044 % 3D cylinder. Radius = 1. 6 circ elecs with elec refinement 0045 % fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]); 0046 % 3D cylinder. Radius = 1. 7 rect elecs with no refinement 0047 % fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]); 0048 % 2D cylinder. Radius = 1. 11 rect elecs with refinement 0049 % fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]); 0050 % 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement 0051 % fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]); 0052 % 2D cylinder. elecs at 0, 90 and 120 degrees 0053 % fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]); 0054 % 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees 0055 % fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]); 0056 % 3D cylinder. elecs at 0, 30, 60, 90 in planes 0057 % fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]); 0058 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes 0059 % el_pos = [0,0.5;30,1;60,1.5;90,2.0]; 0060 % el_sz = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5]; 0061 % fmdl= ng_mk_cyl_models(3,el_pos,el_sz); 0062 % Simple 3D cylinder with a ball 0063 % extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'} 0064 % [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 0065 % img= eidors_obj('image','ball'); img.fwd_model= fmdl; 0066 % img.elem_data(mat_idx{1}) = 1; img.elem_data(mat_idx{2}) = 2; 0067 % 3D cylinder with 8 electrodes and cube 0068 % extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'} 0069 % [fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra); 0070 % 3D cylinder with inner cylinder 0071 % extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'} 0072 % [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 0073 % 2D cylinder with 8 electrodes and hole 0074 % extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'} 0075 % fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra); 0076 % 2D cylinder with 9 electrodes and inner cylinder 0077 % extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'} 0078 % fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra); 0079 % img= eidors_obj('image','ball'); img.fwd_model= fmdl; 0080 % ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2; 0081 % img.elem_data = 1 + 0.1*(ctr<0.2^2); 0082 0083 % (C) Andy Adler, 2009. (C) Alistair Boyle 2013. Licenced under GPL v2 or v3 0084 % $Id: ng_mk_cyl_models.m 6926 2024-05-31 15:34:13Z bgrychtol $ 0085 0086 if ischar(cyl_shape) && strcmp(cyl_shape,'UNIT_TEST'); do_unit_test; return; end 0087 0088 if nargin < 4; extra_ng_code = {'',''}; end 0089 copt.cache_obj = { cyl_shape, elec_pos, elec_shape, extra_ng_code}; 0090 copt.fstr = 'ng_mk_cyl_models'; 0091 copt.cache_on_ng_opt = true; 0092 args = {cyl_shape, elec_pos, elec_shape, extra_ng_code}; 0093 0094 fmdl = eidors_cache(@mk_cyl_model, args, copt); 0095 0096 mat_idx = fmdl.mat_idx; 0097 0098 function fmdl = mk_cyl_model( cyl_shape, elec_pos, elec_shape, extra_ng_code ); 0099 0100 fnstem = tempname; 0101 geofn= [fnstem,'.geo']; 0102 ptsfn= [fnstem,'.msz']; 0103 meshfn= [fnstem,'.vol']; 0104 0105 [tank_height, tank_radius, tank_maxh, is2D] = parse_shape(cyl_shape); 0106 [elecs, centres] = parse_elecs( elec_pos, elec_shape, ... 0107 tank_height, tank_radius, is2D ); 0108 0109 n_pts = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ... 0110 tank_maxh, elecs, extra_ng_code); 0111 if n_pts == 0 0112 call_netgen( geofn, meshfn); 0113 else 0114 call_netgen( geofn, meshfn, ptsfn); 0115 end 0116 0117 fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []); 0118 0119 % delete(geofn); delete(meshfn); delete(ptsfn); % remove temp files 0120 if is2D 0121 fmdl = mdl2d_from3d(fmdl); 0122 end 0123 0124 % convert CEM to PEM if so configured 0125 % TODO shunt model is unsupported 0126 if isfield(fmdl,'electrode'); 0127 fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes); 0128 end 0129 0130 % standard field order 0131 fmdl = eidors_obj('fwd_model',fmdl); 0132 0133 % for the newest netgen, we can't call msz file unless there are 0134 % actually points in it (ie. empty msz files break netgen) 0135 function n_pts_elecs = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ... 0136 tank_maxh, elecs, extra_ng_code); 0137 fid=fopen(geofn,'w'); 0138 write_header(fid,tank_height,tank_radius,tank_maxh,extra_ng_code); 0139 0140 n_elecs = length(elecs); 0141 % elecs(i).pos = [x,y,z] 0142 % elecs(i).shape = 'C' or 'R' 0143 % elecs(i).dims = [radius] or [width,height] 0144 % elecs(i).maxh = '-maxh=#' or ''; 0145 pts_elecs_idx = []; 0146 0147 for i=1:n_elecs 0148 name = sprintf('elec%04d',i); 0149 pos = elecs(i).pos; 0150 switch elecs(i).shape 0151 case 'C' 0152 write_circ_elec(fid,name, pos, pos, ... 0153 elecs(i).dims, tank_radius, elecs(i).maxh); 0154 case 'R' 0155 write_rect_elec(fid,name, pos, pos, ... 0156 elecs(i).dims, tank_radius, elecs(i).maxh); 0157 case 'P' 0158 % I had the good idea of trying to specify points for the point electrodes, 0159 % but netgen doesn't really listen to this. So instead, it only puts points 0160 % close to where you ask. Instead, specifc a rectangular elec where you want it. 0161 if 0 % OLD technique - keep in case we can figure out netgen better 0162 pts_elecs_idx = [ pts_elecs_idx, i]; 0163 continue; % DON'T print solid cyl 0164 else 0165 write_rect_elec(fid,name, pos, pos, ... 0166 elecs(i).dims, tank_radius, elecs(i).maxh); 0167 end 0168 0169 otherwise; error('huh? shouldnt get here'); 0170 end 0171 fprintf(fid,'solid cyl%04d = bigcyl and %s; \n',i,name); 0172 end 0173 0174 % SHOULD tank_maxh go here? - right now it seems to make the whole model refined 0175 fprintf(fid,'tlo bigcyl;\n'); 0176 for i=1:n_elecs 0177 if any(i == pts_elecs_idx); continue; end 0178 fprintf(fid,'tlo cyl%04d cyl -col=[1,0,0];\n ',i); 0179 end 0180 0181 for i=1:length(extra_ng_code)-1 0182 if ~isempty(extra_ng_code{i}) 0183 fprintf(fid,'tlo %s -col=[0,1,0];\n',extra_ng_code{i}); 0184 end 0185 end 0186 0187 fclose(fid); % geofn 0188 % From Documentation: Syntax is 0189 % np 0190 % x1 y1 z1 h1 0191 % x2 y2 z2 h2 0192 n_pts_elecs= length(pts_elecs_idx); 0193 fid=fopen(ptsfn,'w'); 0194 fprintf(fid,'%d\n',n_pts_elecs); 0195 for i = pts_elecs_idx; 0196 posxy = elecs(i).pos(1:2); 0197 fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) ); 0198 end 0199 fclose(fid); % ptsfn 0200 0201 function [tank_height, tank_radius, tank_maxh, is2D] = ... 0202 parse_shape(cyl_shape); 0203 tank_height = cyl_shape(1); 0204 tank_radius = 1; 0205 tank_maxh = 0; 0206 is2D = 0; 0207 0208 if length(cyl_shape)>1; 0209 tank_radius=cyl_shape(2); 0210 end 0211 if length(cyl_shape)>2; 0212 tank_maxh =cyl_shape(3); 0213 end 0214 if tank_height==0; 0215 is2D = 1; 0216 0217 %Need some width to let netgen work, but not too much so 0218 % that it meshes the entire region 0219 tank_height = tank_radius/5; % initial extimate 0220 if tank_maxh>0 0221 tank_height = min(tank_height,2*tank_maxh); 0222 end 0223 end 0224 0225 % ELECTRODE POSITIONS: 0226 % elec_pos = [n_elecs_per_plane,z_planes] 0227 % OR 0228 % elec_pos = [degrees,z] centres of each electrode (N_elecs x 2) 0229 % 0230 % ELECTRODE SHAPES:: 0231 % elec_shape = [width,height, {maxsz}] % Rectangular elecs 0232 % OR 0233 % elec_shape = [radius, {0, maxsz} ] % Circular elecs 0234 % maxsz (OPT) -> max size of mesh elems (default = courase mesh) 0235 % 0236 % OUTPUT: 0237 % elecs(i).pos = [x,y,z] 0238 % elecs(i).shape = 'C' or 'R' 0239 % elecs(i).dims = [radius] or [width,height] 0240 % elecs(i).maxh = '-maxh=#' or ''; 0241 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, hig, rad, is2D ); 0242 0243 n_elecs= size(elec_pos,1); 0244 0245 if n_elecs == 0 0246 elecs= struct([]); % empty 0247 centres= []; 0248 return; 0249 end 0250 0251 if is2D 0252 elec_pos(:,2) = hig/2; 0253 end 0254 0255 % It never makes sense to specify only one elec 0256 % So elec_pos means the number of electrodes in this case 0257 if size(elec_pos,1) == 1 0258 % Parse elec_pos = [n_elecs_per_plane,z_planes] 0259 n_elecs= elec_pos(1); % per plane 0260 th = linspace(0,2*pi, n_elecs+1)'; th(end)=[]; 0261 0262 on_elecs = ones(n_elecs, 1); 0263 el_th = []; 0264 el_z = []; 0265 for i=2:length(elec_pos) 0266 el_th = [el_th; th]; 0267 el_z = [el_z ; on_elecs*elec_pos(i)]; 0268 end 0269 else 0270 el_th = elec_pos(:,1)*2*pi/360; 0271 el_z = elec_pos(:,2); 0272 end 0273 0274 n_elecs= size(el_z,1); 0275 0276 if size(elec_shape,1) == 1 0277 elec_shape = ones(n_elecs,1) * elec_shape; 0278 end 0279 0280 for i= 1:n_elecs 0281 row = elec_shape(i,:); 0282 elecs(i) = elec_spec( row, is2D, hig, rad ); 0283 end 0284 0285 % %MC FIX 05.07.13 - COORDINATE SYSTEM - THETA = 0 AT 3PM AND THETA INCREASE ANTICLOCK 0286 % NO - we will use clockwise coordinate system 0287 centres = [rad*sin(el_th),rad*cos(el_th),el_z]; 0288 0289 for i= 1:n_elecs; elecs(i).pos = centres(i,:); end 0290 0291 if n_elecs == 0 0292 elecs= struct([]); % empty 0293 end 0294 0295 function elec = elec_spec( row, is2D, hig, rad ) 0296 if is2D 0297 if row(1) == 0; 0298 elec.shape = 'P'; 0299 % To create a PEM, we make a square and take the corner. This isn't perfect, since 0300 % the elec isn't quite where we asked for it, but that's as good is I can do. I tried 0301 % asking for two rectangles to touch, but that freaks netgen out. 0302 elec.dims = [rad/20, hig]; 0303 else 0304 elec.shape = 'R'; 0305 elec.dims = [row(1),hig]; 0306 end 0307 else 0308 if row(1) == 0 0309 elec.shape = 'P' 0310 elec.dims = [rad/20, hig/10]; 0311 elseif length(row)<2 || row(2) == 0 % Circular electrodes 0312 elec.shape = 'C'; 0313 elec.dims = row(1); 0314 elseif row(2)>0 % Rectangular electrodes 0315 elec.shape = 'R'; 0316 elec.dims = row(1:2); 0317 else 0318 error('negative electrode width'); 0319 end 0320 end 0321 0322 if length(row)>=3 && row(3) > 0 0323 elec.maxh = sprintf('-maxh=%f', row(3)); 0324 else 0325 elec.maxh = ''; 0326 end 0327 0328 function write_header(fid,tank_height,tank_radius,maxsz,extra); 0329 if maxsz==0; 0330 maxsz = ''; 0331 else 0332 maxsz = sprintf('-maxh=%f',maxsz); 0333 end 0334 0335 extra_ng= ''; 0336 for i=1:length(extra)-1 0337 if ~isempty( extra{i} ) 0338 extra_ng = sprintf(' %s and (not %s) ', ... 0339 extra_ng,extra{i}); 0340 end 0341 end 0342 0343 fprintf(fid,'#Automatically generated by ng_mk_cyl_models\n'); 0344 fprintf(fid,'algebraic3d\n'); 0345 fprintf(fid,'%s\n',extra{end}); % Define extra stuff here 0346 fprintf(fid,'solid cyl=cylinder (0,0,0;0,0,%6.2f;%6.2f); \n', ... 0347 tank_height, tank_radius); 0348 fprintf(fid,['solid bigcyl= plane(0,0,0;0,0,-1)\n' ... 0349 'and plane(0,0,%6.2f;0,0,1)\n' ... 0350 'and cyl %s %s;\n'],tank_height,extra_ng,maxsz); 0351 0352 0353 function write_rect_elec(fid,name,c, dirn,wh,d,maxh) 0354 % writes the specification for a netgen cuboid on fid, named name, centerd on c, 0355 % in the direction given by vector dirn, 0356 % hw = [height, width] and depth d 0357 % direction is in the xy plane 0358 w = wh(1); h= wh(2); 0359 dirn(3) = 0; dirn = dirn/norm(dirn); 0360 dirnp = [-dirn(2),dirn(1),0]; 0361 dirnp = dirnp/norm(dirnp); 0362 0363 bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2]; 0364 tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2]; 0365 fprintf(fid,'solid %s = ', name); 0366 fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ... 0367 bl(1),bl(2),bl(3)); 0368 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ... 0369 bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0); 0370 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ... 0371 bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0); 0372 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ... 0373 tr(1),tr(2),tr(3)); 0374 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ... 0375 tr(1),tr(2),tr(3),dirn(1),dirn(2),0); 0376 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f )%s;\n', ... 0377 tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh); 0378 0379 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh) 0380 % writes the specification for a netgen cylindrical rod on fid, 0381 % named name, centerd on c, 0382 % in the direction given by vector d, radius rd lenght ln 0383 % direction is in the xy plane 0384 % the direction vector 0385 dirn(3) = 0; dirn = dirn/norm(dirn); 0386 0387 % I would divide by 2 here (shorted tube in cyl), but ng doesn't like 0388 % That - it fails for 16 (but no 15 or 17) electrodes 0389 inpt = c - dirn.*(ln/1); 0390 outpt =c + dirn.*(ln/1); 0391 0392 fprintf(fid,'solid %s = ', name); 0393 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ... 0394 inpt(1),inpt(2),inpt(3),-dirn(1),-dirn(2),-dirn(3)); 0395 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ... 0396 outpt(1),outpt(2),outpt(3),dirn(1),dirn(2),dirn(3)); 0397 fprintf(fid,' cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) %s;\n', ... 0398 inpt(1),inpt(2),inpt(3),outpt(1),outpt(2),outpt(3), rd,maxh); 0399 0400 0401 function electrode = pem_from_cem(elecs, electrode, nodes) 0402 % elecs = electrode structure of model, from the parse_elecs function 0403 % electrode = the forward electrode model 0404 % nodes = the coordinates for the nodes 0405 % Can only have one node per electrode so we get a Point Electrode Model. 0406 % Choose the node with the greatest angle, so we atlest pick a consistent 0407 % side of the electrode: NetGen seems to give a random order to the nodes 0408 % in the electrode listing so we can't just pick the first one. 0409 % The nodes aside from those on the edges are not garanteed to be at any 0410 % particular location, so won't be consistent between meshes. 0411 % TODO should probably also adjust contact impedance too: its found later 0412 % by taking the average of the edges around the PEM's node, and those 0413 % will vary for each mesh -- should adjust so all electrodes get a 0414 % consistent effective impedance later. 0415 Ne = length(electrode); 0416 for i = 1:Ne 0417 if elecs(i).shape == 'P' 0418 % find the angles of the nodes for this electrode relative to (0,0) 0419 xy = nodes(electrode(i).nodes,:); 0420 ang = atan2(xy(:,2),xy(:,1)); 0421 % if the angles cover more than 180 degrees, must be an angle 0422 % roll-over from -pi to +pi, so take all the negative angles 0423 % and move them up 0424 if (max(ang) - min(ang)) > pi 0425 ang = ang + (ang <0)*2*pi; 0426 end 0427 % choose the counter-clockwise most node only 0428 % subtract the height so we get bottom left. This is OK, since we have rect elecs 0429 if size(xy,2) == 3 ; ang = ang - xy(:,3); end 0430 [jnk, ind] = max(ang); 0431 electrode(i).nodes = electrode(i).nodes(ind); 0432 end 0433 end 0434 0435 function do_unit_test 0436 for tn = 1:do_test_number(0) 0437 fprintf('>>> ng_mk_cyl_models: TEST #%d\n',tn); 0438 fmdl= do_test_number(tn); 0439 show_fem(fmdl); 0440 end 0441 0442 function fmdl= do_test_number(tn) 0443 eidors_msg('ng_mk_cyl_models: UNIT_TEST #%d',tn,1); 0444 switch tn 0445 case 1; 0446 % Simple 3D cylinder. Radius = 1. No electrodes 0447 fmdl= ng_mk_cyl_models(3,[0],[]); 0448 0449 case 2; 0450 % Simple 2D cylinder. Radius = 2. Set minsize to refine 0451 fmdl= ng_mk_cyl_models([0,2,.2],[0],[]); 0452 0453 case 3; 0454 % 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1 0455 fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]); 0456 0457 case 4; 0458 % 3D cylinder. Radius = 1. 6 circ elecs with elec refinement 0459 fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]); 0460 0461 case 5; 0462 % 3D cylinder. Radius = 1. 7 rect elecs with no refinement 0463 fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]); 0464 0465 case 6; 0466 % 2D cylinder. Radius = 1. 11 rect elecs with refinement 0467 fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]); 0468 0469 case 7; 0470 % 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement 0471 fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]); 0472 0473 case 8; 0474 % 2D cylinder. elecs at 0, 90 and 120 degrees 0475 fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]); 0476 0477 case 9; 0478 % 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees 0479 fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]); 0480 0481 case 10; 0482 % 3D cylinder. elecs at 0, 30, 60, 90 in planes 0483 fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]); 0484 0485 case 11; 0486 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes 0487 el_pos = [0,0.5;30,1;60,1.5;90,2.0]; 0488 el_sz = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5]; 0489 fmdl= ng_mk_cyl_models(3,el_pos,el_sz); 0490 0491 case 12; 0492 % Simple 3D cylinder with a ball 0493 extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'}; 0494 fmdl= ng_mk_cyl_models(3,[0],[],extra); 0495 img= mk_image(fmdl, 1); 0496 img.elem_data(fmdl.mat_idx{2}) = 2; 0497 0498 case 13; 0499 % 3D cylinder with 8 electrodes and cube 0500 extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'}; 0501 [fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra); 0502 0503 case 14; 0504 % 3D cylinder with inner cylinder 0505 extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'}; 0506 [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 0507 0508 case 15; 0509 % 2D cylinder with 8 electrodes and hole 0510 extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'}; 0511 fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra); 0512 0513 case 16; 0514 % 2D cylinder with 9 electrodes and inner cylinder 0515 extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'}; 0516 fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra); 0517 img= eidors_obj('image','ball'); img.fwd_model= fmdl; 0518 ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2; 0519 img.elem_data = 1 + 0.1*(ctr<0.2^2); 0520 0521 case 17; 0522 % Simple 3D cylinder with a ball 0523 extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'}; 0524 ng_write_opt('MSZBRICK',[-1,1,-1,1,2,3,0.1]); 0525 fmdl= ng_mk_cyl_models(3,[0],[],extra); 0526 delete('ng.opt'); 0527 img= mk_image(fmdl, 1); 0528 img.elem_data(fmdl.mat_idx{2}) = 2; 0529 0530 case 0; fmdl = 17; %%%% RETURN MAXIMUM 0531 otherwise; 0532 error('huh?') 0533 end 0534 0535 if exist('img'); fmdl = img; end