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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
