mdl_slice_mapper

 MDL_SLICE_MAPPER: map pixels to FEM elements or nodes
    map = mdl_slice_mapper( fmdl, maptype );

 USAGE:
 fmdl = fwd_model object
     required fields
   fmdl.mdl_slice_mapper.npx   - number of points in horizontal direction
   fmdl.mdl_slice_mapper.npy   - number of points in vertical 
    or
   fmdl.mdl_slice_mapper.x_pts - vector of points in horizontal direction
   fmdl.mdl_slice_mapper.y_pts - vector of points in vertical
     x_pts starts at the left, and y_pts starts at the top, this means
     that y_pts normally would run from max to min
    or
   fmdl.mdl_slice_mapper.npoints    - number of points across the larger
                                      dimension of the model (square)
    or
   fmdl.mdl_slice_mapper.resolution - number of points per unit
   
   fmdl.mdl_slice_mapper.level - any definition accepted by
          LEVEL_MODEL_SLICE for a single slice
   OR (deprecated)
   fmdl.mdl_slice_mapper.centre and .rotate
          are the centre point and rotation matrices around the point

      Optional fields
   fmdl.mdl_slice_mapper.model_2d - indicate a 2D model in 3D space.

 maptype
    for 'elem' map is FEM element nearest the point
    for 'node' map is FEM vertex nearest the point
    for 'nodeinterp' map a npx x npy x (Nd+1) matrix such that for each point i,j
       the nearby nodes are weighted with the corresponding element in the map(i,j).
    for 'get_points' map contains the x an y vectors of points used for
       for mapping as {x, y}. No mapping is performed.


 See also LEVEL_MODEL_SLICE

