Description of ng_mk_extruded

0001 function [fmdl,mat_idx] = ng_mk_extruded_model(shape, elec_pos, elec_shape, ...
0002     extra_ng_code)
0003 % NG_MAKE_EXTRUDED_MODEL: create extruded models using netgen
0004 % [fmdl,mat_idx] = ng_mk_extruded_models(trunk_shape, elec_pos, ...
0005 %                 elec_shape, extra_ng_code);
0006 % INPUT:
0007 % trunk_shape = { height,[x,y],curve_type,maxsz}
0008 %   height      -> if height = 0 calculate a 2D model
0009 %   [x,y]       -> N-by-2 CLOCKWISE list of points defining the 2D shape
0010 %                  NOTE: Use a cell array to specify additional curves for
0011 %                  internal objects
0012 %   curve_type  -> 1 - interpret as vertices (default)
0013 %                  2 - interpret as splines with de Boor points at even
0014 %                  indices (legacy)
0015 %                  3 - interpolate points (piecewise polynomial
0016 %                  interpolation). Syntax [3, N] also specifies the number
0017 %                  of samples to create.
0018 %                  4 - interpolate points with Fourier descriptor. Syntax
0019 %                  [4, N] also specifies the number of samples to create.
0020 %                  NOTE: If additional curves are specified, curve_type can
0021 %                  also be a cell array. Otherwise, curve_type defaults to
0022 %                  1 for internal shapes.
0023 %   maxsz       -> max size of mesh elems (default = course mesh). If there
0024 %                  multiple curves, maxh can be an array specifying the max
0025 %                  size for each internal object.
0026 %
0027 % ELECTRODE POSITIONS:
0028 %  elec_pos = [n_elecs_per_plane,spacing,z_planes] where spacing is either
0029 %             0 for even spacing w.r.t angular positions (0,15,30... deg)
0030 %             or
0031 %             1 for equal distances between electrodes
0032 %             Any fractional part (e.g. 0.15) is interpreted as a starting
0033 %             position -- fraction of 2*pi for values spacing < 1 and
0034 %             fraction of total perimeter for spacing > 1.
0035 %     OR
0036 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0037 %
0038 % ELECTRODE SHAPES::
0039 %  elec_shape = [width,height, [maxsz, pem, discr]]  % Rectangular elecs
0040 %     OR
0041 %  elec_shape = [radius, [0, maxsz, pem, discr] ]  % Circular elecs
0042 %     radius      -> specify -1 for point electrodes
0043 %     maxsz (OPT) -> max size of mesh elems (default = course mesh),
0044 %                    ignored if <= 0
0045 %     pem  (OPT)  -> 1: Point Electrode Model, 0: Complete Electrode Model (default)
0046 %     discr (OPT) -> discretize electrode normals (0 = ignore = default)
0047 %                    Adjusting this value helps Netgen problems with
0048 %                    electrodes facing each other.
0049 %
0050 % Specify either a common electrode shape or for each electrode
0051 %
0052 % CITATION_REQUEST:
0053 % AUTHOR: B Grychtol et al.
0054 % TITLE: Impact of model shape mismatch on reconstruction quality in
0055 % Electrical Impedance Tomography
0056 % JOURNAL: IEEE Trans Med Imag
0057 % YEAR: 2012
0058 % VOL: 31
0059 % NUM: 9
0060 % DOI: 10.1109/TMI.2012.2200904
0061 % PDF: http://www.sce.carleton.ca/faculty/adler/publications/2012/
0062 %      grychtol-2012-model-shape-EIT.pdf
0063 
0064 % (C) Bartlomiej Grychtol, 2010. (C) Alistair Boyle, 2013. Licenced under GPL v2 or v3
0065 % $Id: ng_mk_extruded_model.m 7012 2024-11-26 15:52:57Z bgrychtol $
0066 
0067 % TODO: Implement control segments in the bit that writes the file.
0068 
0069 if ischar(shape) && strcmp(shape,'UNIT_TEST'); do_unit_test; return; end
0070 
0071 citeme(mfilename);
0072 
0073 if nargin < 4; extra_ng_code = {'',''}; end
0074 copt.cache_obj = { shape, elec_pos, elec_shape, extra_ng_code};
0075 copt.fstr = 'ng_mk_extruded_model';
0076 args = { shape, elec_pos, elec_shape, extra_ng_code};
0077 copt.cache_on_ng_opt = true;
0078 
0079 fmdl = eidors_cache(@mk_extruded_model, args, copt);
0080 
0081 mat_idx = fmdl.mat_idx;
0082 
0083 copt.args = {args, 'ng.opt'}; % algo cache on ng.opt
0084 
0085 
0086 function fmdl = mk_extruded_model(shape, elec_pos, elec_shape, ...
0087     extra_ng_code)
0088 
0089 fnstem = tempname;
0090 geofn= [fnstem,'.geo'];
0091 meshfn= [fnstem,'.vol'];
0092 
0093 [tank_height, tank_shape, tank_maxh, is2D] = parse_shape(shape);
0094 [elecs, centres] = parse_elecs(elec_pos, elec_shape, tank_shape, tank_height, is2D );
0095 write_geo_file(geofn, tank_height, tank_shape, ...
0096                tank_maxh, elecs, extra_ng_code);
0097            
0098 if 0% DEBUG SHAPE
0099    plot(tank_shape.vertices(:,1),tank_shape.vertices(:,2));
0100    hold on
0101    plot(centres(:,1),centres(:,2),'sk')
0102    for i = 1:size(elecs,2)
0103        dirn = elecs(i).normal;
0104        quiver(centres(i,1),centres(i,2),dirn(1),dirn(2),'k');
0105    end
0106    hold off
0107    axis equal
0108 end
0109            
0110 call_netgen( geofn, meshfn );
0111 
0112 fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', [],0.01,...
0113     @ng_remove_electrodes);
0114 
0115 delete(geofn); delete(meshfn); % remove temp files
0116 
0117 if is2D
0118     fmdl = mdl2d_from3d(fmdl);
0119 end
0120 
0121 
0122 
0123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0124 % TANK SHAPE (struct):
0125 %         vertices: [Nx2]
0126 %             size: 0.5* length of the diagonal of the containing rectangle
0127 %     edge_normals: [Nx2]
0128 %       vertex_dir: [Nx2] direction of vertex movement when scaling
0129 %         centroid: [x y]
0130 %   vertices_polar: [Nx2] Phi, r
0131 %           convex: [N] boolean array indicating external angle >= 180 deg
0132 %       curve_type: One of three values
0133 %                   1 - Normal, each point is a vertex
0134 %                   2 - Spline, all even points are de Boor points
0135 %                   3 - Same as 1 but will be converted to smooth
0136 %
0137 function [tank_height, tank_shape, tank_maxh, is2D] = parse_shape(shape)
0138     % parses the shape input
0139 
0140     %defaults
0141     is2D = false;
0142     tank_maxh = 0;
0143     tank_shape = [];
0144     tank_shape.curve_type = 1;
0145     curve_info = [];
0146 
0147     if iscell(shape) && length(shape)>2
0148         tank_height = shape{1};
0149         if ~iscell(shape{2})
0150             points = shape{2};
0151         else
0152             c = shape{2};
0153             points = c{1};
0154             if numel(shape{2}) > 1
0155                 tank_shape.additional_shapes = c(2:end);
0156             end
0157         end
0158         
0159         if ~iscell(shape{3})
0160             tank_shape.curve_type = shape{3};
0161             if iscell(tank_shape.curve_type)
0162                 tank_shape.curve_type = tank_shape.curve_type{1};
0163             end
0164         else
0165             c = shape{3};
0166             tank_shape.curve_type = c{1};
0167             if numel(shape{3}) > 1
0168                 tank_shape.additional_curve_type = c(2:end);
0169             end
0170         end
0171         
0172         if max(size(tank_shape.curve_type)) > 1
0173             curve_info = tank_shape.curve_type;
0174             tank_shape.curve_type = curve_info(1);
0175         end
0176 %         if length(shape) > 2
0177 %             tank_height = shape{1};
0178 %         end
0179         if length(shape) > 3
0180             maxh = shape{4};
0181             tank_maxh = maxh(1);
0182             if isfield(tank_shape,'additional_shapes')
0183                 N_curves = 1 + numel(tank_shape.additional_shapes);            
0184                 if numel(maxh) == 1
0185                     tank_maxh(2:N_curves) = maxh(1);
0186                 elseif numel(maxh) == N_curves;
0187                     tank_maxh = maxh;
0188                 else
0189                     error('length of maxh must either be 1 or equal to the total number of curves');
0190                 end
0191             end
0192         else % maxh not specified
0193             if isfield(tank_shape,'additional_shapes')
0194                 N_curves = 1 + numel(tank_shape.additional_shapes);
0195                 tank_maxh(2:N_curves) = tank_maxh;
0196             end
0197         end
0198     else
0199         points = shape;
0200     end
0201     
0202     spln_sgmnts = zeros(size(points)); %default
0203     if tank_shape.curve_type == 2
0204         [points, spln_sgmnts] = remove_linear_control_points(points);
0205     end
0206     
0207     if ~isempty(curve_info)
0208         N = curve_info(2);
0209     else
0210         N = 50;
0211     end
0212     points = interpolate(points,N, tank_shape.curve_type);
0213     spln_sgmnts = zeros(size(points)); % WHY??
0214     
0215     if isfield(tank_shape, 'additional_curve_type')
0216         for i = 1:numel(tank_shape.additional_curve_type)
0217             if numel(tank_shape.additional_curve_type{i}) == 1
0218                 N = 50;
0219             else
0220                 N = tank_shape.additional_curve_type{i}(2);
0221             end
0222             tank_shape.additional_shapes{i} = interpolate(...
0223                 tank_shape.additional_shapes{i},N, tank_shape.additional_curve_type{i}(1));
0224         end
0225     end
0226     
0227     % piecewise polynomial interpolation
0228     if tank_shape.curve_type == 3 
0229         if ~isempty(curve_info)
0230             n_samples = curve_info(2);
0231         else
0232             n_samples = 50;
0233         end
0234         points = interpolate_shape(points, n_samples);
0235         spln_sgmnts = zeros(size(points)); % now needs to be bigger
0236     end
0237 
0238     % Fourier descriptor interpolation
0239     if tank_shape.curve_type == 4
0240         if ~isempty(curve_info)
0241             n_samples = curve_info(2);
0242         else
0243             n_samples = 50;
0244         end
0245         points = fourier_interpolate_shape(points, n_samples);
0246         spln_sgmnts = zeros(size(points)); % now needs to be bigger
0247     end
0248     
0249     tank_shape.centroid = calc_centroid(points);
0250     tank_shape.spln_sgmnts = spln_sgmnts;
0251 
0252     tank_shape.vertices = points;
0253     % diagonal of the containing rectangle:
0254     range_points = max(points) - min(points);
0255     tank_shape.size = 0.5 * sqrt(sum(range_points.^2));
0256     
0257     if tank_height==0
0258         is2D = 1;
0259         % Need some width to let netgen work, but not too much so
0260         % that it meshes the entire region
0261         tank_height = tank_shape.size/5; % initial estimate
0262         if tank_maxh(1)>0
0263             tank_height = min(tank_height,2*tank_maxh(1));
0264         end
0265     end
0266 
0267 
0268     tank_shape.edge_normals = [];
0269     tank_shape.vertex_dir = [];
0270 
0271     tmp = points;
0272     tmp(end+1,:) = tmp(1,:); %duplicate first vertex at the end;
0273 
0274     edges = diff(tmp,1);
0275     tmp = [];
0276     % Normal to vector (x y) is (-y x).
0277     % It points outward for clockwise definition
0278     tmp = circshift(edges, [0 1]); %swap coords
0279     %normalize
0280     lngth = sqrt(sum(tmp.^2, 2));
0281     tmp(:,1) = -tmp(:,1) ./ lngth;
0282     tmp(:,2) = tmp(:,2)  ./ lngth;
0283     tank_shape.edge_normals = tmp;
0284 
0285     tank_shape.vertex_dir = calc_vertex_dir(points, edges, ...
0286         tank_shape.edge_normals);
0287     
0288     
0289     tmp = [];
0290     polar = zeros(size(points));
0291     for i = 1:length(points)
0292         tmp = points(i,:) - tank_shape.centroid;
0293         [polar(i,1) polar(i,2)]  = cart2pol(tmp(1),tmp(2));
0294     end
0295     tank_shape.vertices_polar = polar;
0296     
0297     tank_shape.convex = calc_convex(tank_shape.vertices);
0298     
0299     % debug plot
0300 if 0
0301     pts = edges./2 + points;
0302     plot(tank_shape.vertices(:,1),tank_shape.vertices(:,2),'-o'); hold on;
0303     plot(tank_shape.centroid(:,1),tank_shape.centroid(:,2),'+');
0304     plot(tank_shape.vertices(:,1)+0.05*tank_shape.vertex_dir(:,1),...
0305         tank_shape.vertices(:,2)+0.05*tank_shape.vertex_dir(:,2),'ro-')
0306     quiver(pts(:,1),pts(:,2),tank_shape.edge_normals(:,1),tank_shape.edge_normals(:,2));
0307     hold off
0308 end
0309     
0310     
0311 function new_points = interpolate(points, N, curve_type)
0312 switch curve_type
0313     case 3 
0314         % piecewise polynomial interpolation
0315         new_points = interpolate_shape(points, N);
0316     case 4
0317         % Fourier descriptor interpolation
0318         new_points = fourier_interpolate_shape(points, N);
0319     otherwise 
0320         % do nothing
0321         new_points = points;
0322 end  
0323     
0324     
0325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0326 % INPUT:
0327 % points - [2N x 2] defined vertices (odd) and control points (even)
0328 % OUTPUT:
0329 % points   - same as points but with linear control points removed
0330 % spln_sgmnts - boolean array indicating which segments are splines
0331 function [points, spln_sgmnts] = remove_linear_control_points(points)
0332 n_points = length(points);
0333 points(end+1,:) = points(1,:);
0334 spln_sgmnts(1:(n_points/2)) = 1;
0335 for i = 1:2:n_points
0336     a = (points(i+1,:) - points(i,:));
0337     a = a/norm(a);
0338     b = (points(i+2,:) - points(i,:));
0339     b = b/norm(b); 
0340     if a(1) == b(1) && a(2) == b(2)
0341         spln_sgmnts(i/2 + 0.5) = 0;
0342     end    
0343 end
0344 idx = find(spln_sgmnts==0) * 2;
0345 points(idx,:) = [];
0346 points(end,:) = [];
0347 
0348     
0349 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0350 % INPUT:
0351 % points - [N x 2] defined vertices
0352 % OUTPUT:
0353 % out    - interpolated vertices
0354 function out = interpolate_shape(points, n_points)
0355 % Quadratic spline interpolation of the points provided.
0356 
0357 
0358 [pp m] = piece_poly_fit(points);
0359 p = linspace(0,1,n_points+1)'; p(end) = [];
0360 [th xy] = piece_poly_fit(pp,0,p);
0361 tmp = [th xy];
0362 tmp = sortrows(tmp,-1);% ensure clockwise direction
0363 xy = tmp(:,2:3);
0364 
0365 out = xy + repmat(m, [n_points,1]);
0366 
0367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0368 % INPUT:
0369 % points - [N x 2] defined vertices
0370 % OUTPUT:
0371 % out    - interpolated vertices
0372 function out = fourier_interpolate_shape(points, n_points)
0373 % Quadratic spline interpolation of the points provided.
0374 
0375 
0376 pp = fourier_fit(points, size(points,1)-1); % don't want to overfit
0377 p = linspace(0,1,n_points+1)'; p(end) = [];
0378 xy = fourier_fit(pp,p);
0379 % [th r] = cart2pol(xy);
0380 % tmp = [th xy];
0381 % tmp = sortrows(tmp,-1);% ensure clockwise direction
0382 % xy = tmp(:,2:3);
0383 
0384 out = xy;% + repmat(m, [n_points,1]);
0385 
0386 
0387 function out = calc_vertex_dir(points, edges, edgnrm)
0388 %     calculate the direction of vertex movement if all edges are shifted
0389 %     outwards by 1 unit along their normals:
0390 
0391 %     duplicate last edge at the beginning
0392     edg = [edges(end,:) ; edges];
0393     edgnrm = [edgnrm(end,:) ; edgnrm];
0394 
0395     out = zeros(size(points));
0396     for i = 1:length(points)
0397         p1 = points(i,:) + edgnrm(i,:);
0398         p2 = points(i,:) + edgnrm(i+1,:);
0399 
0400         dir1(1) = edgnrm(i,2); dir1(2) = -edgnrm(i,1);
0401         dir2(1) = edgnrm(i+1,2); dir2(2) = -edgnrm(i+1,1);
0402         % if the edge directions are the same (accounting for round-off
0403         % error), return the edge normal.
0404         if isempty(find(abs(dir1 - dir2) > 1e-14))
0405             out(i,:) = edgnrm(i,:);
0406         else
0407             A = [dir1' , -dir2'];
0408             u = (p2 - p1)';
0409             x = A\u;
0410             out(i,:) = x(1) * dir1 + p1 - points(i,:);
0411         end
0412     end
0413 
0414 function out = calc_centroid(points)
0415 % Calculates the centroid of the shape
0416 % The algorithm identifies a middle point M within the shape and then uses it
0417 % to divide the shape into N triangles (N=number of vertices), calculates
0418 % the area and centroid of each traingle, and finally computes the centroid
0419 % of the shape as a mean of the centroids of the individual traingles
0420 % weighted by their area.
0421 
0422     % it never makes sense to have less than 3 points
0423     n_points = size(points,1);
0424     if  n_points == 3
0425         out = mean(points); % centroid of a triangle
0426         return
0427     end
0428 
0429     out = 0;
0430     pts = [points ; points(1,:)];
0431 
0432     % guess a point in the middle
0433     m = mean(points);
0434 
0435     if ~inpolygon(m(1),m(2),points(:,1),points(:,2))
0436         f1 = figure;
0437         set(f1,'Name', 'Select a point within the shape');
0438         plot(pts(:,1),pts(:,2));
0439         m = ginput(1);
0440         close(f1)
0441     end
0442 
0443     tmp = 0;
0444     tot_area = 0;
0445     for i = 1:n_points
0446         a = pts(i,:);
0447         b = pts(i+1,:);
0448         cntrd = (m + a + b)/3;
0449         area = 0.5 * abs(det([m 1; a 1; b 1]));
0450         tmp = tmp + cntrd*area;
0451         tot_area = tot_area + area;
0452     end
0453 
0454     out = tmp./tot_area;
0455 
0456 function out = calc_convex(verts)
0457 % Returns an array of boolean values for every vertex, true if the external
0458 % angle at this vertex is greater or equal to 180 degrees, false otherwise.
0459 % This marks the vertices which upset the convexity of the polygon and
0460 % require special treatment.
0461 
0462 n_verts = size(verts,1);
0463 tmp = [verts(end,:); verts; verts(1,:)];
0464 verts = tmp;
0465 
0466 for i = 2:n_verts+1
0467     v1 = [verts(i-1,:) - verts(i,:), 0];
0468     v2 = [verts(i+1,:) - verts(i,:), 0];
0469     cp = cross(v1,v2);
0470     out(i-1) = cp(3) >= 0;
0471 end
0472 
0473 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0474 % ELECTRODE POSITIONS:
0475 %  elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z_planes]
0476 %     OR
0477 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0478 %
0479 % ELECTRODE SHAPES::
0480 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0481 %     OR
0482 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0483 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0484 %
0485 % OUTPUT:
0486 %  elecs(i).pos   = [x,y,z]
0487 %  elecs(i).shape = 'C' or 'R'
0488 %  elecs(i).dims  = [radius] or [width,height]
0489 %  elecs(i).maxh  = '-maxh=#' or '';
0490 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, tank_shape, hig, is2D )
0491 
0492     if isempty(elec_pos)
0493         elecs = [];
0494         centres = [];
0495         return; % no electrodes, nothing to do
0496     end
0497     
0498    % temp fix
0499    rad = tank_shape.size;
0500 
0501    % It never makes sense to specify only one elec
0502    % So elec_pos means the number of electrodes in this case
0503    if size(elec_pos,1) == 1
0504        % Parse elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z_planes]
0505        if is2D
0506           % create an electrode at half height
0507           elec_pos(:,3) = hig/2;
0508        end
0509    
0510       n_elecs= elec_pos(1); % per plane
0511       offset = elec_pos(2) - floor(elec_pos(2));
0512       switch floor(elec_pos(2))
0513           case 0
0514               th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0515               th = th + offset*2*pi;
0516               ind = th >= 2*pi;
0517               th(ind) = th(ind) - 2*pi;
0518           case 1
0519               % piece_poly_fit doesn't seem to work very well
0520               if 1%tank_shape.curve_type == 4
0521                   pp = fourier_fit(tank_shape.vertices,...
0522                       size(tank_shape.vertices,1) - 1,tank_shape.vertices(1,:));
0523                   p = linspace(0,1,n_elecs+1)'; p(end) = [];
0524                   p = p + offset;
0525                   xy = fourier_fit(pp,p);
0526                  % NOTE, THIS IS A HACK. Some complicated shapes can't be
0527                  % described by angle alone
0528                   th = atan2(xy(:,2) - tank_shape.centroid(2), ...
0529                              xy(:,1) - tank_shape.centroid(1));
0530 
0531               elseif any( tank_shape.curve_type == [1,2,3] )
0532                   % I can't seem able to get the first electrode exactly on
0533                   % the first vertex
0534                   pp= piece_poly_fit(tank_shape.vertices);
0535                   p = linspace(1,0,n_elecs+1)'; p(end) = [];
0536                   off = offset*2*pi + tank_shape.vertices_polar(1,1);
0537                   th = piece_poly_fit(pp,off,p);
0538               else
0539                   error('curve_type unrecognized');
0540               end
0541           otherwise;
0542              error('Unrecognized value of curve_type specified in floor(elec_pos(2))')
0543       end
0544 
0545       on_elecs = ones(n_elecs, 1);
0546       el_th = []; 
0547       el_z  = []; 
0548       for i=3:length(elec_pos)
0549         el_th = [el_th; th];
0550         el_z  = [el_z ; on_elecs*elec_pos(i)];
0551       end
0552    else
0553       th = elec_pos(:,1)*2*pi/360;
0554       el_th = []; 
0555       el_z  = []; 
0556       for i=2:size(elec_pos,2)
0557         el_th = [el_th; th];
0558         el_z  = [el_z ; elec_pos(:,i)];
0559       end
0560    end
0561       
0562    n_elecs= size(el_z,1); 
0563 
0564    if size(elec_shape,1) == 1
0565       elec_shape = ones(n_elecs,1) * elec_shape;
0566    end
0567 
0568    for i= 1:n_elecs
0569      row = elec_shape(i,:); 
0570      elecs(i) = elec_spec( row, is2D, hig, rad );
0571    end
0572    
0573    
0574    %centres = [rad*sin(el_th),rad*cos(el_th),el_z];
0575    for i= 1:n_elecs; 
0576 %        switch tank_shape.curve_type
0577 %            case 1
0578                [centres(i,1:2), normal] = calc_elec_centre(tank_shape, el_th(i));
0579 %            case{2, 3}
0580 %                [centres(i,1:2), normal] = calc_elec_centre_spline(tank_shape, el_th(i));
0581 %            otherwise
0582 %                error('Unknown curve type');
0583 %        end
0584        centres(i,3) = el_z(i);
0585        elecs(i).pos  = centres(i,:);
0586        if elecs(i).discretize > 0
0587         % this bit is to prevent netgen choking on slightly misalligned electrodes
0588         [th,~] = cart2pol(normal(1),normal(2)); % ~ needed for octave
0589         frac = 2*pi /elecs(i).discretize ;
0590         th = frac * round( th / frac);
0591         [normal(1) normal(2)] = pol2cart(th,1);
0592        end
0593        elecs(i).normal = normal;
0594        
0595    end
0596 
0597    if n_elecs == 0
0598       elecs= struct([]); % empty
0599       centres= []; 
0600    end
0601 
0602    
0603    
0604     function [pos, normal] = calc_elec_centre(tank_shape, th)
0605         % The calculation relies on the theorem that if point D lies on a
0606         % line between B and C, but point A is not on that line, then:
0607         %   |BD|    |AB| sin(<DAB)
0608         %   ---- = ---------------
0609         %   |DC|    |AC| sin(<DAC)
0610         % Thus, B and C are vertices of our shape, A is its centroid and D
0611         % is the sought center of the electrode. All quantities on RHS are
0612         % known.
0613         
0614         % make sure th is between -pi and pi
0615         if th > pi; th = th - 2*pi; end
0616         
0617         
0618         vert_pol = tank_shape.vertices_polar; %[th, r]
0619         
0620      
0621         n_vert = size(vert_pol,1);
0622         vert_pol = [vert_pol , (1:n_vert)'];
0623         % Re-order the vertices -pi to pi. Now counter-clockwise
0624         vert_pol = sortrows(vert_pol,1); 
0625         % find the edge on which the elctrode lies. (Edge 1 is between
0626         % verticies 1 and 2)
0627         idx = find(vert_pol(:,1) > th, 1, 'first');
0628         if isempty(idx); idx = 1; end
0629         edg_no = vert_pol(idx,3);
0630         
0631         
0632         normal = tank_shape.edge_normals(edg_no,:);
0633               
0634         v1 = edg_no;
0635         if edg_no == n_vert % between the last and first vertex
0636             v2 = 1;
0637         else
0638             v2 = v1+1;
0639         end
0640         vert_pol = [];
0641         
0642         
0643         vert_pol = tank_shape.vertices_polar;
0644         vert = tank_shape.vertices;
0645         cntr = tank_shape.centroid;
0646         % position between vertices - see first comment
0647         AB = sqrt(sum( (vert(v1,:) - cntr).^2 ));
0648         AC = sqrt(sum( (vert(v2,:) - cntr).^2 ));
0649         DAB = abs(vert_pol(v1,1)-th); 
0650         if DAB > pi, DAB = abs( DAB - 2*pi); end; 
0651         DAC  = abs(vert_pol(v2,1)-th);
0652         if DAC > pi, DAC = abs( DAC - 2*pi); end;
0653         if DAC ~= 0
0654             ratio = AB * sin(DAB) / (AC * sin(DAC));
0655             pos = vert(v1,:) + ( ratio / (1 + ratio) ) * (vert(v2,:) - vert(v1,:));
0656         else
0657             pos = vert(v2,:);
0658         end
0659 
0660         
0661         
0662    function [pos, normal] = calc_elec_centre_spline(tank_shape, th)
0663         % The calculation proceeds by finding a common point between a line
0664         % from the centroid outwards and the equation of the relevant
0665         % quadratic spline segment defined using 3 control points
0666         
0667         % make sure th is between -pi and pi
0668         if th > pi; th = th - 2*pi; end 
0669         
0670         vert_pol = tank_shape.vertices_polar; %[th, r]
0671         
0672         % The number of vertices must be even, but just in case...
0673         if mod(size(vert_pol,1),2)
0674             error(['The number of points must be even. '...
0675                 'One de Boor control point for every vertex']);
0676         end
0677         
0678         % if the curve is defined as splines, every second point is not
0679         % actually a vertex. We remove them.
0680         if tank_shape.curve_type == 2 || tank_shape.curve_type == 3
0681             vert_pol(2:2:end,:) = [];
0682         end
0683       
0684         n_vert = size(vert_pol,1);
0685    
0686         vert_pol = [vert_pol , (1:n_vert)']; %excludes control points
0687         % Re-order the vertices -pi to pi. Now counter-clockwise
0688         vert_pol = sortrows(vert_pol,1); 
0689         % find the edge on which the electrode lies. Edge 1 is between
0690         % vertices 1 and 2.
0691         idx = find(vert_pol(:,1) > th, 1, 'first');
0692         if isempty(idx); idx = 1; end
0693         edg_no = vert_pol(idx,3);
0694         
0695         v1 = edg_no;
0696         if edg_no == n_vert % between the last and first vertex
0697             v2 = 1;
0698         else
0699             v2 = v1+1;
0700         end
0701         vert_pol = [];
0702         
0703         % correcting for the control points
0704         v1 = 2 * v1 - 1;
0705         v2 = 2 * v2 - 1;
0706         
0707         % the spline goes from point P0 to P2 such that P1-P0 is a tangent
0708         % at P0 and P2-P1 is a tangent at P2
0709         C = tank_shape.centroid;
0710         P0 = tank_shape.vertices(v1,:) - C;
0711         P1 = tank_shape.vertices(v1+1,:) - C; % control point
0712         P2 = tank_shape.vertices(v2,:) - C;
0713         
0714         
0715         % find the gradient of the line from centroid to electrode center:
0716         [x, y] = pol2cart(th, 1);
0717         % FIXME: This doesn't crash only because of round-off errors.
0718         g = y/x;
0719         % (because we subtracted the centroid from the vertices, the line
0720         % passes through the origin now)
0721         
0722         % the spline is f(t) = (1-t)^2 * P0 + 2t(1-t)P1 + t^2 * P2
0723         % which can also be expressed as
0724         f = @(t) (P2 - 2*P1 + P0)*t^2 + 2*(P1 - P0)*t + P0;
0725         % and it's derivative:
0726         df = @(t) 2*(P2 - 2*P1 + P0)*t + 2*(P1 - P0);
0727         % to find the value of t for which the line cross, we substitute
0728         % p0 = y0-ax0 for P0 and so on.
0729         p0 = P0(2) - g * P0(1);
0730         p1 = P1(2) - g * P1(1);
0731         p2 = P2(2) - g * P2(1);
0732         
0733         % thus we have a quadratic equation a*t^2 + b*t + c = 0 where
0734         a = (p2 - 2*p1 + p0);
0735         b = 2* (p1 - p0);
0736         c = p0;
0737         
0738         if abs(a) < 1e-10
0739             t = -c/b;
0740             pos = f(t) + C;
0741             tmp = df(t);
0742             normal = [-tmp(2), tmp(1)] / sqrt(sum(tmp.^2));
0743             return;
0744         end
0745         
0746         % the determinant is
0747         D = b^2 - 4*a*c;
0748         
0749         % find the roots
0750         if D == 0
0751             t = -b / (2 * a);
0752 
0753         elseif D > 0
0754             t1 = (-b - sqrt(D) ) / (2 * a);
0755             t2 = (-b + sqrt(D) ) / (2 * a);
0756             if t1 >= 0 && t1 <= 1
0757                 t = t1;
0758             else
0759                 t = t2;
0760             end
0761         else
0762             error('Something went wrong, cannot place electrode on spline');
0763         end
0764         
0765         pos = f(t) + C;
0766         tmp = df(t);
0767         normal = [-tmp(2), tmp(1)]/ sqrt(sum(tmp.^2));
0768 
0769    
0770    
0771 
0772 function elec = elec_spec( row, is2D, hig, rad )
0773   if     is2D
0774      if length(row)>=2 && row(2) == -1 % Point electrodes
0775         % Create rectangular electrodes with bottom, cw point where we want
0776         elec.shape = 'P' ;
0777         if length(row)>=3 && row(3) > 0
0778            elec.dims  =  row(3);
0779         else
0780            elec.dims  =  rad; % Make big if unspecified
0781         end
0782      else
0783         % create circular electrodes for now, rectangular not yet supported
0784 %         elec.shape = 'C';
0785 %         elec.dims = row(1);
0786         elec.shape = 'R';
0787         elec.dims  = [row(1),hig];
0788      end
0789   else
0790      if length(row)<2 || row(2) == 0 % Circular electrodes
0791         elec.shape = 'C';
0792         elec.dims  = row(1);
0793      elseif row(2) == -1 % Point electrodes
0794         % Create rectangular electrodes with bottom, cw point where we want
0795         elec.shape = 'P'; 
0796         if length(row)>=3 && row(3) > 0
0797            elec.dims  =  row(3);
0798         else
0799            elec.dims  =  rad; % Make big if unspecified
0800         end
0801      elseif row(2)>0      % Rectangular electrodes
0802         elec.shape = 'R';
0803         elec.dims  = row(1:2);
0804      else
0805         error('negative electrode width');
0806      end
0807   end
0808 
0809   if length(row)>=3 && row(3) > 0
0810      elec.maxh = sprintf('-maxh=%f', row(3));
0811   else
0812      elec.maxh = '';
0813   end
0814 
0815   if length(row)<4 || row(4) == 0
0816      elec.model = 'cem'; % Complete Electrode Model (CEM)
0817   else
0818      elec.model = 'pem'; % Point Electrode Model (PEM)
0819   end
0820   %TODO support Shunt Electrode Model (SEM)
0821 
0822   if length(row) < 5 || row(5) == 0
0823       elec.discretize = 0;
0824   else
0825       elec.discretize = row(5);
0826   end
0827   
0828   
0829   
0830   
0831   
0832   
0833   
0834   
0835   
0836   
0837   
0838   
0839   
0840 function write_geo_file(geofn, tank_height, tank_shape, ...
0841                         tank_maxh, elecs, extra_ng_code)
0842     fid=fopen(geofn,'w');
0843     write_header(fid,tank_height,tank_shape,tank_maxh,extra_ng_code);
0844 
0845     n_verts = size(tank_shape.vertices,1);
0846     n_elecs = length(elecs);
0847     %  elecs(i).pos   = [x,y,z]
0848     %  elecs(i).shape = 'C' or 'R'
0849     %  elecs(i).dims  = [radius] or [width,height]
0850     %  elecs(i).maxh  = '-maxh=#' or '';
0851     %  elecs(i).edg_no = i (index of the edge on which the electrode lies)
0852     pts_elecs_idx = [];
0853     for i=1:n_elecs
0854         name = sprintf('elec%04d',i);
0855         pos = elecs(i).pos;
0856         dirn = elecs(i).normal;
0857         switch elecs(i).shape
0858             case 'C'
0859                 write_circ_elec(fid,name, pos, dirn,  ...
0860                     elecs(i).dims, tank_shape.centroid, elecs(i).maxh);
0861             case 'R'
0862                 write_rect_elec(fid,name, pos, dirn,  ...
0863                     elecs(i).dims, tank_shape.size/10, elecs(i).maxh);
0864                 %        case 'P'
0865                 %          pts_elecs_idx = [ pts_elecs_idx, i];
0866                 %          continue; % DON'T print solid cyl
0867 
0868             otherwise; error('unknown electrode shape');
0869         end
0870         %       fprintf(fid,'solid cyl%04d = trunk   and %s; \n',i,name);
0871     end
0872     fprintf(fid,'solid trunk = bound');
0873     if isfield(tank_shape,'additional_shapes')
0874          for i = 1:length(tank_shape.additional_shapes)
0875              fprintf(fid,' and not add_obj%04d',i);
0876          end
0877     end
0878     fprintf(fid,';\n');
0879     
0880     if isfield(tank_shape,'additional_shapes')
0881         for i = 1:length(tank_shape.additional_shapes)
0882             fprintf(fid,'solid add_obj%04dc = add_obj%04d',i,i);
0883             for j = (i+1):length(tank_shape.additional_shapes)
0884                 fprintf(fid,' and not add_obj%04d',j);
0885             end
0886 
0887 % This code was added while trying to debug mixed shapes
0888 %   with solid geometry and extruded shapes. It didn't help
0889 %           if ~isempty(extra_ng_code{1})
0890 %                fprintf(fid,' and not %s',extra_ng_code{1});
0891 %           end
0892 
0893             fprintf(fid,[' and plane(0,0,0;0,0,-1)\n' ...
0894                 '      and  plane(0,0,%6.2f;0,0,1)'],tank_height);
0895             fprintf(fid,';\n');
0896         end
0897     end
0898     
0899     if tank_maxh(1) ~= 0
0900         fprintf(fid,'tlo trunk -transparent -maxh=%f;\n',tank_maxh);
0901     else
0902         fprintf(fid,'tlo trunk -transparent;\n');
0903     end
0904     if ~isempty(extra_ng_code{1})
0905         fprintf(fid,'tlo %s -col=[0,1,0];\n',extra_ng_code{1});
0906     end
0907 
0908     if isfield(tank_shape,'additional_shapes')
0909          for i = 1:length(tank_shape.additional_shapes)
0910              fprintf(fid,'tlo add_obj%04dc -col=[0,1,0];\n',i);
0911          end
0912     end
0913 
0914     for i=1:n_elecs
0915         if any(i == pts_elecs_idx); continue; end
0916         fprintf(fid,'tlo elec%04d -col=[1,0,0] %s;\n',i,elecs(i).maxh);
0917     end
0918 
0919 
0920     fclose(fid); % geofn
0921 
0922    
0923    
0924    function write_header(fid,tank_height,tank_shape,maxh,extra)
0925    if maxh(1)==0; 
0926       maxsz = '';
0927    else
0928       maxsz = sprintf('-maxh=%f',maxh);
0929    end
0930 
0931    if ~isempty( extra{1} )
0932       extra{1} = [' and not ',extra{1}];
0933    end
0934 
0935    
0936    fprintf(fid,'#Automatically generated by ng_mk_extruded_model\n');
0937    fprintf(fid,'algebraic3d\n');
0938    fprintf(fid,'%s\n',extra{2}); % Define extra stuff here
0939    
0940    fprintf(fid,'curve3d extrsncurve=(2; 0,0,0; 0,0,%6.2f; 1; 2,1,2);\n', ...
0941        tank_height+1);
0942 
0943 
0944    write_curve(fid,tank_shape,'outer', 1.15);
0945    write_curve(fid,tank_shape,'inner', 0.99);
0946    write_curve(fid,tank_shape,'surf', 1);
0947    
0948     fprintf(fid,['solid bound= plane(0,0,0;0,0,-1)\n' ...
0949                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0950                 '      and  extrusion(extrsncurve;surf;0,1,0)'...
0951                 '%s %s;\n'],tank_height,extra{1},maxsz);
0952             
0953    fprintf(fid,['solid inner_bound= plane(0,0,0;0,0,-1)\n' ...
0954                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0955                 '      and  extrusion(extrsncurve;inner;0,1,0)'...
0956                 '%s %s;\n'],tank_height,extra{1},maxsz);
0957 
0958    fprintf(fid,['solid outer_bound= plane(0,0,0;0,0,-1)\n' ...
0959                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0960                 '      and  extrusion(extrsncurve;outer;0,1,0)'...
0961                 '%s %s;\n'],tank_height,extra{1},maxsz);
0962            
0963    % EVERYTHING below this line assumes additional shapes are defined
0964    if ~isfield(tank_shape, 'additional_shapes'), return, end
0965    
0966    for i = 1:length(tank_shape.additional_shapes)
0967        name_curve = sprintf('add_curve%04d',i); 
0968        write_curve(fid,tank_shape.additional_shapes{i},name_curve);
0969        name_obj = sprintf('add_obj%04d',i); 
0970        if maxh(1+i)==0
0971            maxsz = '';
0972        else
0973            maxsz = sprintf('-maxh=%f',maxh(1+i));
0974        end
0975        
0976        fprintf(fid,['solid %s= plane(0,0,%6.2f;0,0,-1)\n' ...
0977            '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0978            '      and  extrusion(extrsncurve;%s;0,1,0)'...
0979            '%s %s;\n'],name_obj,-i,tank_height+i,name_curve,extra{1},maxsz);
0980    end
0981                    
0982         
0983    function write_curve(fid, tank_shape, name, scale)
0984         if nargin <4
0985             scale = 1;
0986         end
0987        
0988         is_struct = isstruct(tank_shape);
0989         if ~is_struct
0990             vertices = tank_shape;
0991             STRUCT = false;
0992             if scale ~= 1
0993                 warning('Scale is ignored when second input is an array');
0994                 scale = 1;
0995             end
0996         elseif scale ~= 1
0997             vertices = tank_shape.vertices + ...
0998                 (scale-1)*tank_shape.vertex_dir*tank_shape.size;
0999         else
1000             vertices = tank_shape.vertices;
1001         end
1002        n_vert = size(vertices,1);
1003        
1004        fprintf(fid,'curve2d %s=(%d; \n', name, n_vert);
1005        
1006        for i = 1:n_vert
1007            % because of the definitions of the local axis in extrusion, the
1008            % x coordinate has to be multiplied by -1. This assures the
1009            % object appears at the expected coordinates. To maintain
1010            % clockwise order (required by netget) the vertices are printed
1011            % in the opposite order.
1012            fprintf(fid,'       %6.4f, %6.4f;\n',[-1 1].*vertices(n_vert-i+1,:));
1013            %             fprintf(fid,'       %6.2f, %6.2f;\n',vertices(i,:));
1014        end
1015        if is_struct
1016            spln_sgmnts = tank_shape.spln_sgmnts;
1017        else
1018            spln_sgmnts = zeros(max(size(vertices)));
1019        end
1020        n_sgmnts = length(spln_sgmnts);
1021        fprintf(fid,'       %d;\n',n_sgmnts);
1022        cv = 1; %current vertex
1023        for i = 1:n_sgmnts
1024            if spln_sgmnts(i)
1025                if i == n_sgmnts
1026                   fprintf(fid,'       %d, %d, %d, %d );\n\n\n', 3, cv,cv+1, 1);
1027                else
1028                    fprintf(fid,'       %d, %d, %d, %d; \n', 3, cv, cv+1, cv+2);
1029                end
1030                cv = cv + 2;
1031            else
1032                if i == n_sgmnts
1033                    fprintf(fid,'       %d, %d, %d );\n\n\n', 2, cv, 1);
1034                else
1035                    fprintf(fid,'       %d, %d, %d; \n', 2, cv, cv+1);
1036                end
1037                cv = cv + 1;
1038            end
1039        end
1040        
1041        
1042 function write_circ_elec(fid,name,c, dirn, rd, centroid, maxh)
1043 % writes the specification for a netgen cylindrical rod on fid,
1044 %  named name, centerd on c,
1045 % in the direction given by vector dirn, radius rd
1046 % direction is in the xy plane
1047 
1048     % the direction vector
1049     dirn(3) = 0; dirn = dirn/norm(dirn);
1050 
1051     fprintf(fid,'solid %s  = ', name);
1052     fprintf(fid,['  outer_bound and not inner_bound and '...
1053         'cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) '...
1054         'and plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) '...
1055         'and not bound;\n'], ...
1056         c(1)-dirn(1),c(2)-dirn(2),c(3)-dirn(3),...
1057         c(1)+dirn(1),c(2)+dirn(2),c(3)+dirn(3), rd, ...
1058         centroid(1), centroid(2), 0, -dirn(1), -dirn(2), dirn(3));
1059 
1060 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
1061 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
1062 % in the direction given by vector dirn,
1063 % hw = [height, width]  and depth d
1064 % direction is in the xy plane
1065    w = wh(1); h= wh(2);
1066    dirn(3) = 0; dirn = dirn/norm(dirn);
1067    dirnp = [-dirn(2),dirn(1),0];
1068    dirnp = dirnp/norm(dirnp);
1069 
1070    bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
1071    tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
1072    fprintf(fid,'solid %s  = outer_bound and not inner_bound and', name);
1073    fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
1074            bl(1),bl(2),bl(3));
1075    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1076            bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
1077    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1078            bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
1079    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
1080            tr(1),tr(2),tr(3));
1081    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1082            tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
1083    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f  )\n', ...
1084            tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0);
1085    fprintf(fid,' and not bound;\n');
1086     
1087 % NG_REMOVE_ELECTRODES: cleans up matrices read from a *.vol file
1088 % [srf,vtx,fc,bc,simp,edg,mat_ind]= ng_remove_electrodes...
1089 %     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1090 %
1091 % Used to clean up external objects used to force electrode meshing in
1092 % ng_mk_extruded_model.
1093 %
1094 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes...
1095     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1096 
1097 fc = []; % Unused, and we're not sure what it is;
1098 
1099 % total objects:
1100 N_obj = max(mat_ind);
1101 
1102 % The electodes are the last N_elec simps
1103 elec_ind = mat_ind > (N_obj - N_elec);
1104 
1105 in = unique(simp(~elec_ind,:)); % nodes in real object
1106 out = unique(simp(elec_ind,:)); % nodes in electrodes
1107 boundary = intersect(in,out);   % nodes shared obj/electrodes
1108 out = setdiff(out,boundary);    % nodes only in electrodes
1109 
1110 % remove simps which contain nodes in the "out" list
1111 remove_simp = any( ismember(simp,out), 2);
1112 simp0 = simp;
1113 simp( remove_simp,:) = [];
1114 
1115 % Choose which vertices to keep
1116 vtx_renum = logical( zeros(size(vtx,1),1) );
1117 vtx_renum( in ) = logical(1);
1118 vtx_renum = cumsum(vtx_renum);
1119 
1120 vtx(out,:) = [];
1121 simp =  reshape( vtx_renum(simp), size(simp));
1122 
1123 % recalculate surface
1124 % STUPID MATLAB BUGS MEAN WE CANT allow int32 here
1125 srf= double( find_boundary(simp) );
1126 bc = ones(size(srf,1),1); % Add srf for the electrodes
1127 
1128 % Iterate over electrodes
1129 for i=1:N_elec;
1130   eleci_obj = mat_ind == (N_obj - N_elec + i);
1131   this_elec = unique( simp0( eleci_obj, : ));
1132   eleci_nodes = vtx_renum( intersect( this_elec, in )); 
1133 
1134 % This is the direct way to get electrodes. Instead we need to call the
1135 %   electrode finder function
1136 % elec(i).nodes = eleci_nodes;
1137   
1138   eleci_srf = all( ismember(srf, eleci_nodes), 2);
1139   bc( eleci_srf ) = i+1; % give this elec a surface
1140 end
1141 
1142 mat_ind( remove_simp) = [];
1143 
1144 % Test code:
1145 % fmdl.type='fwd_model'; fmdl.nodes = vtx; fmdl.elems =  simp_obj; fmdl.electrode= elec;
1146 
1147 
1148 
1149 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes_old...
1150     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1151 
1152 % total objects:
1153 N_obj = max(mat_ind);
1154 
1155 % The electodes are the last N_elec simps
1156 e_simp_ind = mat_ind > (N_obj - N_elec);
1157 
1158 in = unique(simp(~e_simp_ind,:));
1159 out = unique(simp(e_simp_ind,:));
1160 boundary = intersect(in,out);
1161 out = setdiff(out,boundary);
1162 
1163 ext_srf_ind = ismember(srf,out);
1164 ext_srf_ind = ext_srf_ind(:,1) | ext_srf_ind(:,2) | ext_srf_ind(:,3);
1165 
1166 srf(ext_srf_ind,:) = [];
1167 bc(ext_srf_ind,:) = [];
1168 fc(ext_srf_ind,:) = [];
1169 simp = simp(~e_simp_ind,:);
1170 mat_ind = mat_ind(~e_simp_ind);
1171 
1172 % fix bc:
1173 n_unique = numel(unique(bc));
1174 missing = setdiff(1:n_unique, unique(bc));
1175 spare = setdiff(unique(bc), 1:n_unique); 
1176 
1177 for i = 1:length(missing)
1178     bc( bc==spare(i) ) = missing(i);
1179 end
1180 
1181 % fix vtx:
1182 v = 1:size(vtx,1);
1183 unused_v = setdiff(v, union(unique(simp),unique(srf))); 
1184 v(unused_v) = [];
1185 for i = 1:size(vtx,1);
1186 %     simp_ind = find(simp == i);
1187 %     srf_ind = find( srf == i);
1188     new_v_ind = find(v == i);
1189     simp( simp == i ) = new_v_ind; 
1190     srf( srf  == i ) = new_v_ind;
1191 end
1192 vtx(unused_v,:) = [];
1193 
1194 
1195 %%
1196 function do_unit_test
1197    %do_unit_test_basic;
1198  do_unit_test_human_thorax;
1199 
1200 
1201 function do_unit_test_basic;
1202     fmdl = [];
1203     mat_idx = [];
1204         a = [
1205        -0.8981   -0.7492   -0.2146    0.3162    0.7935    0.9615    0.6751    0.0565   -0.3635   -0.9745
1206         0.1404    0.5146    0.3504    0.5069    0.2702   -0.2339   -0.8677   -0.6997   -0.8563   -0.4668 ]';
1207     % [fmdl, mat_idx] = ng_mk_extruded_model({2,{a,0.5*a,0.2*a},1},[16,0,1],[0.01]);
1208      
1209     %%
1210     xx=[
1211       -88.5777  -11.4887    4.6893   49.8963  122.7033  150.3033  195.5103 249.7573 ...
1212       258.8013  279.7393  304.9623  309.2443  322.0923  337.7963  340.6503 348.2633 ...
1213       357.3043  358.7333  361.5873  364.9183  365.3943  366.3453  366.3453 365.3943 ...
1214       362.5393  351.5943  343.5053  326.8513  299.2503  288.3073  264.9923 224.0703 ...
1215       206.4633  162.6833  106.5313   92.2543   57.5153    7.0733   -8.6297 -42.4167 ...
1216       -90.9547 -105.7057 -134.2577 -178.0367 -193.2647 -222.7687 -265.5957 -278.9197 ...
1217      -313.1817 -355.5337 -363.6237 -379.3267 -397.8857 -400.7407 -401.6927 -398.8377 ...
1218      -395.0307 -384.0867 -368.3837 -363.6247 -351.7277 -334.1217 -328.4117 -314.1357 ...
1219      -291.2947 -282.7297 -267.0257 -236.5707 -221.8187 -196.5977 -159.4807 -147.5837];
1220 
1221     yy=[
1222      -385.8513 -386.8033 -386.3273 -384.8993 -368.7193 -353.9673 -323.0363 -283.5403 ...
1223      -274.9743 -254.0363 -225.4843 -217.8703 -187.4153 -140.7813 -124.6013  -86.0573 ...
1224       -38.4703  -29.4273   -9.9173   21.0137   32.4347   53.3727   83.8257   93.3437 ...
1225       114.7587  149.0237  161.8717  187.5677  222.3037  231.3447  247.5237  267.5087 ...
1226       271.3177  277.0297  281.3127  279.4097  274.6507  273.2227  276.5547  284.6447 ...
1227       295.1127  297.4927  301.7757  304.1557  302.2537  297.4947  287.5017  282.2667 ...
1228       259.9017  225.6387  213.7427  185.6677  141.4127  125.2337   88.5917   34.8187 ...
1229        17.6897  -22.2803  -73.6723  -85.0923 -117.9263 -163.6083 -176.4573 -205.9613 ...
1230      -245.9343 -256.4023 -275.4373 -304.9403 -315.4083 -332.0623 -352.0473 -355.3783];
1231 
1232     a = [xx; yy]';
1233     a = flipud(a);
1234     % th=linspace(0,2*pi,33)'; th(end)=[];
1235     % a=[sin(th)*0.3,cos(th)];
1236 
1237 
1238     ng_write_opt('MSZBRICK',[-400,400,-400,400,120,180,40]);  
1239     fmdl = ng_mk_extruded_model({300,a,[4,25]},[16,1.11,150],[1]);
1240     show_fem(fmdl);
1241 
1242 function do_unit_test_human_thorax
1243 % load CT2
1244 
1245 %[fmdl, mat_idx] = ng_mk_extruded_model({150,flipud(trunk),1},[16,0,75],[0.01]);
1246 
1247 %[fmdl, mat_idx] = ng_mk_extruded_model({2,{trunk/100, lung_heart_dep/100, heart/100},1},[16,1,1],[0.1]);
1248 %img = mk_image( fmdl, 1);
1249 % img.elem_data(mat_idx{2}) = 1.1;
1250 
1251  trunk = [    -4    -2     2     4     4     2    -2    -4
1252                2     4     4     2    -2    -4    -4    -2]';
1253  heart_lung = [    -2    -1    -0.8  0.8  1     2     2    -2
1254                     1     2     1.8  1.8  2     1    -2    -2]';    
1255  lung = [    -2    -1    -1  -1  1     2     2    -2
1256              1     2     0   0  2     1    -2    -2]';   
1257  heart = [    -1    -1     1     1
1258                0     2     2     0]';
1259 
1260  [fmdl, mat_idx] = ng_mk_extruded_model({2,{trunk, heart_lung, heart},1},[16,1,1],[0.1]);
1261 
1262  
1263   figure, show_fem( fmdl );
ng_mk_extruded_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
