ng_mk_extruded_model

PURPOSE ^

NG_MAKE_EXTRUDED_MODEL: create extruded models using netgen

SYNOPSIS ^

function [fmdl,mat_idx] = ng_mk_extruded_model(shape, elec_pos, elec_shape,extra_ng_code)

DESCRIPTION ^

 NG_MAKE_EXTRUDED_MODEL: create extruded models using netgen
 [fmdl,mat_idx] = ng_mk_extruded_models(trunk_shape, elec_pos, ...
                 elec_shape, extra_ng_code);
 INPUT:
 trunk_shape = { height,[x,y],curve_type,maxsz}
   height      -> if height = 0 calculate a 2D model
   [x,y]       -> N-by-2 CLOCKWISE list of points defining the 2D shape
                  NOTE: Use a cell array to specify additional curves for
                  internal objects
   curve_type  -> 1 - interpret as vertices (default)
                  2 - interpret as splines with de Boor points at even 
                  indices (legacy)
                  3 - interpolate points (piecewise polynomial
                  interpolation). Syntax [3, N] also specifies the number
                  of samples to create.
                  4 - interpolate points with Fourier descriptor. Syntax 
                  [4, N] also specifies the number of samples to create.
                  NOTE: If additional curves are specified, curve_type can
                  also be a cell array. Otherwise, curve_type defaults to
                  1 for internal shapes.
   maxsz       -> max size of mesh elems (default = course mesh)

 ELECTRODE POSITIONS:
  elec_pos = [n_elecs_per_plane,spacing,z_planes] where spacing is either
             0 for even spacing w.r.t angular positions (0,15,30... deg)
             or
             1 for equal distances between electrodes
             Any fractional part (e.g. 0.15) is interpreted as a starting
             position -- fraction of 2*pi for values spacing < 1 and
             fraction of total perimeter for spacing > 1.
     OR
  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)

 ELECTRODE SHAPES::
  elec_shape = [width,height, [maxsz, pem, discr]]  % Rectangular elecs
     OR
  elec_shape = [radius, [0, maxsz, pem, discr] ]  % Circular elecs
     radius      -> specify -1 for point electrodes
     maxsz (OPT) -> max size of mesh elems (default = course mesh),
                    ignored if <= 0 
     pem  (OPT)  -> 1: Point Electrode Model, 0: Complete Electrode Model (default)
     discr (OPT) -> discretize electrode normals (0 = ignore = default)
                    Adjusting this value helps Netgen problems with
                    electrodes facing each other.

 Specify either a common electrode shape or for each electrode

 CITATION_REQUEST:
 AUTHOR: B Grychtol et al.
 TITLE: Impact of model shape mismatch on reconstruction quality in
 Electrical Impedance Tomography
 JOURNAL: IEEE Trans Med Imag
 YEAR: 2012
 VOL: 31
 NUM: 9
 DOI: 10.1109/TMI.2012.2200904
 PDF: http://www.sce.carleton.ca/faculty/adler/publications/2012/
      grychtol-2012-model-shape-EIT.pdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005