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

   fmdl.mdl_slice_mapper.level = Vector [1x3] of intercepts
          of the slice on the x, y, z axis. To specify a z=2 plane
          parallel to the x,y: use levels= [inf,inf,2]
   OR
   fmdl.mdl_slice_mapper.centre and .rotate
          are the centre point and rotation matrices around the point

 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).

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 %
0016 %   fmdl.mdl_slice_mapper.level = Vector [1x3] of intercepts
0017 %          of the slice on the x, y, z axis. To specify a z=2 plane
0018 %          parallel to the x,y: use levels= [inf,inf,2]
0019 %   OR
0020 %   fmdl.mdl_slice_mapper.centre and .rotate
0021 %          are the centre point and rotation matrices around the point
0022 %
0023 % maptype
0024 %    for 'elem' map is FEM element nearest the point
0025 %    for 'node' map is FEM vertex nearest the point
0026 %    for 'nodeinterp' map a npx x npy x (Nd+1) matrix such that for each point i,j
0027 %       the nearby nodes are weighted with the corresponding element in the map(i,j).
0028 
0029 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0030 % $Id: mdl_slice_mapper.m 5112 2015-06-14 13:00:41Z aadler $
0031 
0032 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0033 copt.log_level = 4;
0034 switch maptype
0035    case 'elem';
0036       copt.fstr = 'elem_ptr';
0037       map = eidors_cache(@mdl_elem_mapper,      fmdl,copt);
0038    case 'node';
0039       copt.fstr = 'node_ptr';
0040       map = eidors_cache(@mdl_node_mapper,      fmdl,copt);
0041    case 'nodeinterp';
0042       copt.fstr = 'nodeinterp';
0043       map = eidors_cache(@mdl_nodeinterp_mapper,fmdl,copt);
0044    otherwise;
0045       error('expecting maptype = elem or node');
0046 end
0047 
0048 function elem_ptr = mdl_elem_mapper(fwd_model);
0049    NODE = level_model( fwd_model );
0050    if isfield(fwd_model.mdl_slice_mapper,'model_2d') && ...
0051            fwd_model.mdl_slice_mapper.model_2d && size(NODE,1) == 3
0052        NODE(3,:) = [];
0053    end
0054    ELEM= fwd_model.elems';
0055    if size(NODE,1) ==2 %2D
0056       [x,y] = grid_the_space( fwd_model);
0057       elem_ptr= img_mapper2( NODE, ELEM, x, y);
0058    else
0059       fmdl3 = fwd_model; fmdl3.nodes = NODE'; 
0060       [x,y] = grid_the_space( fmdl3 );
0061       elem_ptr= img_mapper3( NODE, ELEM, x, y);
0062    end
0063 
0064 
0065 function ninterp_ptr = mdl_nodeinterp_mapper(fwd_model);
0066    elem_ptr = mdl_elem_mapper(fwd_model);
0067    NODE = level_model( fwd_model );
0068    fwd_model.nodes = NODE';
0069    [x,y] = grid_the_space( fwd_model);
0070 
0071    ndims = size(NODE,1);
0072    if  ndims == 2;  NODEz = []; else; NODEz= 0; end
0073    ninterp_ptr = zeros(length(x(:)),ndims+1); % reshape later
0074 
0075    for i= find( elem_ptr(:)>0 )'; % look for all x,y inside elements
0076      nodes_i = fwd_model.elems(elem_ptr(i),:);
0077      int_fcn = inv( [ones(1,ndims+1);NODE(:,nodes_i)] );
0078      ninterp_ptr(i,:) = ( int_fcn *[1;x(i);y(i);NODEz] )';
0079    end
0080    ninterp_ptr = reshape( ninterp_ptr, size(x,1), size(x,2), ndims + 1);
0081 
0082 function node_ptr = mdl_node_mapper(fwd_model);
0083    NODE = level_model( fwd_model );
0084    [x,y] = grid_the_space( fwd_model);
0085    node_ptr= node_mapper( NODE, fwd_model.elems', fwd_model.boundary, x, y);
0086 
0087 
0088 
0089 % Search through each element and find the points which
0090 % are in that element
0091 % NPTR is matrix npx x npy with a pointer to the
0092 % node closest to it.
0093 function NPTR= node_mapper( NODE, ELEM, bdy, x, y);
0094   [npy,npx] = size(x);
0095 
0096   NODEx= NODE(1,:);
0097   NODEy= NODE(2,:);
0098   if size(NODE,1) == 2
0099      NODEz2= 0;
0100      bdy= unique(bdy(:));
0101      in = inpolygon(x(:),y(:),NODE(1,bdy)',NODE(2,bdy)');
0102   else
0103      NODEz2= NODE(3,:).^2;
0104      % This is a slow way to get the elems outside the space, but I don't see another
0105      EPTR= img_mapper3(NODE, ELEM, x, y );
0106      in = EPTR>0;
0107   end
0108   NPTR=zeros(npy,npx);
0109 
0110 % This next operation can be vectorized, but we don't
0111 %  do it because that can make really big matrices
0112 
0113   for i= 1: npy
0114      for j= 1: npx
0115         dist2 = (NODEx-x(i,j)).^2 + (NODEy-y(i,j)).^2 + NODEz2;
0116         ff = find(dist2 == min(dist2));
0117         NPTR(i,j) = ff(1);
0118      end
0119   end
0120   NPTR(~in)= 0; % outside
0121 
0122 % Search through each element and find the points which
0123 % are in that element
0124 % EPTR is matrix npx x npy with a pointer to the
0125 % element which contains it.
0126 function EPTR= img_mapper2(NODE, ELEM, x, y );
0127   [npy,npx] = size(x);
0128   v_yx= [-y(:),x(:)];
0129   turn= [0 -1 1;1 0 -1;-1 1 0];
0130   EPTR=zeros(npy,npx);
0131   % for each element j, we get points on the simplex a,b,c
0132   %   area A = abc
0133   %   for each candidate point d,
0134   %      area AA = abd + acd + bcd
0135   %      d is in j if AA = A
0136   for j= 1: size(ELEM,2)
0137     % calculate area of three subtrianges to each candidate point.
0138     xy= NODE(:,ELEM(:,j))';
0139     % come up with a limited set of candidate points which
0140     % may be within the simplex
0141     endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0142              & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0143     % a is determinant of matrix [i,j,k, xy]
0144     a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0145 
0146     aa= sum(abs(ones(length(endr),1)*a'+ ...
0147                 v_yx(endr,:)*xy'*turn)');
0148     endr( abs( (abs(sum(a))-aa) ./ sum(a)) >1e-8)=[];
0149     EPTR(endr)= j;
0150   end %for j=1:ELEM
0151 
0152 % 2D mapper of points to elements. First, we assume that
0153 % The vertex geometry (NODE) has been rotated and translated
0154 % so that the imaging plane is on the z-axis. Then we iterate
0155 % through elements to find the containing each pixel
0156 function EPTR= img_mapper2a(NODE, ELEM, npx, npy );
0157   [x,y] = grid_the_space(npx, npy);
0158 
0159   EPTR=zeros(npy,npx);
0160   % for each element j, we get points on the simplex a,b,c
0161   %   area A = abc
0162   %   for each candidate point d,
0163   %      area AA = abd + acd + bcd
0164   %      d is in j if AA = A
0165   for j= 1: size(ELEM,2)
0166     xyz= NODE(:,ELEM(:,j))';
0167     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0168     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0169 
0170     % Simplex volume is det([v2-v1,v3-v1, ...])
0171     VOL= abs(det(xyz'*[-1,1,0;-1,0,1]'));
0172 
0173     % come up with a limited set of candidate points which
0174     % may be within the simplex
0175     endr=find( y(:)<=max_y & y(:)>=min_y ...
0176              & x(:)<=max_x & x(:)>=min_x );
0177 
0178     nn=  size(ELEM,1); %Simplex vertices
0179     ll=  length(endr);
0180     vol=zeros(ll,nn);
0181     for i=1:nn
0182        i1= i; i2= rem(i,nn)+1;
0183        x1= xyz(i1,1) - x(endr);
0184        y1= xyz(i1,2) - y(endr);
0185        x2= xyz(i2,1) - x(endr);
0186        y2= xyz(i2,2) - y(endr);
0187        vol(:,i)= x1.*y2 - x2.*y1;  % determinant
0188     end
0189 
0190     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0191     EPTR(endr)= j;
0192   end %for j=1:ELEM
0193 
0194 
0195 % 3D mapper of points to elements. First, we assume that
0196 % The vertex geometry (NODE) has been rotated and translated
0197 % so that the imaging plane is on the z-axis. Then we iterate
0198 % through elements to find the containing each pixel
0199 function EPTR= img_mapper3(NODE, ELEM, x, y );
0200   [npy,npx] = size(x);
0201 
0202   EPTR=zeros(npy,npx);
0203   % for each element j, we get points on the simplex a,b,c
0204   %   area A = abc
0205   %   for each candidate point d,
0206   %      area AA = abd + acd + bcd
0207   %      d is in j if AA = A
0208   idx = 1:size(ELEM,2);
0209   z = reshape(NODE(3,ELEM),size(ELEM));
0210   idx(min(z)>0 | max(z)<0) = [];
0211   for j= idx
0212     xyz= NODE(:,ELEM(:,j))';
0213 
0214     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0215     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0216 
0217     % Simplex volume is det([v2-v1,v3-v1, ...])
0218     VOL= abs(det(xyz'*[-1,1,0,0;-1,0,1,0;-1,0,0,1]'));
0219 
0220     % come up with a limited set of candidate points which
0221     % may be within the simplex
0222     endr=find( y(:)<=max_y & y(:)>=min_y ...
0223              & x(:)<=max_x & x(:)>=min_x );
0224 
0225     nn=  size(ELEM,1); %Simplex vertices
0226     ll=  length(endr);
0227     vol=zeros(ll,nn);
0228     for i=1:nn
0229        i1= i; i2= rem(i,nn)+1; i3= rem(i+1,nn)+1;
0230        x1= xyz(i1,1)-x(endr); y1= xyz(i1,2)-y(endr); z1= xyz(i1,3);
0231        x2= xyz(i2,1)-x(endr); y2= xyz(i2,2)-y(endr); z2= xyz(i2,3);
0232        x3= xyz(i3,1)-x(endr); y3= xyz(i3,2)-y(endr); z3= xyz(i3,3);
0233        vol(:,i)= x1.*y2.*z3 - x1.*y3.*z2 - x2.*y1.*z3 + ...
0234                  x3.*y1.*z2 + x2.*y3.*z1 - x3.*y2.*z1;
0235     end
0236 
0237     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0238     EPTR(endr)= j;
0239   end %for j=1:ELEM
0240 
0241 
0242 % Level model: usage
0243 %   NODE= level_model( fwd_model, level );
0244 %
0245 % Level is a 1x3 vector specifying the x,y,z axis intercepts
0246 % NODE describes the vertices in this coord space
0247 
0248 function NODE= level_model( fwd_model )
0249    vtx= fwd_model.nodes;
0250    if mdl_dim(fwd_model) ==2 % 2D case
0251        NODE= vtx';
0252        return;
0253    end
0254 
0255    if     isfield(fwd_model.mdl_slice_mapper,'level')
0256        NODE = level_model_level(vtx, fwd_model.mdl_slice_mapper.level);
0257    elseif isfield(fwd_model.mdl_slice_mapper,'centre')
0258        rotate = fwd_model.mdl_slice_mapper.rotate;
0259        centre = fwd_model.mdl_slice_mapper.centre;
0260        NODE = ctr_norm_model(vtx, rotate, centre);
0261    else   error('mdl_slice_mapper: no field level or centre provided');
0262    end
0263    
0264 function NODE = level_model_level(vtx, level)   
0265    [nn, dims] = size(vtx);
0266    % Infinities tend to cause issues -> replace with realmax
0267    % Don't need to worry about the sign of the inf
0268    level( isinf(level) | isnan(level) ) = realmax;
0269    level( level==0 ) =     1e-10; %eps;
0270 
0271    % Step 1: Choose a centre point in the plane
0272    %  Weight the point by it's inv axis coords
0273    invlev= 1./level;
0274    ctr= invlev / sum( invlev.^2 );
0275 
0276    % Step 2: Choose basis vectors in the plane
0277    %  First is the axis furthest from ctr
0278    [jnk, s_ax]= sort( - abs(level - ctr) );
0279    v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0280    v1= v1 - ctr;
0281    v1= v1 / norm(v1);
0282 
0283    % Step 3: Get off-plane vector, by cross product
0284    v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0285    v2= v2 - ctr;
0286    v2= v2 / norm(v2);
0287    v3= cross(v1,v2);
0288 
0289    % Step 4: Get orthonormal basis. Replace v2
0290    v2= cross(v1,v3);
0291 
0292    % Step 5: Get bases to point in 'positive directions'
0293    v1= v1 * (1-2*(sum(v1)<0));
0294    v2= v2 * (1-2*(sum(v2)<0));
0295    v3= v3 * (1-2*(sum(v3)<0));
0296 
0297    NODE= [v1;v2;v3] * (vtx' - ctr'*ones(1,nn) );
0298    
0299 function NODE = ctr_norm_model(vtx, rotate, centre);
0300    [nn, dims] = size(vtx);
0301    NODE = rotate * (vtx' - centre'*ones(1,nn));
0302 
0303 % Create matrices x y which grid the space of NODE
0304 function  [x,y] = grid_the_space( fmdl);
0305 
0306   xspace = []; yspace = [];
0307   try 
0308      xspace =  fmdl.mdl_slice_mapper.x_pts;
0309      yspace =  fmdl.mdl_slice_mapper.y_pts;
0310   end
0311 
0312   if isempty(xspace)
0313      npx  = fmdl.mdl_slice_mapper.npx;
0314      npy  = fmdl.mdl_slice_mapper.npy;
0315 
0316      xmin = min(fmdl.nodes(:,1));    xmax = max(fmdl.nodes(:,1));
0317      xmean= mean([xmin,xmax]); xrange= xmax-xmin;
0318 
0319      ymin = min(fmdl.nodes(:,2));    ymax = max(fmdl.nodes(:,2));
0320      ymean= mean([ymin,ymax]); yrange= ymax-ymin;
0321 
0322      range= max([xrange, yrange]);
0323      xspace = linspace( xmean - range*0.5, xmean + range*0.5, npx );
0324      yspace = linspace( ymean + range*0.5, ymean - range*0.5, npy );
0325   end
0326 %   if size(xspace,2) == 1
0327       [x,y]=meshgrid( xspace, yspace );
0328 %   else
0329 %       x= xspace;
0330 %       y= yspace;
0331 %   end
0332 %   [x,y]=meshgrid( xspace, yspace );
0333 
0334 function do_unit_test
0335 % 2D NUMBER OF POINTS
0336    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0337    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0338    fmdl.mdl_slice_mapper.npx = 5;
0339    fmdl.mdl_slice_mapper.npy = 5;
0340    eptr = mdl_slice_mapper(fmdl,'elem');
0341    unit_test_cmp('eptr01',eptr,[ 0  0 51  0  0; 0 34 26 30  0;
0342                  62 35  4 29 55; 0 36 32 31  0; 0  0 59  0  0]);
0343 
0344    nptr = mdl_slice_mapper(fmdl,'node');
0345    unit_test_cmp('nptr01',nptr,[ 0  0 28  0  0; 0 14  7 17  0;
0346                  40 13  1  9 32; 0 23 11 20  0; 0  0 36  0  0]);
0347 
0348    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0349    unit_test_cmp('nint01a',nint(2:4,2:4,1),[ 0.8284, 1, 0.8284;1,1,1; 0.8284, 1, 0.8284], 1e-3);
0350 
0351    fmdl.mdl_slice_mapper.npx = 5;
0352    fmdl.mdl_slice_mapper.npy = 3;
0353    eptr = mdl_slice_mapper(fmdl,'elem');
0354    unit_test_cmp('eptr02',eptr,[  0  0 51 0  0;62 35  4 29 55; 0 0 59 0 0]);
0355 
0356    nptr = mdl_slice_mapper(fmdl,'node');
0357    unit_test_cmp('nptr02',nptr,[ 0 0 28 0 0; 40 13 1 9 32; 0 0 36 0 0 ]);
0358 
0359 % DIRECT POINT TESTS
0360    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0361    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0362    fmdl.mdl_slice_mapper.x_pts = linspace(-1,1,5);
0363    fmdl.mdl_slice_mapper.y_pts = [0,0.5];
0364    eptr = mdl_slice_mapper(fmdl,'elem');
0365    unit_test_cmp('eptr03',eptr,[ 62 35 4 29 55; 0 34 26 30 0]);
0366 
0367    nptr = mdl_slice_mapper(fmdl,'node');
0368    unit_test_cmp('nptr03',nptr,[ 40 13 1 9 32; 0 14 7 17 0]);
0369 
0370 % 3D NPOINTS
0371    imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0372    fmdl.mdl_slice_mapper.level = [inf,inf,1];
0373    fmdl.mdl_slice_mapper.npx = 4;
0374    fmdl.mdl_slice_mapper.npy = 4;
0375    eptr = mdl_slice_mapper(fmdl,'elem');
0376    test = zeros(4); test(2:3,2:3) = [512 228;524 533];
0377    unit_test_cmp('eptr04',eptr, test);
0378    nptr = mdl_slice_mapper(fmdl,'node');
0379    test = zeros(4); test(2:3,2:3) = [116 113;118 121];
0380    unit_test_cmp('nptr04',nptr, test);
0381 
0382    fmdl.mdl_slice_mapper.level = [inf,0,inf];
0383    eptr = mdl_slice_mapper(fmdl,'elem');
0384    test = zeros(4); test(1:4,2:3) = [ 792 777; 791 776; 515 500; 239 224];
0385    unit_test_cmp('eptr05',eptr,test);
0386 
0387    nptr = mdl_slice_mapper(fmdl,'node');
0388    test = zeros(4); test(1:2,:) = [ 80, 124, 122, 64; 17, 61, 59, 1];
0389    unit_test_cmp('nptr05',nptr,test);
0390 
0391    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0392    unit_test_cmp('nint05a',nint(2:3,2:3,1),[0,1;0,1],1e-3);
0393    
0394 % Centre and Rotate
0395    fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper,'level');
0396    fmdl.mdl_slice_mapper.centre = [0,0,1];
0397    fmdl.mdl_slice_mapper.rotate = eye(3);
0398    fmdl.mdl_slice_mapper.npx = 4;
0399    fmdl.mdl_slice_mapper.npy = 4;
0400    eptr = mdl_slice_mapper(fmdl,'elem');
0401    test = zeros(4); test(2:3,2:3) = [512 503;524 533];
0402    unit_test_cmp('eptr06',eptr, test);
0403 
0404 % SLOW
0405    imdl = mk_common_model('d3cr',[16,3]); fmdl = imdl.fwd_model;
0406    fmdl.mdl_slice_mapper.level = [inf,inf,1];
0407    fmdl.mdl_slice_mapper.npx = 64;
0408    fmdl.mdl_slice_mapper.npy = 64;
0409    t = cputime;
0410    eptr = mdl_slice_mapper(fmdl,'elem');
0411    txt = sprintf('eptr10 (t=%5.3fs)',cputime - t);
0412    test = [0,122872,122872; 0,122809,122809; 0,122809,122749];
0413    unit_test_cmp(txt,eptr(10:12,10:12),test);  
0414    
0415 
0416 % CHECK ORIENTATION
0417    imdl=mk_common_model('a2c0',16);
0418    img= mk_image(imdl,1); img.elem_data(26)=1.2;
0419    subplot(231);show_fem(img);
0420    subplot(232);show_slices(img);
0421    img.fwd_model.mdl_slice_mapper.npx= 20;
0422    img.fwd_model.mdl_slice_mapper.npy= 30;
0423    img.fwd_model.mdl_slice_mapper.level= [inf,inf,0];
0424    subplot(233);show_slices(img);
0425    img.fwd_model.mdl_slice_mapper.x_pts = [linspace(-1,1,23),.5];
0426    img.fwd_model.mdl_slice_mapper.y_pts = [linspace( 1,-1,34),.5];
0427    subplot(234);show_slices(img);
0428    
0429    imdl = mk_common_model('n3r2',[16,2]);
0430    img = mk_image(imdl,1); vh= fwd_solve(img);
0431    load datacom.mat A B;
0432    img.elem_data(A) = 1.2;
0433    img.elem_data(B) = 0.8;
0434    img.calc_colours.transparency_thresh= 0.25;
0435    show_3d_slices(img);
0436  
0437    cuts = [inf, -2.5, 1.5; inf, 10, 1.5];
0438    subplot(235);  show_3d_slices(img ,[],[],[],cuts );
0439    
0440    cuts = [inf, inf, 0.5; 
0441            1e-10, 2e-10, inf;1e-10, 1e-10, inf;2e-10, 1e-10, inf;
0442            inf  , 1e-10, inf;1e-10, inf  , inf];
0443    subplot(236);  show_3d_slices(img ,[],[],[],cuts );
0444