dicomrt2geometry.m 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155
  1. function Geometry = dicomrt2geometry(dicomrt_dir)
  2. %DICOMRT2GEOMETRY Import DicomRT dataset to Ryan's "Geometry" struct
  3. % Usage:
  4. % [ Geometry ] = dicomrt2geometry ( [input_dir] )
  5. % Input:
  6. % input_dir = directory contains DicomRT data
  7. % [] = prompt to select input_dir
  8. % Output:
  9. % Geometry = Ryan's Geometry struct used with RDX TPS
  10. %
  11. % This function imports DicomRT data set and create RDX TPS compatible Geometry
  12. % struct.
  13. %
  14. % TODO: Switch from filename based Dicom identification to StudyID based.
  15. %
  16. % See also readPinnGeometry.m
  17. %
  18. % Author: Xiaohu Mo
  19. % edits: Peter Ferjancic 2019/01
  20. Geometry = [];
  21. if nargin < 1
  22. % dicomrt_dir = uifile('getdir', 'Select the directory contains the DICOM-RT images');
  23. % !!! Grozomah
  24. load('WiscPlan_preferences.mat')
  25. if isstring(WiscPlan_preferences.inDataPath)
  26. else
  27. WiscPlan_preferences.inDataPath = 'C://';
  28. end
  29. dicomrt_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the directory contains the DICOM-RT images');
  30. WiscPlan_preferences.inDataPath = dicomrt_dir;
  31. thisDir = mfilename('fullpath');
  32. idcs = strfind(thisDir,'\');
  33. prefsdir = thisDir(1:idcs(end-1)-1);
  34. save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences')
  35. end
  36. if exist(dicomrt_dir, 'dir') ~= 7
  37. msgbox('Cannot find input directory');
  38. return;
  39. end
  40. %% -----===== Start Importing =====-----
  41. % set dicom dictionary before any dicominfo or dicomread
  42. dicomdict('factory');
  43. % change path due to the limitation of DicomRT toolbox
  44. % will cd back later
  45. previousPath = pwd;
  46. cd(dicomrt_dir);
  47. % get the list of filenames
  48. dirlist = dir(dicomrt_dir);
  49. filenames = cell(0);
  50. % remove directories
  51. for i = 1:numel(dirlist)
  52. [tilde, tilde, ext] = fileparts(dirlist(i).name);
  53. % if ~dirlist(i).isdir && strcmpi(ext, '.dcm')
  54. if ~dirlist(i).isdir && ~strcmpi(ext, '.txt') % <--- Grozomah: Fix this back to default!
  55. filenames{end+1} = dirlist(i).name;
  56. end
  57. end
  58. % get all dicominfo to a cell array the same size as filelist
  59. for i = 1:numel(filenames)
  60. dcminfo{i} = dicominfo(filenames{i});
  61. end
  62. % TODO integrity check
  63. % UID should match
  64. % should be only one RTSTRUCT file present
  65. %% -----===== Import Dicom CT data set =====-----
  66. % Create DicomRT ct list file
  67. fid = fopen('ctlist.txt', 'w');
  68. for i = 1:numel(dcminfo)
  69. if strcmpi(dcminfo{i}.Modality, 'CT')
  70. fprintf(fid, '%s\n', filenames{i});
  71. end
  72. end
  73. fclose(fid);
  74. [cellCt, ct_xmesh, ct_ymesh, ct_zmesh] = dicomrt_loadct('ctlist.txt');
  75. %% -----===== Set Downsampling flag =====-----
  76. button = questdlg('Do you want to downsample the CT dataset?', 'Downsampling', 'Yes', 'No', 'Yes');
  77. if strcmpi(button, 'yes')
  78. downsample_flag = true;
  79. else
  80. downsample_flag = false;
  81. end
  82. %% -----===== Import DicomRT structure =====-----
  83. % File DicomRT structure file
  84. for i = 1:numel(dcminfo)
  85. if strcmpi(dcminfo{i}.Modality, 'RTSTRUCT')
  86. rsfilename = fullfile(dicomrt_dir, filenames{i});
  87. end
  88. end
  89. % prompt if can not find RTSTRUCT file in this directory
  90. if isempty(rsfilename)
  91. [temp_FileName temp_PathName] = uifile('get', '*.dcm', 'Please select the DicomRT-Struct file');
  92. rsfilename = fullfile(temp_PathName, temp_FileName);
  93. end
  94. % Import the struct
  95. cellVoi = dicomrt_loadvoi(rsfilename);
  96. % Downsampling
  97. if downsample_flag == true
  98. if numel(ct_xmesh) < 128
  99. downsample_factor = 0.5;
  100. dx = ct_xmesh(2)-ct_xmesh(1);
  101. dy = ct_ymesh(2)-ct_ymesh(1);
  102. ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);
  103. ct_ymesh = linspace(ct_ymesh(1)+dy/2, ct_ymesh(end)-dy/2, numel(ct_ymesh)*downsample_factor);
  104. else
  105. downsample_factor = 0.125;
  106. downsample_factor = 1/8;
  107. dx = ct_xmesh(2)-ct_xmesh(1);
  108. dy = ct_ymesh(2)-ct_ymesh(1);
  109. ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);
  110. ct_ymesh = linspace(ct_ymesh(1)+dy/2, ct_ymesh(end)-dy/2, numel(ct_ymesh)*downsample_factor);
  111. end
  112. end
  113. ROIS = dicomrt_rasterizeVoi(cellVoi, cellCt, ct_xmesh, ct_ymesh, ct_zmesh);
  114. %% -----===== Done Importing =====-----
  115. % change the path back
  116. cd(previousPath);
  117. %% -----===== Create Geometry struct =====-----
  118. Geometry.ROIS = ROIS;
  119. Geometry.data = cellCt{2};
  120. % !!! Grozomah - added abs() to prevent errors in later functions
  121. Geometry.voxel_size = abs([ct_xmesh(2)-ct_xmesh(1) ...
  122. ct_ymesh(2)-ct_ymesh(1) ...
  123. ct_zmesh(2)-ct_zmesh(1)]);
  124. Geometry.start = [ct_xmesh(1) ct_ymesh(1) ct_zmesh(1)];
  125. if downsample_flag == true
  126. Geometry.data = resamplegrid(Geometry.data, [downsample_factor downsample_factor 1]);
  127. % Geometry.data = resamplegrid(Geometry.data, 0.5);
  128. % note that ct_xmesh, ct_ymesh, ct_zmesh is already downsampled
  129. Geometry.start(1:2) = Geometry.start(1:2) + Geometry.voxel_size(1:2) / 4;
  130. end
  131. Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
  132. Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
  133. [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);