mdl_slice_mapper

PURPOSE ^

MDL_SLICE_MAPPER: map pixels to FEM elements or nodes

SYNOPSIS ^

function map = mdl_slice_mapper( fmdl, maptype );

DESCRIPTION ^

 MDL_SLICE_MAPPER: map pixels to FEM elements or nodes
    map = mdl_slice_mapper( fmdl, levels, 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

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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, levels, 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 %
0014 %   fmdl.mdl_slice_mapper.level = Vector [1x3] of intercepts
0015 %          of the slice on the x, y, z axis. To specify a z=2 plane
0016 %          parallel to the x,y: use levels= [inf,inf,2]
0017 % maptype
0018 %    for 'elem' map is FEM element nearest the point
0019 %    for 'node' map is FEM vertex nearest the point
0020 %    for 'nodeinterp' map a npx x npy x (Nd+1) matrix such that for each point i,j
0021 %       the nearby nodes are weighted with the corresponding element in the map(i,j).
0022 
0023 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: mdl_slice_mapper.html 2819 2011-09-07 16:43:11Z aadler $
0025 
0026 if isstr(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 switch maptype
0029   case 'elem';       map = mdl_elem_mapper(fmdl);
0030   case 'node';       map = mdl_node_mapper(fmdl);
0031   case 'nodeinterp'; map = mdl_nodeinterp_mapper(fmdl);
0032   otherwise;   error('expecting maptype = elem or node');
0033 end
0034 
0035 function elem_ptr = mdl_elem_mapper(fwd_model);
0036    level= fwd_model.mdl_slice_mapper.level;
0037 
0038    elem_ptr = eidors_obj('get-cache', fwd_model, 'elem_ptr');
0039 
0040    if ~isempty(elem_ptr)
0041       return;
0042    end
0043 
0044    NODE = level_model( fwd_model, level );
0045    ELEM= fwd_model.elems';
0046    if size(NODE,1) ==2 %2D
0047       [x,y] = grid_the_space( fwd_model);
0048       elem_ptr= img_mapper2( NODE, ELEM, x, y);
0049    else
0050       fmdl3 = fwd_model; fmdl3.nodes = NODE'; 
0051       [x,y] = grid_the_space( fmdl3 );
0052       elem_ptr= img_mapper3( NODE, ELEM, x, y);
0053    end
0054 
0055    eidors_obj('set-cache', fwd_model, 'elem_ptr', elem_ptr);
0056    eidors_msg('mdl_slice_mapper: setting cached value', 3);
0057 
0058 function ninterp_ptr = mdl_nodeinterp_mapper(fwd_model);
0059    level= fwd_model.mdl_slice_mapper.level;
0060 
0061    ninterp_ptr = eidors_obj('get-cache', fwd_model, 'ninterp_ptr');
0062    if ~isempty(ninterp_ptr); return; end
0063 
0064 
0065    elem_ptr = mdl_elem_mapper(fwd_model);
0066    NODE = level_model( fwd_model, level );
0067    fwd_model.nodes = NODE';
0068    [x,y] = grid_the_space( fwd_model);
0069 
0070    ndims = size(NODE,1);
0071    if  ndims == 2;  NODEz = []; else; NODEz= 0; end
0072    ninterp_ptr = zeros(length(x(:)),ndims+1); % reshape later
0073 
0074    for i= find( elem_ptr(:)>0 )'; % look for all x,y inside elements
0075      nodes_i = fwd_model.elems(elem_ptr(i),:);
0076      int_fcn = inv( [ones(1,ndims+1);NODE(:,nodes_i)] );
0077      ninterp_ptr(i,:) = ( int_fcn *[1;x(i);y(i);NODEz] )';
0078    end
0079    ninterp_ptr = reshape( ninterp_ptr, size(x,1), size(x,2), ndims + 1);
0080 
0081 
0082    eidors_obj('set-cache', fwd_model, 'ninterp_ptr', ninterp_ptr);
0083    eidors_msg('mdl_slice_mapper: setting cached value', 3);
0084 
0085 function node_ptr = mdl_node_mapper(fwd_model);
0086    level= fwd_model.mdl_slice_mapper.level;
0087 
0088    node_ptr = eidors_obj('get-cache', fwd_model, 'node_ptr');
0089 
0090    if ~isempty(node_ptr)
0091       return;
0092    end
0093 
0094    NODE = level_model( fwd_model, level );
0095    [x,y] = grid_the_space( fwd_model);
0096    node_ptr= node_mapper( NODE, fwd_model.elems', fwd_model.boundary, x, y);
0097 
0098    eidors_obj('set-cache', fwd_model, 'node_ptr', node_ptr);
0099    eidors_msg('mdl_slice_mapper: setting cached value', 3);
0100 
0101 
0102 % Search through each element and find the points which
0103 % are in that element
0104 % NPTR is matrix npx x npy with a pointer to the
0105 % node closest to it.
0106 function NPTR= node_mapper( NODE, ELEM, bdy, x, y);
0107   [npy,npx] = size(x);
0108 
0109   NODEx= NODE(1,:);
0110   NODEy= NODE(2,:);
0111   if size(NODE,1) == 2
0112      NODEz2= 0;
0113      bdy= unique(bdy(:));
0114      in = inpolygon(x(:),y(:),NODE(1,bdy)',NODE(2,bdy)');
0115   else
0116      NODEz2= NODE(3,:).^2;
0117      % This is a slow way to get the elems outside the space, but I don't see another
0118      EPTR= img_mapper3(NODE, ELEM, x, y );
0119      in = EPTR>0;
0120   end
0121   NPTR=zeros(npy,npx);
0122 
0123 % This next operation can be vectorized, but we don't
0124 %  do it because that can make really big matrices
0125 
0126   for i= 1: npy
0127      for j= 1: npx
0128         dist2 = (NODEx-x(i,j)).^2 + (NODEy-y(i,j)).^2 + NODEz2;
0129         ff = find(dist2 == min(dist2));
0130         NPTR(i,j) = ff(1);
0131      end
0132   end
0133   NPTR(~in)= 0; % outside
0134 
0135 % Search through each element and find the points which
0136 % are in that element
0137 % EPTR is matrix npx x npy with a pointer to the
0138 % element which contains it.
0139 function EPTR= img_mapper2(NODE, ELEM, x, y );
0140   [npy,npx] = size(x);
0141   v_yx= [-y(:) x(:)];
0142   turn= [0 -1 1;1 0 -1;-1 1 0];
0143   EPTR=zeros(npy,npx);
0144   % for each element j, we get points on the simplex a,b,c
0145   %   area A = abc
0146   %   for each candidate point d,
0147   %      area AA = abd + acd + bcd
0148   %      d is in j if AA = A
0149   for j= 1: size(ELEM,2)
0150     % calculate area of three subtrianges to each candidate point.
0151     xy= NODE(:,ELEM(:,j))';
0152     % come up with a limited set of candidate points which
0153     % may be within the simplex
0154     endr=find( y(:)<=max(xy(:,2)) & y(:)>=min(xy(:,2)) ...
0155              & x(:)<=max(xy(:,1)) & x(:)>=min(xy(:,1)) );
0156     % a is determinant of matrix [i,j,k, xy]
0157     a= xy([2;3;1],1).*xy([3;1;2],2)- xy([3;1;2],1).*xy([2;3;1],2);
0158 
0159     aa= sum(abs(ones(length(endr),1)*a'+ ...
0160                 v_yx(endr,:)*xy'*turn)');
0161     endr( abs( (abs(sum(a))-aa) ./ sum(a)) >1e-8)=[];
0162     EPTR(endr)= j;
0163   end %for j=1:ELEM
0164 
0165 % 2D mapper of points to elements. First, we assume that
0166 % The vertex geometry (NODE) has been rotated and translated
0167 % so that the imaging plane is on the z-axis. Then we iterate
0168 % through elements to find the containing each pixel
0169 function EPTR= img_mapper2a(NODE, ELEM, npx, npy );
0170   [x,y] = grid_the_space(npx, npy);
0171 
0172   EPTR=zeros(npy,npx);
0173   % for each element j, we get points on the simplex a,b,c
0174   %   area A = abc
0175   %   for each candidate point d,
0176   %      area AA = abd + acd + bcd
0177   %      d is in j if AA = A
0178   for j= 1: size(ELEM,2)
0179     xyz= NODE(:,ELEM(:,j))';
0180     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0181     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0182 
0183     % Simplex volume is det([v2-v1,v3-v1, ...])
0184     VOL= abs(det(xyz'*[-1,1,0;-1,0,1]'));
0185 
0186     % come up with a limited set of candidate points which
0187     % may be within the simplex
0188     endr=find( y(:)<=max_y & y(:)>=min_y ...
0189              & x(:)<=max_x & x(:)>=min_x );
0190 
0191     nn=  size(ELEM,1); %Simplex vertices
0192     ll=  length(endr);
0193     vol=zeros(ll,nn);
0194     for i=1:nn
0195        i1= i; i2= rem(i,nn)+1;
0196        x1= xyz(i1,1) - x(endr);
0197        y1= xyz(i1,2) - y(endr);
0198        x2= xyz(i2,1) - x(endr);
0199        y2= xyz(i2,2) - y(endr);
0200        vol(:,i)= x1.*y2 - x2.*y1;  % determinant
0201     end
0202 
0203     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0204     EPTR(endr)= j;
0205   end %for j=1:ELEM
0206 
0207 
0208 % 3D mapper of points to elements. First, we assume that
0209 % The vertex geometry (NODE) has been rotated and translated
0210 % so that the imaging plane is on the z-axis. Then we iterate
0211 % through elements to find the containing each pixel
0212 function EPTR= img_mapper3(NODE, ELEM, x, y );
0213   [npy,npx] = size(x);
0214 
0215   EPTR=zeros(npy,npx);
0216   % for each element j, we get points on the simplex a,b,c
0217   %   area A = abc
0218   %   for each candidate point d,
0219   %      area AA = abd + acd + bcd
0220   %      d is in j if AA = A
0221   for j= 1: size(ELEM,2)
0222     xyz= NODE(:,ELEM(:,j))';
0223     min_z= min(xyz(:,3)); max_z= max(xyz(:,3));
0224     if (min_z>0 || max_z<0)
0225         continue;
0226     end
0227     min_x= min(xyz(:,1)); max_x= max(xyz(:,1));
0228     min_y= min(xyz(:,2)); max_y= max(xyz(:,2));
0229 
0230     % Simplex volume is det([v2-v1,v3-v1, ...])
0231     VOL= abs(det(xyz'*[-1,1,0,0;-1,0,1,0;-1,0,0,1]'));
0232 
0233     % come up with a limited set of candidate points which
0234     % may be within the simplex
0235     endr=find( y(:)<=max_y & y(:)>=min_y ...
0236              & x(:)<=max_x & x(:)>=min_x );
0237 
0238     nn=  size(ELEM,1); %Simplex vertices
0239     ll=  length(endr);
0240     vol=zeros(ll,nn);
0241     for i=1:nn
0242        i1= i; i2= rem(i,nn)+1; i3= rem(i+1,nn)+1;
0243        x1= xyz(i1,1)-x(endr); y1= xyz(i1,2)-y(endr); z1= xyz(i1,3);
0244        x2= xyz(i2,1)-x(endr); y2= xyz(i2,2)-y(endr); z2= xyz(i2,3);
0245        x3= xyz(i3,1)-x(endr); y3= xyz(i3,2)-y(endr); z3= xyz(i3,3);
0246        vol(:,i)= x1.*y2.*z3 - x1.*y3.*z2 - x2.*y1.*z3 + ...
0247                  x3.*y1.*z2 + x2.*y3.*z1 - x3.*y2.*z1;
0248     end
0249 
0250     endr( sum(abs(vol),2) - VOL >1e-8 )=[];
0251     EPTR(endr)= j;
0252   end %for j=1:ELEM
0253 
0254 
0255 % Level model: usage
0256 %   NODE= level_model( fwd_model, level );
0257 %
0258 % Level is a 1x3 vector specifying the x,y,z axis intercepts
0259 % NODE describes the vertices in this coord space
0260 
0261 function NODE= level_model( fwd_model, level )
0262 
0263    vtx= fwd_model.nodes;
0264    [nn, dims] = size(vtx);
0265    if dims ==2 % 2D case
0266        NODE= vtx';
0267        return;
0268    end
0269 
0270    % Infinities tend to cause issues -> replace with realmax
0271    % Don't need to worry about the sign of the inf
0272    level( isinf(level) | isnan(level) ) = realmax;
0273    level( level==0 ) =     1e-10; %eps;
0274 
0275    % Step 1: Choose a centre point in the plane
0276    %  Weight the point by it's inv axis coords
0277    invlev= 1./level;
0278    ctr= invlev / sum( invlev.^2 );
0279 
0280    % Step 2: Choose basis vectors in the plane
0281    %  First is the axis furthest from ctr
0282    [jnk, s_ax]= sort( - abs(level - ctr) );
0283    v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0284    v1= v1 - ctr;
0285    v1= v1 / norm(v1);
0286 
0287    % Step 3: Get off-plane vector, by cross product
0288    v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0289    v2= v2 - ctr;
0290    v2= v2 / norm(v2);
0291    v3= cross(v1,v2);
0292 
0293    % Step 4: Get orthonormal basis. Replace v2
0294    v2= cross(v1,v3);
0295 
0296    % Step 5: Get bases to point in 'positive directions'
0297    v1= v1 * (1-2*(sum(v1)<0));
0298    v2= v2 * (1-2*(sum(v2)<0));
0299    v3= v3 * (1-2*(sum(v3)<0));
0300 
0301    NODE= [v1;v2;v3] * (vtx' - ctr'*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 
0327   [x,y]=meshgrid( xspace, yspace );
0328 
0329 function do_unit_test
0330 % 2D NUMBER OF POINTS
0331    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0332    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0333    fmdl.mdl_slice_mapper.npx = 5;
0334    fmdl.mdl_slice_mapper.npy = 5;
0335    eptr = mdl_slice_mapper(fmdl,'elem');
0336    unit_test_cmp('eptr01',eptr,[ 0  0 51  0  0; 0 34 26 30  0;
0337                  62 35  4 29 55; 0 36 32 31  0; 0  0 59  0  0]);
0338 
0339    nptr = mdl_slice_mapper(fmdl,'node');
0340    unit_test_cmp('nptr01',nptr,[ 0  0 28  0  0; 0 14  7 17  0;
0341                  40 13  1  9 32; 0 23 11 20  0; 0  0 36  0  0]);
0342 
0343    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0344    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);
0345 
0346    fmdl.mdl_slice_mapper.npx = 5;
0347    fmdl.mdl_slice_mapper.npy = 3;
0348    eptr = mdl_slice_mapper(fmdl,'elem');
0349    unit_test_cmp('eptr02',eptr,[  0  0 51 0  0;62 35  4 29 55; 0 0 59 0 0]);
0350 
0351    nptr = mdl_slice_mapper(fmdl,'node');
0352    unit_test_cmp('nptr02',nptr,[ 0 0 28 0 0; 40 13 1 9 32; 0 0 36 0 0 ]);
0353 
0354 % DIRECT POINT TESTS
0355    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0356    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0357    fmdl.mdl_slice_mapper.x_pts = linspace(-1,1,5);
0358    fmdl.mdl_slice_mapper.y_pts = [0,0.5];
0359    eptr = mdl_slice_mapper(fmdl,'elem');
0360    unit_test_cmp('eptr03',eptr,[ 62 35 4 29 55; 0 34 26 30 0]);
0361 
0362    nptr = mdl_slice_mapper(fmdl,'node');
0363    unit_test_cmp('nptr03',nptr,[ 40 13 1 9 32; 0 14 7 17 0]);
0364 
0365 % 3D NPOINTS
0366    imdl = mk_common_model('n3r2',8); fmdl = imdl.fwd_model;
0367    fmdl.mdl_slice_mapper.level = [inf,inf,1];
0368    fmdl.mdl_slice_mapper.npx = 4;
0369    fmdl.mdl_slice_mapper.npy = 4;
0370    eptr = mdl_slice_mapper(fmdl,'elem');
0371    test = zeros(4); test(2:3,2:3) = [512 228;524 533];
0372    unit_test_cmp('eptr04',eptr, test);
0373    nptr = mdl_slice_mapper(fmdl,'node');
0374    test = zeros(4); test(2:3,2:3) = [116 113;118 121];
0375    unit_test_cmp('nptr04',nptr, test);
0376 
0377    fmdl.mdl_slice_mapper.level = [inf,0,inf];
0378    eptr = mdl_slice_mapper(fmdl,'elem');
0379    test = zeros(4); test(1:4,2:3) = [ 792 777; 791 776; 515 500; 239 224];
0380    unit_test_cmp('eptr05',eptr,test);
0381 
0382    nptr = mdl_slice_mapper(fmdl,'node');
0383    test = zeros(4); test(1:2,:) = [ 80, 124, 122, 64; 17, 61, 59, 1];
0384    unit_test_cmp('nptr05',nptr,test);
0385 
0386    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0387    unit_test_cmp('nint05a',nint(2:3,2:3,1),[0,1;0,1],1e-3);
0388 
0389
Generated on Tue 09-Aug-2011 11:38:31 by m2html © 2005