vb_voximg2surimg

 MAP VOXEL IMAGES TO SURFACE IMAGE USING INVERSE DISTANCE WEIGHTING
 INTERPOLATION OR NEAREST NEIGHBOR SELECTION

 [Input]
 XYZvox : voxel image coordinates with unit mm [Nvox,3]
 Tvox   : voxel value [Nvox,1]
 XYZsurf : coordinates of surface mesh points with unit mm [Nv,3] 
 r2thres : distance threshold with unit mm^2 
 
 [Optional]
 mode (=1)  : flag to switch mapping mode
              1 = inverse distance weighting interpolation
              2 = nearest neighbor selection

 [Output]
 Tsurf  : mapped surface image


 2016/07/06 O.Yamashita
 * add eps to distance calculation to prevent inverse of distance from
 being infinity
 2016/07/01 O.Yamashita

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function Tsurf = vb_voximg2surimg(XYZvox,Tvox,XYZsurf, r2thres, mode)
0002 % MAP VOXEL IMAGES TO SURFACE IMAGE USING INVERSE DISTANCE WEIGHTING
0003 % INTERPOLATION OR NEAREST NEIGHBOR SELECTION
0004 %
0005 % [Input]
0006 % XYZvox : voxel image coordinates with unit mm [Nvox,3]
0007 % Tvox   : voxel value [Nvox,1]
0008 % XYZsurf : coordinates of surface mesh points with unit mm [Nv,3]
0009 % r2thres : distance threshold with unit mm^2
0010 %
0011 % [Optional]
0012 % mode (=1)  : flag to switch mapping mode
0013 %              1 = inverse distance weighting interpolation
0014 %              2 = nearest neighbor selection
0015 %
0016 % [Output]
0017 % Tsurf  : mapped surface image
0018 %
0019 %
0020 % 2016/07/06 O.Yamashita
0021 % * add eps to distance calculation to prevent inverse of distance from
0022 % being infinity
0023 % 2016/07/01 O.Yamashita
0024 %
0025 % Copyright (C) 2011, ATR All Rights Reserved.
0026 % License : New BSD License(see VBMEG_LICENSE.txt)
0027 
0028 if nargin < 5
0029     mode = 1;
0030 end
0031 
0032 Nv = size(XYZsurf,1);
0033 Tsurf = zeros(Nv,1);
0034 
0035 switch mode
0036     case 1,
0037         fprintf('  inverse distance weighting interpolation within sphere of radius %2.2f mm ...\n', sqrt(r2thres));
0038 
0039         for vv = 1: size(XYZsurf,1),
0040             XYZref = XYZsurf(vv,:);
0041 
0042             % square distance and inverse square distance
0043             r2 = (XYZvox(:,1)-XYZref(1)).^2+ (XYZvox(:,2)-XYZref(2)).^2+ (XYZvox(:,3)-XYZref(3)).^2 +eps;
0044             ir2 = 1./r2;
0045 
0046             % inverse distance weighting
0047             ixvox = find(r2 < r2thres);
0048             if ~isempty(ixvox)  % voxels found within sphere of radius r2thres
0049                 Tsurf(vv) = Tvox(ixvox)' * ir2(ixvox)/sum(ir2(ixvox));
0050             end
0051 
0052             if mod(vv,1000)==0, fprintf('%d..', vv); end
0053 
0054         end
0055         
0056     case 2
0057          fprintf('  nearest neighbor selection within sphere of radius %2.2f mm ...\n', sqrt(r2thres));
0058          
0059          for vv = 1: size(XYZsurf,1),
0060              XYZref = XYZsurf(vv,:);
0061              
0062              % square distance and inverse square distance
0063              r2 = (XYZvox(:,1)-XYZref(1)).^2+ (XYZvox(:,2)-XYZref(2)).^2+ (XYZvox(:,3)-XYZref(3)).^2;
0064             
0065              % nearest neighbor
0066              [r2min,jx] = min(r2);
0067              if r2min < r2thres
0068                  Tsurf(vv) = Tvox(jx);
0069              end
0070              if mod(vv,1000)==0, fprintf('%d..', vv); end
0071          end
0072    
0073     otherwise 
0074         error('mode must be 1 or 2 !');
0075             
0076 end