0001 function map = mdl_slice_mapper( fmdl, maptype )
0002 % MDL_SLICE_MAPPER: map pixels to FEM elements or nodes
0003 %    map = mdl_slice_mapper( fmdl, maptype );
0004 %
0005 % USAGE:
0006 % fmdl = fwd_model object
0007 %     required fields
0008 %   fmdl.mdl_slice_mapper.npx   - number of points in horizontal direction
0009 %   fmdl.mdl_slice_mapper.npy   - number of points in vertical
0010 %    or
0011 %   fmdl.mdl_slice_mapper.x_pts - vector of points in horizontal direction
0012 %   fmdl.mdl_slice_mapper.y_pts - vector of points in vertical
0013 %     x_pts starts at the left, and y_pts starts at the top, this means
0014 %     that y_pts normally would run from max to min
0015 %    or
0016 %   fmdl.mdl_slice_mapper.npoints    - number of points across the larger
0017 %                                      dimension of the model (square)
0018 %    or
0019 %   fmdl.mdl_slice_mapper.resolution - number of points per unit
0020 %
0021 %   fmdl.mdl_slice_mapper.level - any definition accepted by
0022 %          LEVEL_MODEL_SLICE for a single slice
0023 %   OR (deprecated)
0024 %   fmdl.mdl_slice_mapper.centre and .rotate
0025 %          are the centre point and rotation matrices around the point
0026 %
0027 %      Optional fields
0028 %   fmdl.mdl_slice_mapper.model_2d - indicate a 2D model in 3D space.
0029 %
0030 % maptype
0031 %    for 'elem' map is FEM element nearest the point
0032 %    for 'node' map is FEM vertex nearest the point
0033 %    for 'nodeinterp' map a npx x npy x (Nd+1) matrix such that for each point i,j
0034 %       the nearby nodes are weighted with the corresponding element in the map(i,j).
0035 %    for 'get_points' map contains the x an y vectors of points used for
0036 %       for mapping as {x, y}. No mapping is performed.
0037 %
0038 %
0039 % See also LEVEL_MODEL_SLICE
0040 
0041 % (C) 2006-2022 Andy Adler and Bartek Grychtol.
0042 % License: GPL version 2 or version 3
0043 % $Id: mdl_slice_mapper.m 6513 2022-12-30 19:20:08Z aadler $
0044 
0045 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0046 
0047 if nargin < 3, lev_no = 1; end
0048 
0049 copt.log_level = 4;
0050 params = {fmdl, maptype};
0051 copt.cache_obj = {fmdl.nodes, fmdl.elems, fmdl.mdl_slice_mapper, maptype};
0052 copt.fstr = 'mdl_slice_mapper';
0053 switch maptype
0054    case {'elem', 'node', 'nodeinterp'}
0055       map = eidors_cache(@do_mdl_slice_mapper,      params, copt);
0056    case 'get_points'
0057       map = get_points(fmdl);  
0058    otherwise
0059       error('expecting maptype = elem or node');
0060 end
0061 
0062 function map = do_mdl_slice_mapper(fmdl, maptype)
0063 
0064     [NODE, ELEM] = level_model( fmdl);
0065     if isfield(fmdl.mdl_slice_mapper,'model_2d') && ...
0066           fmdl.mdl_slice_mapper.model_2d && ...
0067           size(fmdl.nodes,2) == 3
0068        NODE(3,:) = [];
0069     end
0070     fmdl.nodes = NODE';
0071     [x, y] = grid_the_space( fmdl);
0072    
0073     switch maptype
0074        case 'elem'
0075           map = mdl_elem_mapper(NODE, ELEM, x, y);
0076        case 'node'
0077           bnd = [];
0078           try bnd = fmdl.boundary; end
0079           map = mdl_node_mapper(NODE, ELEM, bnd, x, y);
0080        case 'nodeinterp'
0081           map = mdl_nodeinterp_mapper(NODE, ELEM, x, y);
0082     end
0083 
0084 function elem_ptr = mdl_elem_mapper(NODE, ELEM, x, y)
0085    if size(NODE,1) ==2 %2D
0086       elem_ptr= img_mapper2( NODE, ELEM, x, y);
0087    else
0088       elem_ptr= img_mapper3( NODE, ELEM, x, y);
0089    end
0090 
0091 function val = use_triangulation
0092    ver = eidors_obj('interpreter_version');
0093    val = ~ver.isoctave && ver.ver >= 9.013; % pointLocation was slow/buggy before
0094    
0095    
0096 function ninterp_ptr = mdl_nodeinterp_mapper(NODE, ELEM, x, y)
0097    if use_triangulation
0098      ninterp_ptr = mdl_nodeinterp_mapper_triangulation( ...
0099         double(NODE), double(ELEM), x, y);
0100      return
0101    end
0102    elem_ptr = mdl_elem_mapper(NODE, ELEM, x, y);
0103    
0104    ndims = size(NODE,1);
0105    if  ndims == 2;  NODEz = []; else; NODEz= 0; end
0106    ninterp_ptr = zeros(length(x(:)),ndims+1); % reshape later
0107    
0108    elems = ELEM';
0109    for i= find( elem_ptr(:)>0 )' % look for all x,y inside elements
0110      nodes_i = elems(elem_ptr(i),:);
0111      ninterp_ptr(i,:) = ( [ones(1,ndims+1);NODE(:,nodes_i)] \ [1;x(i);y(i);NODEz] )';
0112    end
0113    ninterp_ptr = reshape( ninterp_ptr, size(x,1), size(x,2), ndims + 1);
0114 
0115 function ninterp_ptr = mdl_nodeinterp_mapper_triangulation(NODE, ELEM, x, y)
0116    ndims = size(NODE,1);
0117    pts = [x(:),y(:)];
0118    if size(NODE,1) == 3, pts(:,3) = 0; end
0119       
0120    s = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');
0121    TR = triangulation(ELEM', NODE');
0122    warning(s.state, 'MATLAB:triangulation:PtsNotInTriWarnId');
0123    [el, bc] =   TR.pointLocation(pts);
0124    bc(isnan(el),:) = 0;
0125    ninterp_ptr = reshape(bc,size(x,1), size(x,2), ndims + 1);
0126    
0127 function node_ptr = mdl_node_mapper(NODE, ELEM, bnd, x, y)
0128 
0129    if use_triangulation
0130        node_ptr = node_mapper_triangulation( NODE, ELEM, x, y);
0131    else
0132        ndims = size(NODE,1);
0133        if ndims == 2
0134          % old code
0135          if isempty(bnd), bnd = find_boundary(ELEM'); end
0136          node_ptr= node_mapper( NODE, ELEM, bnd, x, y);
0137        else
0138          node_ptr= node_mapper_dsearchn( NODE, ELEM, x, y);
0139        end
0140    end
0141    
0142 % in 3D this is somewhat faster in matlab. Perfomance in octave varies with size
0143 function node_ptr = node_mapper_dsearchn( NODE, ELEM, x, y)
0144    if size(NODE,1) == 2
0145       node_ptr = dsearchn(NODE', ELEM', [x(:),y(:)], 0);
0146    else
0147       pts = [x(:),y(:)];
0148       pts(:,3) = 0;
0149       [NODE, ELEM, use_nodes] = limit_3dmodel_to_slice(NODE,ELEM);
0150       node_ptr = dsearchn(NODE', ELEM', pts , 0);
0151       in = node_ptr>0;
0152       node_ptr(in) = use_nodes(node_ptr(in));
0153    end
0154    node_ptr = reshape(node_ptr, size(x));
0155    
0156 function node_ptr = node_mapper_triangulation( NODE, ELEM, x, y)
0157     s = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');
0158     TR = triangulation(ELEM', NODE');
0159     warning(s.state, 'MATLAB:triangulation:PtsNotInTriWarnId');
0160     pts = [x(:), y(:)];
0161     if size(NODE,1) == 3
0162         pts(:,3) = 0;
0163     end
0164     id = TR.pointLocation(pts);
0165     in = ~isnan(id);
0166     node_ptr = zeros(size(in));
0167     node_ptr(in) = TR.nearestNeighbor(pts(in,:));
0168     node_ptr = reshape(node_ptr, size(x));
0169     
0170 
0171    
0172 % Search through each element and find the points which
0173 % are in that element
0174 % NPTR is matrix npx x npy with a pointer to the
0175 % node closest to it.
0176 function NPTR= node_mapper( NODE, ELEM, bdy, x, y);
0177   [npy,npx] = size(x);
0178 
0179   NODEx= NODE(1,:);
0180   NODEy= NODE(2,:);
0181   if size(NODE,1) == 2
0182      NODEz2= 0;
0183      bdy= unique(bdy(:));
0184      in = inpolygon(x(:),y(:),NODE(1,bdy)',NODE(2,bdy)');
0185   else
0186      NODEz2= NODE(3,:).^2;
0187      % This is a slow way to get the elems outside the space, but I don't see another
0188      EPTR= img_mapper3(NODE, ELEM, x, y );
0189      in = EPTR>0;
0190   end
0191   NPTR=zeros(npy,npx);
0192 
0193 % This next operation can be vectorized, but we don't
0194 %  do it because that can make really big matrices
0195 
0196   for i= 1: npy
0197      for j= 1: npx
0198         dist2 = (NODEx-x(i,j)).^2 + (NODEy-y(i,j)).^2 + NODEz2;
0199         [~, ff] = min(dist2);
0200         NPTR(i,j) = ff;
0201      end
0202   end
0203   NPTR(~in)= 0; % outside
0204 
0205 % Search through each element and find the points which
0206 % are in that element
0207 % EPTR is matrix npx x npy with a pointer to the
0208 % element which contains it.
0209 function EPTR= img_mapper2(NODE, ELEM, x, y );
0210   ver = eidors_obj('interpreter_version');
0211   if ver.isoctave
0212     id = tsearch(NODE(1,:),NODE(2,:), ELEM', x(:),y(:));
0213   else
0214     if use_triangulation 
0215        s = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');
0216        TR = triangulation(ELEM', NODE');
0217        warning(s.state, 'MATLAB:triangulation:PtsNotInTriWarnId');
0218        id = TR.pointLocation([x(:), y(:)]);
0219     else
0220        EPTR = img_mapper2_old(NODE, ELEM, x, y );
0221        return
0222     end
0223   end
0224   id(isnan(id)) = 0;
0225   EPTR = reshape(id,size(x));
0226 
0227 % Search through each element and find the points which
0228 % are in that element
0229 % EPTR is matrix npx x npy with a pointer to the
0230 % element which contains it.
0231 function EPTR= img_mapper2_old(NODE, ELEM, x, y );
0232   [npy,npx] = size(x);
0233   v_yx= [-y(:),x(:)];
0234   turn= [0 -1 1;1 0 -1;-1 1 0];
0235   EPTR=zeros(npy,npx);
0236   % for each element j, we get points on the simplex a,b,c
0237   %   area A = abc
0238   %   for each candidate point d,
0239   %      area AA = abd + acd + bcd
0240   %      d is in j if AA = A
0241   for j= 1: size(ELEM,2)
0242     % calculate area of three subtrianges to each candidate point.
0243     xy= NODE(:,ELEM(:,j))';
0244     % come up with a limited set of candidate points which
0245     % may be within the simplex
0246     endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0247              & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0248     % a is determinant of matrix [i,j,k, xy]
0249     a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0250 
0251     aa= sum(abs(ones(length(endr),1)*a'+ ...
0252                 v_yx(endr,:)*xy'*turn)');
0253     endr( abs( (abs(sum(a))-aa) ./ sum(a)) >1e-8)=[];
0254     EPTR(endr)= j;
0255   end %for j=1:ELEM
0256   
0257   
0258 % 2D mapper of points to elements. First, we assume that
0259 % The vertex geometry (NODE) has been rotated and translated
0260 % so that the imaging plane is on the z-axis. Then we iterate
0261 % through elements to find the containing each pixel
0262 function EPTR= img_mapper2a(NODE, ELEM, npx, npy );
0263   [x,y] = grid_the_space(npx, npy);
0264 
0265   EPTR=zeros(npy,npx);
0266   % for each element j, we get points on the simplex a,b,c
0267   %   area A = abc
0268   %   for each candidate point d,
0269   %      area AA = abd + acd + bcd
0270   %      d is in j if AA = A
0271   for j= 1: size(ELEM,2)
0272     xyz= NODE(:,ELEM(:,j))';
0273     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0274     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0275 
0276     % Simplex volume is det([v2-v1,v3-v1, ...])
0277     VOL= abs(det(xyz'*[-1,1,0;-1,0,1]'));
0278 
0279     % come up with a limited set of candidate points which
0280     % may be within the simplex
0281     endr=find( y(:)<=max_y & y(:)>=min_y ...
0282              & x(:)<=max_x & x(:)>=min_x );
0283 
0284     nn=  size(ELEM,1); %Simplex vertices
0285     ll=  length(endr);
0286     vol=zeros(ll,nn);
0287     for i=1:nn
0288        i1= i; i2= rem(i,nn)+1;
0289        x1= xyz(i1,1) - x(endr);
0290        y1= xyz(i1,2) - y(endr);
0291        x2= xyz(i2,1) - x(endr);
0292        y2= xyz(i2,2) - y(endr);
0293        vol(:,i)= x1.*y2 - x2.*y1;  % determinant
0294     end
0295 
0296     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0297     EPTR(endr)= j;
0298   end %for j=1:ELEM
0299 
0300 
0301 % 3D mapper of points to elements. First, we assume that
0302 % The vertex geometry (NODE) has been rotated and translated
0303 % so that the imaging plane is on the z-axis. Then we iterate
0304 % through elements to find the containing each pixel
0305 function EPTR= img_mapper3(NODE, ELEM, x, y );
0306 
0307   ver = eidors_obj('interpreter_version');
0308   if ver.isoctave 
0309       img2d = mdl_3d_to_2d(NODE, ELEM);  
0310       id = tsearch(img2d.fwd_model.nodes(:,1),img2d.fwd_model.nodes(:,2), ...
0311                    img2d.fwd_model.elems, x(:),y(:));
0312     %  id = tsearchn(NODE(:,use_nodes)', map(ELEM(:,use_elem))', [x(:),y(:),zeros(numel(x),1)]);
0313       in = ~isnan(id);
0314       id(in) = img2d.elem_data(id(in));
0315       id(~in) = 0;
0316   else
0317       if ver.ver <  9.013 % pointLocation was slow/buggy before
0318           EPTR = img_mapper3_old(NODE, ELEM, x, y);
0319           return
0320       end
0321       s = warning('off', 'MATLAB:triangulation:PtsNotInTriWarnId');
0322       TR = triangulation(double(ELEM'), double(NODE'));
0323       warning(s.state, 'MATLAB:triangulation:PtsNotInTriWarnId');
0324       pts = [x(:),y(:)]; pts(:,3) = 0;
0325       id = pointLocation(TR, pts);
0326       id(isnan(id)) = 0;     
0327   end 
0328   EPTR = reshape(id,size(x));
0329 
0330 function EPTR= img_mapper3_old(NODE, ELEM, x, y );
0331   [npy,npx] = size(x);
0332 
0333   EPTR=zeros(npy,npx);
0334   % for each element j, we get points on the simplex a,b,c
0335   %   area A = abc
0336   %   for each candidate point d,
0337   %      area AA = abd + acd + bcd
0338   %      d is in j if AA = A
0339   idx = 1:size(ELEM,2);
0340   z = reshape(NODE(3,ELEM),size(ELEM));
0341   idx(min(z)>0 | max(z)<0) = [];
0342   for j= idx
0343     xyz= NODE(:,ELEM(:,j))';
0344 
0345     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0346     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0347 
0348     % Simplex volume is det([v2-v1,v3-v1, ...])
0349     VOL= abs(det(xyz'*[-1,1,0,0;-1,0,1,0;-1,0,0,1]'));
0350 
0351     % come up with a limited set of candidate points which
0352     % may be within the simplex
0353     endr=find( y(:)<=max_y & y(:)>=min_y ...
0354              & x(:)<=max_x & x(:)>=min_x );
0355 
0356     nn=  size(ELEM,1); %Simplex vertices
0357     ll=  length(endr);
0358     vol=zeros(ll,nn);
0359     xendr = x(endr); yendr = y(endr);
0360     for i=1:nn
0361        i1= i; i2= rem(i,nn)+1; i3= rem(i+1,nn)+1;
0362        x1= xyz(i1,1)-xendr; y1= xyz(i1,2)-yendr; z1= xyz(i1,3);
0363        x2= xyz(i2,1)-xendr; y2= xyz(i2,2)-yendr; z2= xyz(i2,3);
0364        x3= xyz(i3,1)-xendr; y3= xyz(i3,2)-yendr; z3= xyz(i3,3);
0365        vol(:,i)= x1.*y2.*z3 - x1.*y3.*z2 - x2.*y1.*z3 + ...
0366                  x3.*y1.*z2 + x2.*y3.*z1 - x3.*y2.*z1;
0367     end
0368 
0369     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0370     EPTR(endr)= j;
0371   end %for j=1:ELEM
0372   
0373 function [NODE, ELEM, use_nodes, use_elem] = limit_3dmodel_to_slice(NODE,ELEM)
0374     use_elem = 1:size(ELEM,2);
0375     z = reshape(NODE(3,ELEM),size(ELEM));
0376     rng = max(z(:)) - min(z(:));
0377     use_elem(min(z)>0.01*rng | max(z)<-0.01*rng) = [];
0378     use_nodes = unique(ELEM(:,use_elem));
0379     map = zeros(size(NODE,2),1); map(use_nodes) = 1:numel(use_nodes);
0380     NODE = NODE(:,use_nodes);
0381     ELEM = map(ELEM(:,use_elem));
0382   
0383 % function to be used on a leveled model (slice at z=0)
0384 function [img2d, use_nodes] = mdl_3d_to_2d(NODE, ELEM)  
0385   [NODE, ELEM, use_nodes, use_elem] = limit_3dmodel_to_slice(NODE,ELEM);
0386   fmdl.nodes = NODE';
0387   fmdl.elems = ELEM';
0388   fmdl.type = 'fwd_model';
0389   fmdl.mdl_slice_mesher.interp_elems = false;
0390   img.fwd_model = fmdl;
0391   img.elem_data = use_elem;
0392   img.type = 'image';
0393   img2d = mdl_slice_mesher(img,[inf inf 0]);
0394   
0395 % Level model: usage
0396 %   NODE= level_model( fwd_model, level );
0397 %
0398 % Level is a 1x3 vector specifying the x,y,z axis intercepts
0399 % NODE describes the vertices in this coord space
0400 
0401 function [NODE, ELEM] = level_model( fwd_model )
0402    vtx= fwd_model.nodes;
0403    ELEM = fwd_model.elems';
0404    
0405    if mdl_dim(fwd_model) ==2 || ... % 2D case
0406          isfield(fwd_model, 'mdl_slice_mapper') && ...
0407          isfield(fwd_model.mdl_slice_mapper, 'model_2d') && ...
0408          fwd_model.mdl_slice_mapper.model_2d && ...
0409          ~isfield(fwd_model.mdl_slice_mapper, 'level')
0410      
0411        NODE= vtx';
0412        return;
0413    end
0414          
0415    if     isfield(fwd_model.mdl_slice_mapper,'level')
0416        N_slices = level_model_slice(fwd_model.mdl_slice_mapper.level);
0417        if N_slices > 1
0418            eidors_msg(['Multiple slices defined on forward model. '...
0419                'Using the first slice.'], 2);
0420        end
0421        NODE = level_model_slice(vtx, fwd_model.mdl_slice_mapper.level, 1);
0422    elseif isfield(fwd_model.mdl_slice_mapper,'centre')
0423        % just in case somebody was using that interface
0424        rotate = fwd_model.mdl_slice_mapper.rotate;
0425        centre = fwd_model.mdl_slice_mapper.centre;
0426        warning('EIDORS:MDL_SLICE_MAPPER:DeprecatedInterface',...
0427             'Specifying mdl_slice_mapper.rotate and .centre is deprecated')
0428        level = struct('centre', centre,'rotation_matrix',rotate);
0429        N_slices = level_model_slice(level);
0430        if N_slices > 1
0431            warning(['Multiple slices defined on forward model. '...
0432                'Using the first slice.'])
0433        end
0434        NODE = level_model_slice(vtx, level , 1);       
0435    else   error('mdl_slice_mapper: no field level or centre provided');
0436    end
0437    
0438    NODE = NODE'; 
0439    
0440 function pts = get_points(fwd_model);
0441    NODE = level_model( fwd_model );
0442    if isfield(fwd_model.mdl_slice_mapper,'model_2d') && ...
0443            fwd_model.mdl_slice_mapper.model_2d && size(NODE,1) == 3
0444        NODE(3,:) = [];
0445    end
0446    if size(NODE,1) ==2 %2D
0447       [x,y] = grid_the_space( fwd_model, 'only_get_points');
0448    else
0449       fmdl3 = fwd_model; fmdl3.nodes = NODE'; 
0450       [x,y] = grid_the_space( fmdl3 , 'only_get_points');
0451    end
0452    pts = {x, y};
0453   
0454    
0455 % Create matrices x y which grid the space of NODE
0456 function  [x,y] = grid_the_space( fmdl, flag);
0457 
0458   if nnz(isfield(fmdl.mdl_slice_mapper, {'x_pts', 'npoints', 'npx', 'resolution'}))>1
0459       warning('Multiple specifications of points to map provided. Precedence not guaranteed')
0460   end
0461 
0462   xspace = []; yspace = [];
0463   try 
0464      xspace =  fmdl.mdl_slice_mapper.x_pts;
0465      yspace =  fmdl.mdl_slice_mapper.y_pts;
0466   end
0467   
0468   
0469 
0470   if isempty(xspace)
0471 
0472      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0473      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0474 
0475      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0476      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0477      
0478     
0479      npx = []; npy = [];
0480      if all(isfield(fmdl.mdl_slice_mapper,{'npx','npy'}))
0481          npx  = fmdl.mdl_slice_mapper.npx;
0482          npy  = fmdl.mdl_slice_mapper.npy;
0483      elseif isfield(fmdl.mdl_slice_mapper, 'npoints')
0484          maxrange = max(xrange,yrange); 
0485          xrange = maxrange; yrange = maxrange;
0486          npx  = fmdl.mdl_slice_mapper.npoints;
0487          npy  = fmdl.mdl_slice_mapper.npoints;
0488      else
0489          res = 1;
0490          try res =  fmdl.mdl_slice_mapper.resolution; end
0491          npx = ceil(xrange*res);
0492          npy = ceil(yrange*res);
0493      end
0494      
0495      xdiv = xrange/npx; ydiv = yrange/npy;
0496      xspace = linspace( xmean - xrange/2 - xdiv/2, xmean + xrange*0.5 + xdiv/2, npx + 2);
0497      yspace = linspace( ymean + yrange/2 + ydiv/2, ymean - yrange*0.5 - ydiv/2, npy + 2);
0498      xspace = xspace(2:end-1);
0499      yspace = yspace(2:end-1);
0500 
0501   end
0502   if nargin > 1 && ischar(flag) && strcmp(flag, 'only_get_points')
0503       x = xspace; y = yspace;
0504       return
0505   end
0506   [x,y]=meshgrid( xspace, yspace );
0507 
0508 function do_unit_test
0509 % 2D NUMBER OF POINTS
0510    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0511    fmdl.nodes = 1e-15 * round(1e15*fmdl.nodes);
0512    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0513    fmdl.mdl_slice_mapper.npx = 5;
0514    fmdl.mdl_slice_mapper.npy = 5;
0515    eptr = mdl_slice_mapper(fmdl,'elem');
0516    
0517    si = @(x,y) sub2ind([5,5],x, y);
0518    % points lie on nodes and edges, we allow any associated element
0519    els = {si(2,2), [25,34];
0520           si(2,4), [27,30];
0521           si(3,3), [1,2,3,4];
0522           si(4,2), [33,36];
0523           si(4,4), [28,31];
0524           };
0525    
0526    tst = eptr == [    0   37   38   39    0
0527                      46   25    5   27   42
0528                      47    8    1    6   41
0529                      48   33    7   28   40
0530                       0   45   44   43    0];
0531                       
0532    for i = 1:size(els,1)
0533     tst(els{i,1}) = ismember(eptr(els{i,1}), els{i,2});
0534    end
0535    unit_test_cmp('eptr01',tst,true(5,5));
0536 
0537    nptr = mdl_slice_mapper(fmdl,'node');
0538    test = [   0   27   28   29    0
0539              41    6    7    8   31
0540              40   13    1    9   32
0541              39   12   11   10   33
0542               0   37   36   35    0];
0543    unit_test_cmp('nptr01',nptr,test);
0544 
0545    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0546    unit_test_cmp('nint01a',nint(2:4,2:4,1),[ 0.2627, 0.6909, 0.2627; 
0547                                              0.6906, 1,      0.6906;
0548                                              0.2627, 0.6909, 0.2627], 1e-3);
0549 
0550    fmdl.mdl_slice_mapper.npx = 6;
0551    fmdl.mdl_slice_mapper.npy = 4;
0552    eptr = mdl_slice_mapper(fmdl,'elem');
0553    res = [  0   49   38   38   52    0
0554            62   23    9   10   20   55
0555            63   24   14   13   19   54
0556             0   60   44   44   57    0];
0557    unit_test_cmp('eptr02',eptr,res);
0558 
0559    nptr = mdl_slice_mapper(fmdl,'node');
0560    res = [  0   27   15   16   29    0
0561            25    6    2    3    8   18
0562            24   12    5    4   10   19
0563             0   37   22   21   35    0];
0564    unit_test_cmp('nptr02',nptr,res);
0565 
0566 % DIRECT POINT TESTS
0567    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0568    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0569    fmdl.mdl_slice_mapper.x_pts = linspace(-.95,.95,4);
0570    fmdl.mdl_slice_mapper.y_pts = [0,0.5];
0571    eptr = mdl_slice_mapper(fmdl,'elem');
0572    unit_test_cmp('eptr03',eptr,[ 47 8 6 41; 0 25 27 0]);
0573 
0574    nptr = mdl_slice_mapper(fmdl,'node');
0575    unit_test_cmp('nptr03',nptr,[ 40 13 9 32; 0 6 8 0]);
0576 
0577 % 3D NPOINTS
0578    imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0579    fmdl.mdl_slice_mapper.level = [inf,inf,1]; % co-planar with faces
0580    fmdl.mdl_slice_mapper.npx = 4;
0581    fmdl.mdl_slice_mapper.npy = 4;
0582    eptr = mdl_slice_mapper(fmdl,'elem');
0583    si = @(x,y) sub2ind([4,4],x, y);
0584    % points lie on nodes and edges, we allow any associated element
0585    els = {si(1,2), [ 42,317];
0586           si(1,3), [ 30,305];
0587           si(2,1), [ 66,341];
0588           si(2,2), [243,518];
0589           si(2,3), [231,506];
0590           si(2,4), [  6,281];
0591           si(3,1), [ 78,353];
0592           si(3,2), [252,527];
0593           si(3,3), [264,539];
0594           si(3,4), [138,413];
0595           si(4,2), [102,377];
0596           si(4,3), [114,389];};
0597    tst = eptr == zeros(4); 
0598    for i = 1:size(els,1)
0599     tst(els{i,1}) = ismember(eptr(els{i,1}), els{i,2});
0600    end
0601    unit_test_cmp('eptr04',tst, true(4));
0602    nptr = mdl_slice_mapper(fmdl,'node');
0603    test = [   0   101    99     0
0604             103   116   113    97
0605             105   118   121   111
0606               0   107   109     0];
0607    unit_test_cmp('nptr04',nptr, test);
0608 
0609    fmdl.mdl_slice_mapper.level = [inf,0.01,inf];
0610    eptr = mdl_slice_mapper(fmdl,'elem'); 
0611    test = [621   792   777   555
0612            345   516   501   279
0613            343   515   499   277
0614             69   240   225     3];
0615    unit_test_cmp('eptr05',eptr,test);
0616 
0617    nptr = mdl_slice_mapper(fmdl,'node');
0618    test = [230   250   248   222
0619            167   187   185   159
0620            104   124   122    96
0621             41    61    59    33];
0622    unit_test_cmp('nptr05',nptr,test);
0623 
0624    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0625    unit_test_cmp('nint05a',nint(2:3,2:3,4),[.1250,.1250;.8225,.8225],1e-3);
0626    
0627 % Centre and Rotate
0628    fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper,'level');
0629    fmdl.mdl_slice_mapper.centre = [0,0,0.9];
0630    fmdl.mdl_slice_mapper.rotate = eye(3);
0631    fmdl.mdl_slice_mapper.npx = 4;
0632    fmdl.mdl_slice_mapper.npy = 4;
0633    eptr = mdl_slice_mapper(fmdl,'elem');
0634    test = [  0    42    30     0
0635             66   243   229     6
0636             78   250   264   138
0637              0   102   114     0];
0638    unit_test_cmp('eptr06',eptr, test);
0639 
0640 % SLOW
0641    imdl = mk_common_model('d3cr',[16,3]); fmdl = imdl.fwd_model;
0642    fmdl.nodes = 1e-15*round(1e15*fmdl.nodes);
0643    fmdl.mdl_slice_mapper.level = [inf,inf,1];
0644    fmdl.mdl_slice_mapper.npx = 64;
0645    fmdl.mdl_slice_mapper.npy = 64;
0646    t = cputime;
0647    eptr = mdl_slice_mapper(fmdl,'elem');
0648    txt = sprintf('eptr10 (t=%5.3fs)',cputime - t);
0649 % Note that triangulation gives different
0650 % results near the edge
0651    test = [122872   122872   122748
0652            122809   122749   122689
0653            122749   122749   122689];
0654    unit_test_cmp(txt,eptr(10:12,11:13),test);  
0655    
0656 
0657 % CHECK ORIENTATION
0658    imdl=mk_common_model('a2c0',16);
0659    img= mk_image(imdl,1); img.elem_data(26)=1.2;
0660    subplot(231);show_fem(img);
0661    subplot(232);show_slices(img);
0662    img.fwd_model.mdl_slice_mapper.npx= 20;
0663    img.fwd_model.mdl_slice_mapper.npy= 30;
0664    img.fwd_model.mdl_slice_mapper.level= [inf,inf,0];
0665    subplot(233);show_slices(img);
0666    img.fwd_model.mdl_slice_mapper.x_pts = [linspace(-1,1,23),.5];
0667    img.fwd_model.mdl_slice_mapper.y_pts = [linspace( 1,-1,34),.5];
0668    subplot(234);show_slices(img);
0669    
0670    imdl = mk_common_model('n3r2',[16,2]);
0671    img = mk_image(imdl,1); vh= fwd_solve(img);
0672    load datacom.mat A B;
0673    img.elem_data(A) = 1.2;
0674    img.elem_data(B) = 0.8;
0675    img.calc_colours.transparency_thresh= 0.25;
0676    show_3d_slices(img);
0677  
0678    cuts = [inf, -2.5, 1.5; inf, 10, 1.5];
0679    subplot(235);  show_3d_slices(img ,[],[],[],cuts );
0680    
0681    cuts = [inf, inf, 0.5; 
0682            1e-10, 2e-10, inf;1e-10, 1e-10, inf;2e-10, 1e-10, inf;
0683            inf  , 1e-10, inf;1e-10, inf  , inf];
0684    subplot(236);  show_3d_slices(img ,[],[],[],cuts );
0685