resample_geometry.m 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  1. function resample_geometry(Geometry, Nx, Ny, Nz)
  2. Geometry_old = Geometry;
  3. str = inputdlg({'Enter X voxel dimension (cm):', 'Enter Y voxel dimension (cm):', ...
  4. 'Enter Z voxel dimension (cm):'}, 'input', [1,8], {'0.2', '0.2', '0.2'});
  5. % dimensions in mm
  6. X_dim = str2double(str{1});
  7. Y_dim = str2double(str{2});
  8. Z_dim = str2double(str{3});
  9. disp(['Desired voxel dimensions: ' num2str(X_dim) ' ' num2str(Y_dim) ' ' num2str(Z_dim) ' mm'])
  10. % round to appropriate number of voxels
  11. Nx = round(size(Geometry_old.data,1)*Geometry_old.voxel_size(1)/X_dim);
  12. Ny = round(size(Geometry_old.data,2)*Geometry_old.voxel_size(2)/Y_dim);
  13. Nz = round(size(Geometry_old.data,3)*Geometry_old.voxel_size(3)/Z_dim);
  14. % find the correct scale factor
  15. fx = Nx/size(Geometry_old.data,1);
  16. fy = Ny/size(Geometry_old.data,2);
  17. fz = Nz/size(Geometry_old.data,3);
  18. %% start converting the Geometry file
  19. Geometry.voxel_size = Geometry_old.voxel_size./[fx, fy, fz];
  20. disp(['New voxel sizes: ' num2str(Geometry.voxel_size(1)) ' ' num2str(Geometry.voxel_size(2)) ...
  21. ' ' num2str(Geometry.voxel_size(3)) ' mm'])
  22. load('WiscPlan_preferences.mat')
  23. if ~isfield(WiscPlan_preferences,'inDataPath')
  24. Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
  25. elseif WiscPlan_preferences.inDataPath ~= 0
  26. Geometry.patient_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the new patient directory');
  27. else
  28. Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
  29. end
  30. Geometry.data = imresize3(Geometry_old.data, [Nx, Ny, Nz]);
  31. Geometry.rhomw = imresize3(Geometry_old.rhomw, [Nx, Ny, Nz]);
  32. Geometry.Smw = imresize3(Geometry_old.Smw, [Nx, Ny, Nz]);
  33. Geometry.Fmw2 = imresize3(Geometry_old.Fmw2, [Nx, Ny, Nz]);
  34. Geometry.BTV = logical(round(imresize3(double(Geometry_old.BTV), [Nx, Ny, Nz])));
  35. Geometry.Ring = logical(round(imresize3(double(Geometry_old.Ring), [Nx, Ny, Nz])));
  36. Geometry.x = Geometry.start(1) +linspace(0, Nx-1, Nx)*Geometry.voxel_size(1);
  37. Geometry.y = Geometry.start(2) +linspace(0, Ny-1, Ny)*Geometry.voxel_size(2);
  38. Geometry.z = Geometry.start(3) +linspace(0, Nz-1, Nz)*Geometry.voxel_size(3);
  39. Geometry.ROIS = {};
  40. for ROI_i = 1:numel(Geometry_old.ROIS)
  41. Geometry.ROIS{ROI_i}.name = Geometry_old.ROIS{ROI_i}.name;
  42. tabula = zeros(size(Geometry_old.data));
  43. tabula(Geometry_old.ROIS{ROI_i}.ind) = 1;
  44. tabula = imresize3(tabula, [Nx, Ny, Nz]);
  45. tabula = imgaussfilt3(tabula, 0.5);
  46. Geometry.ROIS{ROI_i}.ind = find(tabula>0.3);
  47. end
  48. %% create export folder and start populating with files
  49. mkdir([Geometry.patient_dir '\matlab_files']);
  50. save([Geometry.patient_dir '\matlab_files\Geometry.mat'], 'Geometry');
  51. % resample beamlets
  52. if false
  53. % if isfile([Geometry_old.patient_dir '\batch_dose_backup.bin'])
  54. disp('calculated beamlet file found - resampling')
  55. % resample the beamlets file
  56. f = waitbar(0, 'Resampling beamlets');
  57. beamlets_in = read_ryan_beamlets([Geometry_old.patient_dir '\batch_dose_backup.bin'], 'ryan');
  58. for beam_i = 1:numel(beamlets_in)
  59. waitbar(beam_i/numel(beamlets_in),f, ['Correcting beamlet shift: ' num2str(beam_i)])
  60. tabula = zeros(beamlets_in{beam_i}.x_count, beamlets_in{beam_i}.y_count, beamlets_in{beam_i}.z_count);
  61. beamlets_out{beam_i}.num = beamlets_in{beam_i}.num-1; %reading beamlets adds +1 to number, this fixes it
  62. beamlets_out{beam_i}.x_count = Nx;
  63. beamlets_out{beam_i}.y_count = Ny;
  64. beamlets_out{beam_i}.z_count = Nz;
  65. tabula(beamlets_in{beam_i}.non_zero_indices) = beamlets_in{beam_i}.non_zero_values;
  66. tabula = imresize3(tabula, [Nx, Ny, Nz]);
  67. beamlets_out{beam_i}.non_zero_indices = find(tabula>0.005*max(tabula(:)));
  68. beamlets_out{beam_i}.non_zero_values = tabula( beamlets_out{beam_i}.non_zero_indices);
  69. end
  70. write_ryan_beamlets([Geometry.patient_dir '\batch_dose.bin'],beamlets_out);
  71. write_ryan_beamlets([Geometry.patient_dir '\batch_dose_backup.bin'],beamlets_out);
  72. close(f)
  73. else
  74. disp('no beamlet file found')
  75. end
  76. disp('done!')
  77. end