resamplegrid.m 2.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657
  1. function [ img_out ] = resamplegrid ( img_inDat, resample_factor )
  2. %RESAMPLEGRID resample a matrix by interpolation
  3. % Usage:
  4. % [ M_resample ] = resamplegrid ( M, resample_factor )
  5. % M_resample: resampled matrix
  6. % M: input matrix
  7. % resample_factor: floating point vector, size of resampled grid is
  8. % round(resample_factor*dimX), round(resample_factor*dimY), dimZ
  9. %
  10. % Development status: under development
  11. % Copyleft (c) Xiaohu Mo
  12. % Version: 0.1
  13. %
  14. % TODO
  15. % 3d resample
  16. % different resample factor in each dimension
  17. % interpolation method
  18. if isinteger(img_inDat)
  19. warning('Converted integer input array to float');
  20. img_inDat = double(img_inDat);
  21. end
  22. is_input_flat = false;
  23. if ~isstruct(img_inDat)
  24. is_input_flat = true;
  25. img_in.data = img_inDat;
  26. img_in.start = [0 0 0];
  27. img_in.spacing = [1 1 1];
  28. end
  29. if numel(resample_factor) == 1
  30. % resample_factor = repmat(resample_factor, size(size(img_in.data)));
  31. resample_factor = [resample_factor resample_factor 1];
  32. end
  33. xg1 = img_in.start(1) + img_in.spacing(1) * (0:size(img_in.data,1)-1);
  34. xg2 = img_in.start(2) + img_in.spacing(2) * (0:size(img_in.data,2)-1);
  35. xg3 = img_in.start(3) + img_in.spacing(3) * (0:size(img_in.data,3)-1);
  36. ext_out = round(size(img_in.data) .* resample_factor);
  37. img_out.data = zeros(ext_out);
  38. img_out.spacing = img_in.spacing .* size(img_in.data) ./ size(img_out.data);
  39. img_out.start = img_in.start - img_in.spacing/2 + img_out.spacing/2;
  40. xq1 = img_out.start(1) + img_out.spacing(1) * (0:size(img_out.data,1)-1);
  41. xq2 = img_out.start(2) + img_out.spacing(2) * (0:size(img_out.data,2)-1);
  42. xq3 = img_out.start(3) + img_out.spacing(3) * (0:size(img_out.data,3)-1);
  43. [Xg1, Xg2, Xg3] = ndgrid(xg1, xg2, xg3);
  44. [Xq1, Xq2, Xq3] = ndgrid(xq1, xq2, xq3);
  45. % note: 'linear' results in all NaN in case of extrapolation
  46. img_out.data = interpn(Xg1,Xg2,Xg3,img_in.data, Xq1, Xq2, Xq3, 'spline');
  47. % F = griddedInterpolant(Xg1, Xg2, Xg3, img_in.data, 'spline');
  48. % img_out.data = F(Xq1, Xq2, Xq3);
  49. if is_input_flat
  50. img_out = img_out.data;
  51. end