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

 See also: GET_IMG_DATA, POINT_IN_TET, MK_GRID_MODEL, NDGRID

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