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 6033 2019-12-30 16:50:55Z 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 electrodes
0565         [th,~] = cart2pol(normal(1),normal(2)); % ~ needed for octave
0566         frac = 2*pi /elecs(i).discretize ;
0567         th = frac * round( th / frac);
0568         [normal(1) normal(2)] = pol2cart(th,1);
0569        end
0570        elecs(i).normal = normal;
0571        
0572    end
0573 
0574    if n_elecs == 0
0575       elecs= struct([]); % empty
0576       centres= []; 
0577    end
0578 
0579    
0580    
0581     function [pos, normal] = calc_elec_centre(tank_shape, th)
0582         % The calculation relies on the theorem that if point D lies on a
0583         % line between B and C, but point A is not on that line, then:
0584         %   |BD|    |AB| sin(<DAB)
0585         %   ---- = ---------------
0586         %   |DC|    |AC| sin(<DAC)
0587         % Thus, B and C are vertices of our shape, A is its centroid and D
0588         % is the sought center of the electrode. All quantities on RHS are
0589         % known.
0590         
0591         % make sure th is between -pi and pi
0592         if th > pi; th = th - 2*pi; end
0593         
0594         
0595         vert_pol = tank_shape.vertices_polar; %[th, r]
0596         
0597      
0598         n_vert = size(vert_pol,1);
0599         vert_pol = [vert_pol , (1:n_vert)'];
0600         % Re-order the vertices -pi to pi. Now counter-clockwise
0601         vert_pol = sortrows(vert_pol,1); 
0602         % find the edge on which the elctrode lies. (Edge 1 is between
0603         % verticies 1 and 2)
0604         idx = find(vert_pol(:,1) > th, 1, 'first');
0605         if isempty(idx); idx = 1; end
0606         edg_no = vert_pol(idx,3);
0607         
0608         
0609         normal = tank_shape.edge_normals(edg_no,:);
0610               
0611         v1 = edg_no;
0612         if edg_no == n_vert % between the last and first vertex
0613             v2 = 1;
0614         else
0615             v2 = v1+1;
0616         end
0617         vert_pol = [];
0618         
0619         
0620         vert_pol = tank_shape.vertices_polar;
0621         vert = tank_shape.vertices;
0622         cntr = tank_shape.centroid;
0623         % position between vertices - see first comment
0624         AB = sqrt(sum( (vert(v1,:) - cntr).^2 ));
0625         AC = sqrt(sum( (vert(v2,:) - cntr).^2 ));
0626         DAB = abs(vert_pol(v1,1)-th); 
0627         if DAB > pi, DAB = abs( DAB - 2*pi); end; 
0628         DAC  = abs(vert_pol(v2,1)-th);
0629         if DAC > pi, DAC = abs( DAC - 2*pi); end;
0630         if DAC ~= 0
0631             ratio = AB * sin(DAB) / (AC * sin(DAC));
0632             pos = vert(v1,:) + ( ratio / (1 + ratio) ) * (vert(v2,:) - vert(v1,:));
0633         else
0634             pos = vert(v2,:);
0635         end
0636 
0637         
0638         
0639    function [pos, normal] = calc_elec_centre_spline(tank_shape, th)
0640         % The calculation proceeds by finding a common point between a line
0641         % from the centroid outwards and the equation of the relevant
0642         % quadratic spline segment defined using 3 control points
0643         
0644         % make sure th is between -pi and pi
0645         if th > pi; th = th - 2*pi; end 
0646         
0647         vert_pol = tank_shape.vertices_polar; %[th, r]
0648         
0649         % The number of vertices must be even, but just in case...
0650         if mod(size(vert_pol,1),2)
0651             error(['The number of points must be even. '...
0652                 'One de Boor control point for every vertex']);
0653         end
0654         
0655         % if the curve is defined as splines, every second point is not
0656         % actually a vertex. We remove them.
0657         if tank_shape.curve_type == 2 || tank_shape.curve_type == 3
0658             vert_pol(2:2:end,:) = [];
0659         end
0660       
0661         n_vert = size(vert_pol,1);
0662    
0663         vert_pol = [vert_pol , (1:n_vert)']; %excludes control points
0664         % Re-order the vertices -pi to pi. Now counter-clockwise
0665         vert_pol = sortrows(vert_pol,1); 
0666         % find the edge on which the electrode lies. Edge 1 is between
0667         % vertices 1 and 2.
0668         idx = find(vert_pol(:,1) > th, 1, 'first');
0669         if isempty(idx); idx = 1; end
0670         edg_no = vert_pol(idx,3);
0671         
0672         v1 = edg_no;
0673         if edg_no == n_vert % between the last and first vertex
0674             v2 = 1;
0675         else
0676             v2 = v1+1;
0677         end
0678         vert_pol = [];
0679         
0680         % correcting for the control points
0681         v1 = 2 * v1 - 1;
0682         v2 = 2 * v2 - 1;
0683         
0684         % the spline goes from point P0 to P2 such that P1-P0 is a tangent
0685         % at P0 and P2-P1 is a tangent at P2
0686         C = tank_shape.centroid;
0687         P0 = tank_shape.vertices(v1,:) - C;
0688         P1 = tank_shape.vertices(v1+1,:) - C; % control point
0689         P2 = tank_shape.vertices(v2,:) - C;
0690         
0691         
0692         % find the gradient of the line from centroid to electrode center:
0693         [x, y] = pol2cart(th, 1);
0694         % FIXME: This doesn't crash only because of round-off errors.
0695         g = y/x;
0696         % (because we subtracted the centroid from the vertices, the line
0697         % passes through the origin now)
0698         
0699         % the spline is f(t) = (1-t)^2 * P0 + 2t(1-t)P1 + t^2 * P2
0700         % which can also be expressed as
0701         f = @(t) (P2 - 2*P1 + P0)*t^2 + 2*(P1 - P0)*t + P0;
0702         % and it's derivative:
0703         df = @(t) 2*(P2 - 2*P1 + P0)*t + 2*(P1 - P0);
0704         % to find the value of t for which the line cross, we substitute
0705         % p0 = y0-ax0 for P0 and so on.
0706         p0 = P0(2) - g * P0(1);
0707         p1 = P1(2) - g * P1(1);
0708         p2 = P2(2) - g * P2(1);
0709         
0710         % thus we have a quadratic equation a*t^2 + b*t + c = 0 where
0711         a = (p2 - 2*p1 + p0);
0712         b = 2* (p1 - p0);
0713         c = p0;
0714         
0715         if abs(a) < 1e-10
0716             t = -c/b;
0717             pos = f(t) + C;
0718             tmp = df(t);
0719             normal = [-tmp(2), tmp(1)] / sqrt(sum(tmp.^2));
0720             return;
0721         end
0722         
0723         % the determinant is
0724         D = b^2 - 4*a*c;
0725         
0726         % find the roots
0727         if D == 0
0728             t = -b / (2 * a);
0729 
0730         elseif D > 0
0731             t1 = (-b - sqrt(D) ) / (2 * a);
0732             t2 = (-b + sqrt(D) ) / (2 * a);
0733             if t1 >= 0 && t1 <= 1
0734                 t = t1;
0735             else
0736                 t = t2;
0737             end
0738         else
0739             error('Something went wrong, cannot place electrode on spline');
0740         end
0741         
0742         pos = f(t) + C;
0743         tmp = df(t);
0744         normal = [-tmp(2), tmp(1)]/ sqrt(sum(tmp.^2));
0745 
0746    
0747    
0748 
0749 function elec = elec_spec( row, is2D, hig, rad )
0750   if     is2D
0751      if length(row)>=2 && row(2) == -1 % Point electrodes
0752         % Create rectangular electrodes with bottom, cw point where we want
0753         elec.shape = 'P' ;
0754         if length(row)>=3 && row(3) > 0
0755            elec.dims  =  row(3);
0756         else
0757            elec.dims  =  rad; % Make big if unspecified
0758         end
0759      else
0760         % create circular electrodes for now, rectangular not yet supported
0761 %         elec.shape = 'C';
0762 %         elec.dims = row(1);
0763         elec.shape = 'R';
0764         elec.dims  = [row(1),hig];
0765      end
0766   else
0767      if length(row)<2 || row(2) == 0 % Circular electrodes
0768         elec.shape = 'C';
0769         elec.dims  = row(1);
0770      elseif row(2) == -1 % Point electrodes
0771         % Create rectangular electrodes with bottom, cw point where we want
0772         elec.shape = 'P'; 
0773         if length(row)>=3 && row(3) > 0
0774            elec.dims  =  row(3);
0775         else
0776            elec.dims  =  rad; % Make big if unspecified
0777         end
0778      elseif row(2)>0      % Rectangular electrodes
0779         elec.shape = 'R';
0780         elec.dims  = row(1:2);
0781      else
0782         error('negative electrode width');
0783      end
0784   end
0785 
0786   if length(row)>=3 && row(3) > 0
0787      elec.maxh = sprintf('-maxh=%f', row(3));
0788   else
0789      elec.maxh = '';
0790   end
0791 
0792   if length(row)<4 || row(4) == 0
0793      elec.model = 'cem'; % Complete Electrode Model (CEM)
0794   else
0795      elec.model = 'pem'; % Point Electrode Model (PEM)
0796   end
0797   %TODO support Shunt Electrode Model (SEM)
0798 
0799   if length(row) < 5 || row(5) == 0
0800       elec.discretize = 0;
0801   else
0802       elec.discretize = row(5);
0803   end
0804   
0805   
0806   
0807   
0808   
0809   
0810   
0811   
0812   
0813   
0814   
0815   
0816   
0817 function write_geo_file(geofn, tank_height, tank_shape, ...
0818                         tank_maxh, elecs, extra_ng_code)
0819     fid=fopen(geofn,'w');
0820     write_header(fid,tank_height,tank_shape,tank_maxh,extra_ng_code);
0821 
0822     n_verts = size(tank_shape.vertices,1);
0823     n_elecs = length(elecs);
0824     %  elecs(i).pos   = [x,y,z]
0825     %  elecs(i).shape = 'C' or 'R'
0826     %  elecs(i).dims  = [radius] or [width,height]
0827     %  elecs(i).maxh  = '-maxh=#' or '';
0828     %  elecs(i).edg_no = i (index of the edge on which the electrode lies)
0829     pts_elecs_idx = [];
0830     %^keyboard
0831     for i=1:n_elecs
0832         name = sprintf('elec%04d',i);
0833         pos = elecs(i).pos;
0834         dirn = elecs(i).normal;
0835         switch elecs(i).shape
0836             case 'C'
0837                 write_circ_elec(fid,name, pos, dirn,  ...
0838                     elecs(i).dims, tank_shape.centroid, elecs(i).maxh);
0839             case 'R'
0840                 write_rect_elec(fid,name, pos, dirn,  ...
0841                     elecs(i).dims, tank_shape.size/10, elecs(i).maxh);
0842                 %        case 'P'
0843                 %          pts_elecs_idx = [ pts_elecs_idx, i];
0844                 %          continue; % DON'T print solid cyl
0845 
0846             otherwise; error('unknown electrode shape');
0847         end
0848         %       fprintf(fid,'solid cyl%04d = trunk   and %s; \n',i,name);
0849     end
0850     fprintf(fid,'solid trunk = bound');
0851     if isfield(tank_shape,'additional_shapes')
0852          for i = 1:length(tank_shape.additional_shapes)
0853              fprintf(fid,' and not add_obj%04d',i);
0854          end
0855     end
0856     fprintf(fid,';\n');
0857     
0858     if isfield(tank_shape,'additional_shapes')
0859         for i = 1:length(tank_shape.additional_shapes)
0860             fprintf(fid,'solid add_obj%04dc = add_obj%04d',i,i);
0861             for j = (i+1):length(tank_shape.additional_shapes)
0862                 fprintf(fid,' and not add_obj%04d',j);
0863             end
0864 
0865 % This code was added while trying to debug mixed shapes
0866 %   with solid geometry and extruded shapes. It didn't help
0867 %           if ~isempty(extra_ng_code{1})
0868 %                fprintf(fid,' and not %s',extra_ng_code{1});
0869 %           end
0870 
0871             fprintf(fid,[' and plane(0,0,0;0,0,-1)\n' ...
0872                 '      and  plane(0,0,%6.2f;0,0,1)'],tank_height);
0873             fprintf(fid,';\n');
0874         end
0875     end
0876     
0877     if tank_maxh ~= 0
0878         fprintf(fid,'tlo trunk -transparent -maxh=%f;\n',tank_maxh);
0879     else
0880         fprintf(fid,'tlo trunk -transparent;\n');
0881     end
0882     if ~isempty(extra_ng_code{1})
0883         fprintf(fid,'tlo %s -col=[0,1,0];\n',extra_ng_code{1});
0884     end
0885 
0886     if isfield(tank_shape,'additional_shapes')
0887          for i = 1:length(tank_shape.additional_shapes)
0888              fprintf(fid,'tlo add_obj%04dc -col=[0,1,0];\n',i);
0889          end
0890     end
0891 
0892     for i=1:n_elecs
0893         if any(i == pts_elecs_idx); continue; end
0894         fprintf(fid,'tlo elec%04d -col=[1,0,0] %s;\n',i,elecs(i).maxh);
0895     end
0896 
0897 
0898     fclose(fid); % geofn
0899 
0900    
0901    
0902    function write_header(fid,tank_height,tank_shape,maxsz,extra)
0903    if maxsz==0; 
0904       maxsz = '';
0905    else
0906       maxsz = sprintf('-maxh=%f',maxsz);
0907    end
0908 
0909    if ~isempty( extra{1} )
0910       extra{1} = [' and not ',extra{1}];
0911    end
0912 
0913    
0914    fprintf(fid,'#Automatically generated by ng_mk_extruded_model\n');
0915    fprintf(fid,'algebraic3d\n');
0916    fprintf(fid,'%s\n',extra{2}); % Define extra stuff here
0917    
0918    fprintf(fid,'curve3d extrsncurve=(2; 0,0,0; 0,0,%6.2f; 1; 2,1,2);\n', ...
0919        tank_height+1);
0920 
0921 
0922    write_curve(fid,tank_shape,'outer', 1.15);
0923    write_curve(fid,tank_shape,'inner', 0.99);
0924    write_curve(fid,tank_shape,'surf', 1);
0925    
0926     fprintf(fid,['solid bound= plane(0,0,0;0,0,-1)\n' ...
0927                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0928                 '      and  extrusion(extrsncurve;surf;0,1,0)'...
0929                 '%s %s;\n'],tank_height,extra{1},maxsz);
0930             
0931    fprintf(fid,['solid inner_bound= plane(0,0,0;0,0,-1)\n' ...
0932                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0933                 '      and  extrusion(extrsncurve;inner;0,1,0)'...
0934                 '%s %s;\n'],tank_height,extra{1},maxsz);
0935 
0936    fprintf(fid,['solid outer_bound= plane(0,0,0;0,0,-1)\n' ...
0937                 '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0938                 '      and  extrusion(extrsncurve;outer;0,1,0)'...
0939                 '%s %s;\n'],tank_height,extra{1},maxsz);
0940            
0941    % EVERYTHING below this line assumes additional shapes are defined
0942    if ~isfield(tank_shape, 'additional_shapes'), return, end
0943    
0944    for i = 1:length(tank_shape.additional_shapes)
0945        name_curve = sprintf('add_curve%04d',i); 
0946        write_curve(fid,tank_shape.additional_shapes{i},name_curve);
0947        name_obj = sprintf('add_obj%04d',i); 
0948        fprintf(fid,['solid %s= plane(0,0,%6.2f;0,0,-1)\n' ...
0949            '      and  plane(0,0,%6.2f;0,0,1)\n' ...
0950            '      and  extrusion(extrsncurve;%s;0,1,0)'...
0951            '%s %s;\n'],name_obj,-i,tank_height+i,name_curve,extra{1},maxsz);
0952    end
0953                    
0954         
0955    function write_curve(fid, tank_shape, name, scale)
0956         if nargin <4
0957             scale = 1;
0958         end
0959        
0960         is_struct = isstruct(tank_shape);
0961         if ~is_struct
0962             vertices = tank_shape;
0963             STRUCT = false;
0964             if scale ~= 1
0965                 warning('Scale is ignored when second input is an array');
0966                 scale = 1;
0967             end
0968         elseif scale ~= 1
0969             vertices = tank_shape.vertices + ...
0970                 (scale-1)*tank_shape.vertex_dir*tank_shape.size;
0971         else
0972             vertices = tank_shape.vertices;
0973         end
0974        n_vert = size(vertices,1);
0975        
0976        fprintf(fid,'curve2d %s=(%d; \n', name, n_vert);
0977        
0978        for i = 1:n_vert
0979            % because of the definitions of the local axis in extrusion, the
0980            % x coordinate has to be multiplied by -1. This assures the
0981            % object appears at the expected coordinates. To maintain
0982            % clockwise order (required by netget) the vertices are printed
0983            % in the opposite order.
0984            fprintf(fid,'       %6.4f, %6.4f;\n',[-1 1].*vertices(n_vert-i+1,:));
0985            %             fprintf(fid,'       %6.2f, %6.2f;\n',vertices(i,:));
0986        end
0987        if is_struct
0988            spln_sgmnts = tank_shape.spln_sgmnts;
0989        else
0990            spln_sgmnts = zeros(max(size(vertices)));
0991        end
0992        n_sgmnts = length(spln_sgmnts);
0993        fprintf(fid,'       %d;\n',n_sgmnts);
0994        cv = 1; %current vertex
0995        for i = 1:n_sgmnts
0996            if spln_sgmnts(i)
0997                if i == n_sgmnts
0998                   fprintf(fid,'       %d, %d, %d, %d );\n\n\n', 3, cv,cv+1, 1);
0999                else
1000                    fprintf(fid,'       %d, %d, %d, %d; \n', 3, cv, cv+1, cv+2);
1001                end
1002                cv = cv + 2;
1003            else
1004                if i == n_sgmnts
1005                    fprintf(fid,'       %d, %d, %d );\n\n\n', 2, cv, 1);
1006                else
1007                    fprintf(fid,'       %d, %d, %d; \n', 2, cv, cv+1);
1008                end
1009                cv = cv + 1;
1010            end
1011        end
1012        
1013        
1014 function write_circ_elec(fid,name,c, dirn, rd, centroid, maxh)
1015 % writes the specification for a netgen cylindrical rod on fid,
1016 %  named name, centerd on c,
1017 % in the direction given by vector dirn, radius rd
1018 % direction is in the xy plane
1019 
1020     % the direction vector
1021     dirn(3) = 0; dirn = dirn/norm(dirn);
1022 
1023     fprintf(fid,'solid %s  = ', name);
1024     fprintf(fid,['  outer_bound and not inner_bound and '...
1025         'cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) '...
1026         'and plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) '...
1027         'and not bound;\n'], ...
1028         c(1)-dirn(1),c(2)-dirn(2),c(3)-dirn(3),...
1029         c(1)+dirn(1),c(2)+dirn(2),c(3)+dirn(3), rd, ...
1030         centroid(1), centroid(2), 0, -dirn(1), -dirn(2), dirn(3));
1031 
1032 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
1033 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
1034 % in the direction given by vector dirn,
1035 % hw = [height, width]  and depth d
1036 % direction is in the xy plane
1037    w = wh(1); h= wh(2);
1038    dirn(3) = 0; dirn = dirn/norm(dirn);
1039    dirnp = [-dirn(2),dirn(1),0];
1040    dirnp = dirnp/norm(dirnp);
1041 
1042    bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
1043    tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
1044    fprintf(fid,'solid %s  = outer_bound and not inner_bound and', name);
1045    fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
1046            bl(1),bl(2),bl(3));
1047    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1048            bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
1049    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1050            bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
1051    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
1052            tr(1),tr(2),tr(3));
1053    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1054            tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
1055    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f  )\n', ...
1056            tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0);
1057    fprintf(fid,' and not bound;\n');
1058     
1059 % NG_REMOVE_ELECTRODES: cleans up matrices read from a *.vol file
1060 % [srf,vtx,fc,bc,simp,edg,mat_ind]= ng_remove_electrodes...
1061 %     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1062 %
1063 % Used to clean up external objects used to force electrode meshing in
1064 % ng_mk_extruded_model.
1065 %
1066 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes...
1067     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1068 
1069 fc = []; % Unused, and we're not sure what it is;
1070 
1071 % total objects:
1072 N_obj = max(mat_ind);
1073 
1074 % The electodes are the last N_elec simps
1075 elec_ind = mat_ind > (N_obj - N_elec);
1076 
1077 in = unique(simp(~elec_ind,:)); % nodes in real object
1078 out = unique(simp(elec_ind,:)); % nodes in electrodes
1079 boundary = intersect(in,out);   % nodes shared obj/electrodes
1080 out = setdiff(out,boundary);    % nodes only in electrodes
1081 
1082 % remove simps which contain nodes in the "out" list
1083 remove_simp = any( ismember(simp,out), 2);
1084 simp0 = simp;
1085 simp( remove_simp,:) = [];
1086 
1087 % Choose which vertices to keep
1088 vtx_renum = logical( zeros(size(vtx,1),1) );
1089 vtx_renum( in ) = logical(1);
1090 vtx_renum = cumsum(vtx_renum);
1091 
1092 vtx(out,:) = [];
1093 simp =  reshape( vtx_renum(simp), size(simp));
1094 
1095 % recalculate surface
1096 % STUPID MATLAB BUGS MEAN WE CANT allow int32 here
1097 srf= double( find_boundary(simp) );
1098 bc = ones(size(srf,1),1); % Add srf for the electrodes
1099 
1100 % Iterate over electrodes
1101 for i=1:N_elec;
1102   eleci_obj = mat_ind == (N_obj - N_elec + i);
1103   this_elec = unique( simp0( eleci_obj, : ));
1104   eleci_nodes = vtx_renum( intersect( this_elec, in )); 
1105 
1106 % This is the direct way to get electrodes. Instead we need to call the
1107 %   electrode finder function
1108 % elec(i).nodes = eleci_nodes;
1109   
1110   eleci_srf = all( ismember(srf, eleci_nodes), 2);
1111   bc( eleci_srf ) = i+1; % give this elec a surface
1112 end
1113 
1114 mat_ind( remove_simp) = [];
1115 
1116 % Test code:
1117 % fmdl.type='fwd_model'; fmdl.nodes = vtx; fmdl.elems =  simp_obj; fmdl.electrode= elec;
1118 
1119 
1120 
1121 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes_old...
1122     (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1123 
1124 % total objects:
1125 N_obj = max(mat_ind);
1126 
1127 % The electodes are the last N_elec simps
1128 e_simp_ind = mat_ind > (N_obj - N_elec);
1129 
1130 in = unique(simp(~e_simp_ind,:));
1131 out = unique(simp(e_simp_ind,:));
1132 boundary = intersect(in,out);
1133 out = setdiff(out,boundary);
1134 
1135 ext_srf_ind = ismember(srf,out);
1136 ext_srf_ind = ext_srf_ind(:,1) | ext_srf_ind(:,2) | ext_srf_ind(:,3);
1137 
1138 srf(ext_srf_ind,:) = [];
1139 bc(ext_srf_ind,:) = [];
1140 fc(ext_srf_ind,:) = [];
1141 simp = simp(~e_simp_ind,:);
1142 mat_ind = mat_ind(~e_simp_ind);
1143 
1144 % fix bc:
1145 n_unique = numel(unique(bc));
1146 missing = setdiff(1:n_unique, unique(bc));
1147 spare = setdiff(unique(bc), 1:n_unique); 
1148 
1149 for i = 1:length(missing)
1150     bc( bc==spare(i) ) = missing(i);
1151 end
1152 
1153 % fix vtx:
1154 v = 1:size(vtx,1);
1155 unused_v = setdiff(v, union(unique(simp),unique(srf))); 
1156 v(unused_v) = [];
1157 for i = 1:size(vtx,1);
1158 %     simp_ind = find(simp == i);
1159 %     srf_ind = find( srf == i);
1160     new_v_ind = find(v == i);
1161     simp( simp == i ) = new_v_ind; 
1162     srf( srf  == i ) = new_v_ind;
1163 end
1164 vtx(unused_v,:) = [];
1165 
1166 
1167 %%
1168 function [fmdl, mat_idx] = do_unit_test
1169 fmdl = [];
1170 mat_idx = [];
1171     a = [
1172    -0.8981   -0.7492   -0.2146    0.3162    0.7935    0.9615    0.6751    0.0565   -0.3635   -0.9745
1173     0.1404    0.5146    0.3504    0.5069    0.2702   -0.2339   -0.8677   -0.6997   -0.8563   -0.4668 ]';
1174 % [fmdl, mat_idx] = ng_mk_extruded_model({2,{a,0.5*a,0.2*a},1},[16,0,1],[0.01]);
1175 % load CT2
1176 
1177 % [fmdl, mat_idx] = ng_mk_extruded_model({150,flipud(trunk),1},[16,0,75],[0.01]);
1178 
1179 % [fmdl, mat_idx] = ng_mk_extruded_model({2,{trunk/100, lung_heart_dep/100, heart/100},1},[16,1,1],[0.1]);
1180 % img = mk_image( fmdl, 1);
1181 %  img.elem_data(mat_idx{2}) = 1.1;
1182 
1183 % trunk = [    -4    -2     2     4     4     2    -2    -4
1184 %               2     4     4     2    -2    -4    -4    -2]';
1185 % heart_lung = [    -2    -1    -0.8  0.8  1     2     2    -2
1186 %                    1     2     1.8  1.8  2     1    -2    -2]';
1187 % lung = [    -2    -1    -1  -1  1     2     2    -2
1188 %             1     2     0   0  2     1    -2    -2]';
1189 % heart = [    -1    -1     1     1
1190 %               0     2     2     0]';
1191 
1192 % [fmdl, mat_idx] = ng_mk_extruded_model({2,{trunk, heart_lung, heart},1},[16,1,1],[0.1]);
1193 
1194 %  figure, show_fem( fmdl );
1195  
1196 %%
1197 xx=[
1198   -88.5777  -11.4887    4.6893   49.8963  122.7033  150.3033  195.5103 249.7573 ...
1199   258.8013  279.7393  304.9623  309.2443  322.0923  337.7963  340.6503 348.2633 ...
1200   357.3043  358.7333  361.5873  364.9183  365.3943  366.3453  366.3453 365.3943 ...
1201   362.5393  351.5943  343.5053  326.8513  299.2503  288.3073  264.9923 224.0703 ...
1202   206.4633  162.6833  106.5313   92.2543   57.5153    7.0733   -8.6297 -42.4167 ...
1203   -90.9547 -105.7057 -134.2577 -178.0367 -193.2647 -222.7687 -265.5957 -278.9197 ...
1204  -313.1817 -355.5337 -363.6237 -379.3267 -397.8857 -400.7407 -401.6927 -398.8377 ...
1205  -395.0307 -384.0867 -368.3837 -363.6247 -351.7277 -334.1217 -328.4117 -314.1357 ...
1206  -291.2947 -282.7297 -267.0257 -236.5707 -221.8187 -196.5977 -159.4807 -147.5837];
1207 
1208 yy=[
1209  -385.8513 -386.8033 -386.3273 -384.8993 -368.7193 -353.9673 -323.0363 -283.5403 ...
1210  -274.9743 -254.0363 -225.4843 -217.8703 -187.4153 -140.7813 -124.6013  -86.0573 ...
1211   -38.4703  -29.4273   -9.9173   21.0137   32.4347   53.3727   83.8257   93.3437 ...
1212   114.7587  149.0237  161.8717  187.5677  222.3037  231.3447  247.5237  267.5087 ...
1213   271.3177  277.0297  281.3127  279.4097  274.6507  273.2227  276.5547  284.6447 ...
1214   295.1127  297.4927  301.7757  304.1557  302.2537  297.4947  287.5017  282.2667 ...
1215   259.9017  225.6387  213.7427  185.6677  141.4127  125.2337   88.5917   34.8187 ...
1216    17.6897  -22.2803  -73.6723  -85.0923 -117.9263 -163.6083 -176.4573 -205.9613 ...
1217  -245.9343 -256.4023 -275.4373 -304.9403 -315.4083 -332.0623 -352.0473 -355.3783];
1218 
1219 a = [xx; yy]';
1220 a = flipud(a);
1221 % th=linspace(0,2*pi,33)'; th(end)=[];
1222 % a=[sin(th)*0.3,cos(th)];
1223 
1224 
1225      fmdl = ng_mk_extruded_model({300,a,[4,25]},[16,1.11,150],[1]);
1226      show_fem(fmdl);
ng_mk_extruded_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
