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.

 NOTES:
  - There are numerical issues before Matlab R2022b.
  - May be slow in Octave for large models.

 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 % NOTES:
0016 %  - There are numerical issues before Matlab R2022b.
0017 %  - May be slow in Octave for large models.
0018 %
0019 % See also: GET_IMG_DATA, MDL_SLICE_MAPPER
0020 
0021 % (C) 2021-2024 Bartek Grychtol. License: GPL version 2 or version 3
0022 % $Id: img_point_mapper.m 6712 2024-03-29 22:22:27Z bgrychtol $
0023 
0024 
0025 if nargin == 1 && ischar(img) && strcmp(img, 'UNIT_TEST')
0026    do_unit_test;
0027    return
0028 end
0029 
0030 data = get_img_data(img);
0031 fmdl = img.fwd_model;
0032 if nargin < 3
0033     if size(data,1) == size(img.fwd_model.nodes,1)
0034         maptype = 'nodeinterp';
0035     else
0036         maptype = 'elem';
0037     end
0038 end
0039 
0040 if ~exist('OCTAVE_VERSION')
0041     val = img_point_mapper_matlab(fmdl, pts, maptype, data);
0042 else
0043     val = img_point_mapper_octave(fmdl, pts, maptype, data);
0044 end
0045 
0046 
0047 function val = img_point_mapper_matlab(fmdl, pts, maptype, data)
0048 
0049 TR = eidors_cache(@triangulation,{fmdl.elems, fmdl.nodes});
0050 
0051 switch maptype
0052     case 'elem'
0053         id = pointLocation(TR, pts);
0054         val = NaN(size(id));
0055         valid = ~isnan(id);
0056         val(valid,:) = data(id(valid),:);
0057     case 'node'
0058         id = nearestNeighbor(TR, pts);
0059         val = data(id,:);
0060     case 'nodeinterp'
0061         [id, bc] = pointLocation(TR, pts);
0062         n_pts = size(pts,1);
0063         map = builtin('sparse', repelem((1:n_pts)',1,4), fmdl.elems(id,:), bc, n_pts, size(fmdl.nodes,1));
0064         val = map * data;
0065     otherwise
0066         error('maptype must be ''elem'', ''node'', or ''nodeinterp''.')
0067 end
0068 
0069 function val = img_point_mapper_octave(fmdl, pts, maptype, data)
0070 switch maptype
0071     case 'elem'
0072         id = tsearchn(fmdl.nodes, fmdl.elems, pts);
0073         val = NaN(size(id));
0074         valid = ~isnan(id);
0075         val(valid,:) = data(id(valid),:);
0076     case 'node'
0077         id = dsearchn(fmdl.nodes, fmdl.elems, pts);
0078         val = data(id,:);
0079     case 'nodeinterp'
0080         [id, bc] = tsearchn(fmdl.nodes, fmdl.elems, pts);
0081         n_pts = size(pts,1);
0082         map = builtin('sparse', repelem((1:n_pts)',1,4), fmdl.elems(id,:), bc, n_pts, size(fmdl.nodes,1));
0083         val = map * data;
0084     otherwise
0085         error('maptype must be ''elem'', ''node'', or ''nodeinterp''.')
0086 end
0087 
0088 
0089 
0090 function do_unit_test
0091 
0092     slc = mk_grid_model([], 0:10,0:10, 0:1);
0093     img = mk_image(slc, 1:100);
0094     unit_test_cmp('3D c2f elem',img_point_mapper(img, [5.5 5.5 .5], 'elem'), 56);
0095     img = rmfield(img, 'elem_data');
0096     img.node_data = 1:242;
0097     unit_test_cmp('3D c2f node',img_point_mapper(img, [5.6 5.6 .6], 'node'), 194);
0098     unit_test_cmp('3D c2f nodeinterp',img_point_mapper(img, [5.6 5.6 .6], 'nodeinterp'), 140.8, 1e-9);
0099 
0100