0001 function [fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0002 elec_shape, extra_ng_code);
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074 if ischar(ellip_shape) && strcmp(ellip_shape,'UNIT_TEST'), do_unit_test, return, end
0075 if nargin < 4; extra_ng_code = {'',''}; end
0076
0077 copt.cache_obj = { ellip_shape, elec_pos, elec_shape, extra_ng_code};
0078 copt.fstr = 'ng_mk_ellip_models';
0079 args = {ellip_shape, elec_pos, elec_shape, extra_ng_code};
0080
0081 fmdl = eidors_cache(@mk_ellip_model,args,copt);
0082
0083 mat_idx = fmdl.mat_idx;
0084
0085 function fmdl = mk_ellip_model( ellip_shape, elec_pos, elec_shape, extra_ng_code );
0086
0087 fnstem = tempname;
0088 geofn= [fnstem,'.geo'];
0089 ptsfn= [fnstem,'.msz'];
0090 meshfn= [fnstem,'.vol'];
0091
0092 [tank_height, tank_radius, tank_maxh, is2D] = parse_shape(ellip_shape);
0093 [elecs, centres] = parse_elecs( elec_pos, elec_shape, ...
0094 tank_height, tank_radius, is2D );
0095
0096 n_pts = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0097 tank_maxh, elecs, extra_ng_code);
0098 if n_pts == 0
0099 call_netgen( geofn, meshfn);
0100 else
0101 call_netgen( geofn, meshfn, ptsfn);
0102 end
0103
0104 fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0105
0106 delete(geofn); delete(meshfn); delete(ptsfn);
0107 if is2D
0108 fmdl = mdl2d_from3d(fmdl);
0109 end
0110
0111
0112
0113 if isfield(fmdl,'electrode');
0114 fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0115 end
0116
0117
0118 function n_pts_elecs = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0119 tank_maxh, elecs, extra_ng_code);
0120 fid=fopen(geofn,'w');
0121 write_header(fid,tank_height,tank_radius,tank_maxh,extra_ng_code);
0122
0123 n_elecs = length(elecs);
0124
0125
0126
0127
0128 pts_elecs_idx = [];
0129
0130 for i=1:n_elecs
0131 name = sprintf('elec%04d',i);
0132 pos = elecs(i).pos;
0133
0134 ab = tank_radius(1)/tank_radius(2);
0135 dirn= pos.*[inv(ab), ab, 0 ];
0136 switch elecs(i).shape
0137 case 'C'
0138 write_circ_elec(fid,name, pos, dirn, ...
0139 elecs(i).dims, tank_radius, elecs(i).maxh);
0140 case 'R'
0141 write_rect_elec(fid,name, pos, dirn, ...
0142 elecs(i).dims, tank_radius, elecs(i).maxh);
0143 case 'P'
0144 if 0
0145 pts_elecs_idx = [ pts_elecs_idx, i];
0146 continue;
0147 end
0148 write_rect_elec(fid,name, pos, dirn, ...
0149 elecs(i).dims, tank_radius, elecs(i).maxh);
0150
0151 otherwise; error('huh? shouldnt get here');
0152 end
0153 fprintf(fid,'solid cyl%04d = %s and not bigcyl; \n',i,name);
0154 end
0155
0156
0157 fprintf(fid,'tlo bigcyl;\n');
0158 for i=1:n_elecs
0159 if any(i == pts_elecs_idx); continue; end
0160 fprintf(fid,'tlo cyl%04d cyl -col=[1,0,0];\n ',i);
0161 end
0162
0163 for i=1:length(extra_ng_code)-1
0164 if ~isempty(extra_ng_code{i})
0165 fprintf(fid,'tlo %s -col=[0,1,0];\n',extra_ng_code{i});
0166 end
0167 end
0168
0169 fclose(fid);
0170
0171
0172
0173
0174 n_pts_elecs= length(pts_elecs_idx);
0175 fid=fopen(ptsfn,'w');
0176 fprintf(fid,'%d\n',n_pts_elecs);
0177 for i = pts_elecs_idx;
0178 posxy = elecs(i).pos(1:2);
0179 fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0180 end
0181 fclose(fid);
0182
0183 function [tank_height, tank_radius, tank_maxh, is2D] = ...
0184 parse_shape(cyl_shape);
0185 tank_height = cyl_shape(1);
0186 tank_radius = [1,1];
0187 tank_maxh = 0;
0188 is2D = 0;
0189 lcs = length(cyl_shape);
0190
0191 if lcs == 2
0192 tank_radius(1)=cyl_shape(2);
0193 elseif lcs >= 3
0194 tank_radius=cyl_shape(2:3);
0195 if diff(tank_radius) == 0;
0196 warning(['Using ng_mk_ellip_models to create cylinder. This may fail for '...
0197 'even electrode numbers. Recommend use ng_mk_cyl_models']);
0198 end
0199 end
0200 if length(cyl_shape)>=4;
0201 tank_maxh =cyl_shape(4);
0202 end
0203 if tank_height==0;
0204 is2D = 1;
0205
0206
0207
0208 tank_height = min(tank_radius)/5;
0209 if tank_maxh>0
0210 tank_height = min(tank_height,2*tank_maxh);
0211 end
0212 end
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, hig, rad, is2D );
0231
0232 if is2D
0233 elec_pos(:,2) = hig/2;
0234 end
0235
0236
0237
0238 if size(elec_pos,1) == 1
0239
0240 n_elecs= elec_pos(1);
0241 th = ellip_space_elecs( n_elecs, rad );
0242
0243 on_elecs = ones(n_elecs, 1);
0244 el_th = [];
0245 el_z = [];
0246 for i=2:length(elec_pos)
0247 el_th = [el_th; th];
0248 el_z = [el_z ; on_elecs*elec_pos(i)];
0249 end
0250 else
0251 el_th = elec_pos(:,1)*2*pi/360;
0252 el_z = elec_pos(:,2);
0253 end
0254
0255 n_elecs= size(el_z,1);
0256
0257 if size(elec_shape,1) == 1
0258 elec_shape = ones(n_elecs,1) * elec_shape;
0259 end
0260
0261 for i= 1:n_elecs
0262 row = elec_shape(i,:);
0263 elecs(i) = elec_spec( row, is2D, hig, rad );
0264 end
0265
0266
0267 centres = [rad(1)*sin(el_th),rad(2)*cos(el_th),el_z];
0268 for i= 1:n_elecs; elecs(i).pos = centres(i,:); end
0269
0270 if n_elecs == 0
0271 elecs= struct([]);
0272 end
0273
0274
0275 function th = ellip_space_elecs( n_elecs, rad )
0276
0277
0278
0279
0280 if n_elecs==0; th=[]; return; end
0281
0282 th = linspace(0,2*pi, 100*(n_elecs)); th(1)=[];
0283 len = cumsum( sqrt( rad(1)*cos(th).^2 + rad(2)*sin(th).^2 ) );
0284 len = len/max(len);
0285 xi = linspace(0,1,n_elecs+1); xi(1)= []; xi(end)=[];
0286 yi = interp1(len,th,xi);
0287
0288 th= [0;yi(:)];
0289 for exact = 0:3;
0290 eth = exact/2*pi;
0291 ff = abs(th-eth)<1e-3;
0292 th(ff) = eth;
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
0300
0301
0302 elec.dims = [min(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 = [min(rad)/20, hig/10];
0311 elseif length(row)<2 || row(2) == 0
0312 elec.shape = 'C';
0313 elec.dims = row(1);
0314 elseif row(2)>0
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
0329 function write_header(fid,tank_height,tank_radius,maxsz,extra);
0330 if maxsz==0;
0331 maxsz = '';
0332 else
0333 maxsz = sprintf('-maxh=%f',maxsz);
0334 end
0335
0336 extra_ng= '';
0337 for i=1:length(extra)-1
0338 if ~isempty( extra{i} )
0339 extra_ng = sprintf(' %s and (not %s) ', ...
0340 extra_ng,extra{i});
0341 end
0342 end
0343
0344 fprintf(fid,'#Automatically generated by ng_mk_ellip_models\n');
0345 fprintf(fid,'algebraic3d\n');
0346 fprintf(fid,'%s\n',extra{end});
0347 fprintf(fid,'solid cyl=ellipticcylinder (0,0,0;%.4f,0,0;0,%.2f8,0); \n', ...
0348 tank_radius);
0349 fprintf(fid,['solid bigcyl= plane(0,0,0;0,0,-1)\n' ...
0350 'and plane(0,0,%.4f;0,0,1)\n' ...
0351 'and cyl %s %s;\n'],tank_height,extra_ng,maxsz);
0352
0353
0354 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0355
0356
0357
0358
0359 d= min(d);
0360 w = wh(1); h= wh(2);
0361 dirn(3) = 0; dirn = dirn/norm(dirn);
0362 dirnp = [-dirn(2),dirn(1),0];
0363 dirnp = dirnp/norm(dirnp);
0364
0365 bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
0366 tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
0367 fprintf(fid,'solid %s = ', name);
0368 fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
0369 bl(1),bl(2),bl(3));
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),-dirn(1),-dirn(2),0);
0372 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0373 bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
0374 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
0375 tr(1),tr(2),tr(3));
0376 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0377 tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
0378 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f )%s;\n', ...
0379 tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh);
0380
0381 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0382
0383
0384
0385
0386
0387 dirn(3) = 0; dirn = dirn/norm(dirn);
0388
0389 ln = min(ln);
0390
0391
0392
0393 inpt = c - dirn.*(ln/6);
0394 outpt =c + dirn.*(ln/6);
0395
0396 fprintf(fid,'solid %s = ', name);
0397 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0398 inpt(1),inpt(2),inpt(3),-dirn(1),-dirn(2),-dirn(3));
0399 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0400 outpt(1),outpt(2),outpt(3),dirn(1),dirn(2),dirn(3));
0401 fprintf(fid,' cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) %s;\n', ...
0402 inpt(1),inpt(2),inpt(3),outpt(1),outpt(2),outpt(3), rd,maxh);
0403
0404
0405 function electrode = pem_from_cem(elecs, electrode, nodes)
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419 Ne = length(electrode);
0420 for i = 1:Ne
0421 if elecs(i).shape == 'P'
0422
0423 xy = nodes(electrode(i).nodes,:);
0424 ang = atan2(xy(:,2),xy(:,1));
0425
0426
0427
0428 if (max(ang) - min(ang)) > pi
0429 ang = ang + (ang <0)*2*pi;
0430 end
0431
0432 if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0433 [jnk, ind] = max(ang);
0434 electrode(i).nodes = electrode(i).nodes(ind);
0435 end
0436 end
0437
0438
0439 function do_unit_test
0440
0441 fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]); show_fem(fmdl);
0442
0443
0444 fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
0445
0446
0447 fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
0448
0449
0450 fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
0451
0452
0453 fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]);
0454 th = 45* [2*pi/360];
0455 fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
0456
0457
0458 fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
0459
0460
0461 el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0462 el_sz = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0463 fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
0464
0465
0466 extra={'ball','solid ball = sphere(0.5,0.5,1;0.4);'};
0467 [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra);
0468 img= mk_image(fmdl, 1);
0469 img.elem_data(mat_idx{2}) = 2; show_fem(img);
0470
0471 b1 = 'solid ball1= sphere(0.5, 0.5,1;0.2);';
0472 b2 = 'solid ball2= sphere(0.5,-0.5,1;0.2);';
0473 extra = {'ball1','ball2',[b1,b2]};
0474 [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra);
0475 img= mk_image(fmdl, 1);
0476 img.elem_data(mat_idx{2}) = 2;
0477 img.elem_data(mat_idx{3}) = 0.5;
0478 show_fem(img);
0479