img_point_mapper

IMG_POINT_MAPPER - image values at points
 V = IMG_POINT_MAPPER(IMG, P, MAPTYPE ) returns the matrix of values at
 points P (N x 3) in EIDORS image IMG.
   IMG     - EIDORS image structure
   P       - a list of points (N x 3)
   MAPTYPE - specifies the value returned for each point:
     'elem'      - value of the enclosing elements
     'node'      - value of the nearest node
     'nodeinterp'- interpoloation of node values of the enclosing element
   IMG_POINT_MAPPER uses GET_IMG_DATA to obtain image values. Data must be
   defined on elements or nodes as appropriate for the chosen MAPTYPE.
   If MAPTYPE is not specified, 'elem' or 'nodeinterp' is selected.

 NOTE: Only works in Matlab. There are numerical issues before R2022b.

 See also: GET_IMG_DATA, MDL_SLICE_MAPPER

0001 function val = img_point_mapper(img, pts, maptype )
0002 %IMG_POINT_MAPPER - image values at points
0003 % V = IMG_POINT_MAPPER(IMG, P, MAPTYPE ) returns the matrix of values at
0004 % points P (N x 3) in EIDORS image IMG.
0005 %   IMG     - EIDORS image structure
0006 %   P       - a list of points (N x 3)
0007 %   MAPTYPE - specifies the value returned for each point:
0008 %     'elem'      - value of the enclosing elements
0009 %     'node'      - value of the nearest node
0010 %     'nodeinterp'- interpoloation of node values of the enclosing element
0011 %   IMG_POINT_MAPPER uses GET_IMG_DATA to obtain image values. Data must be
0012 %   defined on elements or nodes as appropriate for the chosen MAPTYPE.
0013 %   If MAPTYPE is not specified, 'elem' or 'nodeinterp' is selected.
0014 %
0015 % NOTE: Only works in Matlab. There are numerical issues before R2022b.
0016 %
0017 % See also: GET_IMG_DATA, MDL_SLICE_MAPPER
0018 
0019 % (C) 2021 Bartek Grychtol. License: GPL version 2 or version 3
0020 % $Id: img_point_mapper.m 6404 2022-11-21 20:25:12Z bgrychtol $
0021 
0022 
0023 if nargin == 1 && ischar(img) && strcmp(img, 'UNIT_TEST')
0024    do_unit_test;
0025    return
0026 end
0027 
0028 data = get_img_data(img);
0029 fmdl = img.fwd_model;
0030 if nargin < 3
0031     if size(data,1) == size(img.fwd_model.nodes,1)
0032         maptype = 'nodeinterp';
0033     else
0034         maptype = 'elem';
0035     end
0036 end
0037 
0038 TR = eidors_cache(@triangulation,{fmdl.elems, fmdl.nodes});
0039 
0040 switch maptype
0041     case 'elem'
0042         id = pointLocation(TR, pts);
0043         val = data(id,:);
0044     case 'node'
0045         id = nearestNeighbor(TR, pts);
0046         val = data(id,:);
0047     case 'nodeinterp'
0048         [id, bc] = pointLocation(TR, pts);
0049         n_pts = size(pts,1);
0050         map = builtin('sparse', repelem((1:n_pts)',1,4), fmdl.elems(id,:), bc, n_pts, size(fmdl.nodes,1));
0051         val = map * data;
0052     otherwise
0053         error('maptype must be ''elem'', ''node'', or ''nodeinterp''.')
0054 end
0055 
0056 
0057 function do_unit_test
0058 
0059     slc = mk_grid_model([], 0:10,0:10, 0:1);
0060     img = mk_image(slc, 1:100);
0061     unit_test_cmp('3D c2f elem',img_point_mapper(img, [5.5 5.5 .5], 'elem'), 56)
0062