1
0

load2geometry.m 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226
  1. function [Geometry] = load2geometry
  2. %% -----===== Get CT files =====-----
  3. load('WiscPlan_preferences.mat')
  4. if ~isfield(WiscPlan_preferences,'inDataPath')
  5. Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
  6. elseif WiscPlan_preferences.inDataPath ~= 0
  7. Geometry.data_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the directory with patient data');
  8. else
  9. Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
  10. end
  11. WiscPlan_preferences.inDataPath = Geometry.data_dir;
  12. thisDir = mfilename('fullpath');
  13. idcs = strfind(thisDir,'\');
  14. prefsdir = thisDir(1:idcs(end)-1);
  15. save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences')
  16. % Geometry.data_dir = uigetdir('C:\Box Sync\Optimal Stopping\Data\CDP', 'Select folder with patient data');
  17. % Geometry.data_dir = uigetdir('\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2', 'Select folder with patient data');
  18. [IMG_in, IMG_path, filterIdx] = uigetfile([{'*.nrrd'; '*.am'}], 'Select CT image', Geometry.data_dir);
  19. % 'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\'
  20. total_num_steps = 4;
  21. hWaitbar = waitbar(0);
  22. % -- geometry explained --
  23. % geometry_file important parts:
  24. % Geometry.ROIS - actually only to get the binary mask of the target
  25. % Geometry.data - CT image in HU
  26. % Geometry.voxel_size - size of voxels in cm
  27. % Geometry.start - image coordinate origin
  28. % these are derived from above:
  29. % +Geometry.rhomw - CT image with values in relative water density
  30. % +Geometry.Smw - very similair to rhomw
  31. % +Geometry.Fmw2
  32. % +Geometry.BTV
  33. % +Geometry.ring - a ring around PTV used to pull down IMRT optimization
  34. % --------------------------
  35. switch filterIdx
  36. case 0
  37. warning('No file selected, aborting!')
  38. return
  39. case 1
  40. [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
  41. % offset the CT, because that is what WiscPlan wants
  42. Geometry.data = Geometry.data + 1000;
  43. disp('NRRD file loaded!')
  44. Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
  45. Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
  46. [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
  47. disp(['current voxel size: ' num2str(meta.spacedirections(1)) ' ' ...
  48. num2str(meta.spacedirections(2)) ' ' num2str(meta.spacedirections(3)) ' cm.'])
  49. x = input('Is this correct? (y/n)', 's')
  50. switch x
  51. case 'n'
  52. rescale_f = input('specify resampling factor')
  53. case 'y'
  54. rescale_f = 1;
  55. end
  56. Geometry.voxel_size = meta.spacedirections * rescale_f;
  57. Geometry.start=meta.spaceorigin;
  58. case 2
  59. imgIn=am2mat([IMG_path, IMG_in]);
  60. Geometry.data = permute(imgIn.data, [2,1,3]);
  61. % offset the CT, because that is what WiscPlan wants
  62. Geometry.data = Geometry.data + 1000;
  63. disp('NRRD file loaded!')
  64. Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
  65. Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
  66. [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
  67. Geometry.voxel_size = imgIn.voxel_size;
  68. Geometry.start=imgIn.start;
  69. case 3
  70. % something selected
  71. error('please do not select *all* in file selection')
  72. otherwise
  73. % error('You should never see this.')
  74. end
  75. % Geometry.ROIS
  76. %% -----===== Get mask/contour files =====-----
  77. [TRGT_in, TRGT_path, filterIdx2] = uigetfile([{'*.nrrd'; '*.am'}], 'Select file with Target', Geometry.data_dir);
  78. switch filterIdx2
  79. case 0
  80. warning('No file selected!')
  81. case 1
  82. [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
  83. disp('NRRD file loaded!')
  84. Geometry.ROIS{1}.ind = find(TRGT_img>0);
  85. Geometry.ROIS{1}.name= 'Target';
  86. case 2
  87. TRGT_in=am2mat([TRGT_path, TRGT_in]);
  88. disp('AM file loaded!')
  89. TRGT_img = permute(TRGT_in.data, [2,1,3]);
  90. Geometry.ROIS{1}.ind = find(TRGT_img>0);
  91. Geometry.ROIS{1}.name= 'Target';
  92. case 3
  93. error('please do not select *all* in file selection')
  94. % something selected
  95. otherwise
  96. % error('You should never see this.')
  97. end
  98. %% -----===== Save geometry files =====-----
  99. patient_dir = uifile('getdir', 'Save the patient data to directory');
  100. Geometry.patient_dir = patient_dir;
  101. % patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
  102. % make directories
  103. mkdir(Geometry.patient_dir);
  104. mkdir(fullfile(Geometry.patient_dir, 'beamlet_batch_files'));
  105. mkdir(fullfile(Geometry.patient_dir, 'geometry_files'));
  106. mkdir(fullfile(Geometry.patient_dir, 'matlab_files'));
  107. mkdir(fullfile(Geometry.patient_dir, 'opt_input'));
  108. mkdir(fullfile(Geometry.patient_dir, 'opt_output'));
  109. %% -----===== BTV and Ring creation =====-----
  110. waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
  111. ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
  112. [target_idx okay] = listdlg('ListString', ROI_names, ...
  113. 'SelectionMode', 'single', 'Name', 'Target Selection', ...
  114. 'PromptString', 'Please select the target ROI. ');
  115. if okay ~= 1
  116. msgbox('Plan creation aborted');
  117. delete(hWaitbar);
  118. return;
  119. end
  120. [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
  121. 'BTV margin specification', 1, {'0.6'});
  122. if isempty(BTV_margin_answer)
  123. BTV_margin = 0.6;
  124. else
  125. BTV_margin = str2double(BTV_margin_answer{1});
  126. end
  127. % target_idx and BTV_margin are set. Expand PTV to BTV
  128. PTVmask = false(size(Geometry.rhomw));
  129. % for target_idx = target_indices
  130. PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
  131. % end
  132. if BTV_margin > 0
  133. if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
  134. bwD = bwdistsc(PTVmask, Geometry.voxel_size);
  135. Geometry.BTV = bwD <= BTV_margin;
  136. elseif BTV_margin < min(Geometry.voxel_size)
  137. warning('BTV margin too small!')
  138. Geometry.BTV = zeros(size(Geometry.data));
  139. end
  140. end
  141. % Create btv
  142. Geometry.ROIS{end+1} = Geometry.ROIS{end};
  143. Geometry.ROIS{end}.name = 'BTV';
  144. Geometry.ROIS{end}.num_curves = 0;
  145. Geometry.ROIS{end}.curves = {};
  146. Geometry.ROIS{end}.ind = find(Geometry.BTV);
  147. Geometry.ROIS{end}.visible = false;
  148. % Create ring
  149. [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
  150. 'Ring margin specification', 1, {'1'});
  151. if isempty(ring_margin_answer)
  152. ring_margin = 1;
  153. else
  154. ring_margin = str2double(ring_margin_answer{1});
  155. end
  156. bwD = bwdistsc(PTVmask, Geometry.voxel_size);
  157. Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
  158. Geometry.Ring = xor(Geometry.Ring, PTVmask);
  159. Geometry.ROIS{end+1} = Geometry.ROIS{end};
  160. Geometry.ROIS{end}.name = 'Ring';
  161. Geometry.ROIS{end}.num_curves = 0;
  162. Geometry.ROIS{end}.curves = {};
  163. Geometry.ROIS{end}.ind = find(Geometry.Ring);
  164. Geometry.ROIS{end}.visible = false;
  165. %% -----===== Save the matlab files =====-----
  166. waitbar(2.33/total_num_steps, hWaitbar, 'Step 3: Save matlab geometry files');
  167. % Save matlab geometry file
  168. save(fullfile(Geometry.patient_dir, 'matlab_files', 'Geometry.mat'), 'Geometry');
  169. waitbar(2.66/total_num_steps, hWaitbar, 'Step 3: Save raw geometry files');
  170. % Write binary geometry files
  171. fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'rhomw.bin'), 'w');
  172. fwrite(fid, Geometry.rhomw, 'single');
  173. fclose(fid);
  174. fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'Smw.bin'), 'w');
  175. fwrite(fid, Geometry.Smw, 'single');
  176. fclose(fid);
  177. fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'Fmw2.bin'), 'w');
  178. fwrite(fid, Geometry.Fmw2, 'single');
  179. fclose(fid);
  180. delete(hWaitbar);
  181. waitfor(msgbox(['Plan geometry created successfully in ' '"' Geometry.patient_dir '"']));
  182. end