calc_grid_points

CALC_GRID_POINTS - image values at points in grid
 VAL = CALC_GRID_POINTS(img, xpts, ypts, zpts) returns the matrix of image
 values at points on a rectangular grid as defined by
 NDGRID(xpts,ypts,zpts).

 [VAL, C2F] = CALC_GRID_POINTS( ... ) also returns a mapping of points to
 img.fwd_model.elems (logical [n_elems x n_points])

 See also: GET_IMG_DATA, POINT_IN_TET, MK_GRID_MODEL, NDGRID

0001 function [val, c2f] = calc_grid_points(img, xpts, ypts, zpts)
0002 %CALC_GRID_POINTS - image values at points in grid
0003 % VAL = CALC_GRID_POINTS(img, xpts, ypts, zpts) returns the matrix of image
0004 % values at points on a rectangular grid as defined by
0005 % NDGRID(xpts,ypts,zpts).
0006 %
0007 % [VAL, C2F] = CALC_GRID_POINTS( ... ) also returns a mapping of points to
0008 % img.fwd_model.elems (logical [n_elems x n_points])
0009 %
0010 % See also: GET_IMG_DATA, POINT_IN_TET, MK_GRID_MODEL, NDGRID
0011 
0012 %TODO:
0013 % - consider adding support for leveling the model first
0014 
0015 % (C) 2021 Bartek Grychtol and Andy Adler.
0016 % License: GPL version 2 or version 3
0017 % $Id: calc_grid_points.m 6903 2024-05-28 18:38:34Z bgrychtol $
0018 
0019 
0020 if nargin == 1 && ischar(img) && strcmp(img, 'UNIT_TEST')
0021     do_unit_test;
0022     return
0023 end
0024 
0025 fmdl = img.fwd_model;
0026 data = get_img_data(img); 
0027 
0028 
0029 [fmdl, use_elem] = crop_model(fmdl, xpts, ypts, zpts);
0030 [x,y,z]= ndgrid( xpts, ypts, zpts);
0031 pts = [x(:),y(:),z(:)];
0032 point2tet = point_in_tet(fmdl, pts, eps);
0033 if nargout == 2
0034    c2f = logical(sparse(size(pts,1), size(use_elem,1)));
0035    c2f(:, use_elem) = point2tet;
0036    c2f = c2f'; % the usual orientation
0037 end
0038 
0039 if isfield(img, 'elem_data')
0040     val =  point2tet * data(use_elem,:);
0041     val = bsxfun(@rdivide,val,sum(point2tet,2));
0042 elseif isfield(img, 'node_data')
0043     copt.fstr = 'calc_grid_points';
0044     copt.cache_obj = {fmdl.nodes, fmdl.elems, pts};
0045     point2node = eidors_cache(@calc_point_node_interpolation, {fmdl, pts, point2tet}, copt);
0046     val =  point2node * data;
0047     pts_out = ~any(point2tet,2);
0048     val(pts_out) = NaN;
0049 end
0050 
0051 % val is sparse sometimes in octave ... not sure why?
0052 val = reshape(full(val), [size(x) size(data,2)]);
0053 return
0054 
0055 
0056 %-------------------------------------------------------------------------%
0057 % Calculate node interpolation matrix
0058 function p2n = calc_point_node_interpolation(fmdl, pts, point2tet)
0059     [pts_idx,els_idx] = find(point2tet);
0060     
0061     %deal with nodes shared by multiple tets
0062     [pts_idx, idx] = unique(pts_idx, 'stable');
0063     point2tet = sparse(pts_idx, els_idx(idx), 1, size(point2tet,1), size(point2tet,2));
0064     
0065     elem_ptr = find(any(point2tet,1));
0066     ndims = 3;
0067 
0068     
0069     pts_map = zeros(size(pts,1),1); pts_map(pts_idx) = 1:length(pts_idx);
0070     pts_idx = repmat(pts_idx,1,ndims+1);
0071 
0072     nds_idx = zeros(size(pts_idx));
0073     int_ptr = zeros(size(pts_idx));
0074 
0075     for e= elem_ptr % loop over elements containing points
0076         nodes_i = fmdl.elems(e,:);
0077         pts_i = find(point2tet(:,e));
0078         nds_idx(pts_map(pts_i),:) = repmat(nodes_i,length(pts_i),1);
0079         
0080 %         int_fcn = inv( [ones(1,ndims+1);fmdl.nodes(nodes_i,:)'] );
0081 %         int_ptr(pts_map(pts_i),:) = ( int_fcn *[ones(numel(pts_i),1),pts(pts_i,:)]' )';
0082         
0083         % matlab claims this is better:
0084         int_ptr(pts_map(pts_i),:) = ...
0085             ([ones(1,ndims+1);fmdl.nodes(nodes_i,:)'] \ ...
0086             [ones(numel(pts_i),1),pts(pts_i,:)]')';
0087     end
0088     p2n = sparse(pts_idx,nds_idx,int_ptr, size(pts,1),size(fmdl.nodes,1));
0089    
0090 
0091 %-------------------------------------------------------------------------%
0092 % Crop parts of the model outside the grid
0093 function [mdl, use_elem] = crop_model(mdl, varargin)
0094     use_elem = true([size(mdl.elems,1),1]);
0095     for i = 1:3
0096         coord = mdl.nodes(:,i);
0097         use_elem = use_elem & any(coord(mdl.elems) >= min(varargin{i}),2);
0098         use_elem = use_elem & any(coord(mdl.elems) <= max(varargin{i}),2);
0099     end
0100     mdl.elems = mdl.elems(use_elem,:);
0101     % is it worth removing unused nodes?
0102             
0103 
0104 
0105 function do_unit_test
0106    imdl = mk_common_model('n3r2',[16,2]);
0107    img = mk_image(imdl,1);
0108    load datacom.mat A B;
0109    img.elem_data(A) = 1.2;
0110    img.elem_data(B) = 0.8;
0111    xvec = -.5:.05:.5;
0112    yvec = -1:.05:1;
0113    zvec = 1:.11:2;  
0114    rmdl = mk_grid_model([],xvec,yvec,zvec);
0115    middle_point = @(v) (v(1:end-1) + v(2:end))/2;
0116    xpts = middle_point(xvec);
0117    ypts = middle_point(yvec);
0118    zpts = middle_point(zvec);
0119    
0120    % elem data
0121    val = calc_grid_points(img, xpts, ypts, zpts);
0122    rimg = mk_image(rmdl, val(:));
0123    subplot(221)
0124    show_fem(img);
0125    hold on
0126    show_fem(rimg);
0127    hold off
0128    subplot(222)
0129    show_slices(val);
0130    
0131    % node data
0132    img.fwd_model = fix_model(img.fwd_model,struct('node2elem',1));
0133    N2E = img.fwd_model.node2elem;
0134    nimg = rmfield(img,'elem_data');
0135    nimg.node_data = N2E * img.elem_data ./ sum(N2E,2);
0136    val = calc_grid_points(nimg, xpts, ypts, zpts);
0137    
0138    subplot(223)
0139    show_fem(nimg);
0140    hold on
0141    show_fem(rimg);
0142    hold off
0143    subplot(224)
0144    show_slices(val);
0145