loadOptResults.m 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294
  1. function [optSettings, optResults] = loadOptResults(varargin)
  2. % Loads the weights and dose distribution corresponding to a given input
  3. % filename. The input filename has the following format:
  4. %bemal
  5. % [S,R] = loadOptResults(optimizerFolder,inputFileName) loads all of the
  6. % files from
  7. % the linlsq optimization associated with inputFileName. DVHs are
  8. % calculated for all of the associated dose files.
  9. % [S,R] = loadOptResults(optimizerFolder,inputFileName,'last') loads only the last
  10. % dose and beamlet weights files, and calculates the DVHs for each
  11. % tissue type only for those files.
  12. %
  13. % inputFileName has the following format:
  14. % Niterations
  15. % 2000
  16. % Nperbatch
  17. % 50
  18. % prescription_filename
  19. % input0/prescription.txt
  20. % initial_beam_weights_filename
  21. % input0/init_beam_weights.img
  22. % beamlet_header_file
  23. % beamletbatches0/beamlet_header.txt
  24. % dose_batch_base_name
  25. % output0/dosebatch
  26. % dose_batch_extension
  27. % img
  28. % weight_batch_base_name
  29. % output0/weightbatch
  30. % weight_batch_extension
  31. % img
  32. % obj_func_name
  33. % output0/objFunc.img
  34. %
  35. % Where all paths are relative to the directory containing the
  36. % inputFileName argument.
  37. %
  38. % RTF 1/6/07
  39. warning off;
  40. % constants
  41. Ndvhbins = 1000; % number of bins to use for the cumulative DVH calculation
  42. if length(varargin) == 1
  43. error('Too few input arguments');
  44. elseif length(varargin) == 2
  45. optimizerFolder = varargin{1};
  46. optInputFile = varargin{2};
  47. readAllOutputFiles = 1;
  48. elseif length(varargin) == 3
  49. optimizerFolder = varargin{1};
  50. optInputFile = varargin{2};
  51. if strcmp(deblank(lower(varargin{3})),'last')
  52. readAllOutputFiles = 0;
  53. else
  54. readAllOutputFiles = 1;
  55. end
  56. else
  57. error('Too many input arguments');
  58. end
  59. % Extract the folder name that contains the input file
  60. inputFileNameRev = fliplr(optInputFile); % flip the input filename around
  61. % pop off the reversed file name
  62. [fileNameRev,inputFolderRev] = strtok(inputFileNameRev,{'/','\'});
  63. inputFolder = fliplr(inputFolderRev);
  64. if isempty(inputFolder) % input folder is the current folder
  65. inputFolder = [pwd '/'];
  66. end
  67. % get the patient header
  68. headerRev = fliplr(optimizerFolder);
  69. % pop off patient name
  70. [headerNameRev,optFolderRev] = strtok(headerRev,{'/','\'});
  71. patientHeader = fliplr(headerNameRev);
  72. [optSettings,missingInfo] = loadOptSettings(optimizerFolder,optInputFile); % load the optimization settings
  73. if numel(missingInfo)
  74. fprintf('Information missing from optimization files:\n');
  75. for k=1:numel(missingInfo)
  76. fprintf('%s\n',missingInfo{k});
  77. end
  78. fprintf('\n');
  79. end
  80. % find the dose grid size
  81. if isfield(optSettings,'beamletInfo') && isfield(optSettings.beamletInfo,'beamletDim')
  82. siz = optSettings.beamletInfo.beamletDim;
  83. elseif isfield(optSettings,'prescInfo') && isfield(optSettings.prescInfo,'presc') ...
  84. && isfield(optSettings.prescInfo.presc,'siz')
  85. siz = optSettings.prescInfo.presc.siz;
  86. else
  87. siz = [];
  88. end
  89. if isfield(optSettings,'optInfo') && isfield(optSettings.optInfo,'Niterations') ...
  90. && isfield(optSettings.optInfo,'Nperbatch') && isfield(optSettings.optInfo,'outputFolder')...
  91. && isfield(optSettings.optInfo,'doseBatchBaseName') && isfield(optSettings.optInfo,'doseBatchExtension')...
  92. && isfield(optSettings.optInfo,'weightBatchBaseName') && isfield(optSettings.optInfo,'weightBatchExtension')
  93. % find which dose and weights files are present
  94. % must 'ls' the opt_output folder. One file contains the objective
  95. % function, while the others are either dose batches or weight batch
  96. % files.
  97. % xmo: MUST use mod verision of ls in *nix
  98. atmp = lsl(optSettings.optInfo.outputFolder);
  99. btmp = size(atmp);
  100. if btmp(1)>4
  101. Nbatches = ceil((btmp(1) - 3)/2);
  102. elseif btmp(1)==4
  103. Nbatches = 1;
  104. else
  105. error('No dose batches found');
  106. end
  107. % Nbatches = ceil(optSettings.optInfo.Niterations/optSettings.optInfo.Nperbatch)+1; % extra '1' is for initial guess
  108. doseFileNames = cell(1,Nbatches);
  109. doseFileExist = zeros(1,Nbatches); % flags to test for existence of files
  110. weightFileNames = cell(1,Nbatches);
  111. weightFileExist = zeros(1,Nbatches);
  112. % Dave's new implementation
  113. k=0;
  114. for t=1:btmp(1)
  115. if ~isempty(findstr(optSettings.optInfo.doseBatchBaseName,mat2str(atmp(t,:))));
  116. k=k+1;
  117. batchname = strtrim(strrep(mat2str(atmp(t,:)),'''',''));
  118. str = strrep(strtok(batchname,'.'),optSettings.optInfo.doseBatchBaseName,'');
  119. batchNames{k} = str;
  120. doseFileNames{k} = [optSettings.optInfo.outputFolder '/' ...
  121. strtrim(strrep(mat2str(atmp(t,:)),'''','')) ];
  122. weightFileNames{k} = strrep(doseFileNames{k},optSettings.optInfo.doseBatchBaseName,optSettings.optInfo.weightBatchBaseName);
  123. end
  124. end
  125. for k=1:Nbatches
  126. % ryans implementation
  127. % doseFileNames{k} = [optSettings.optInfo.outputFolder '/' ...
  128. % optSettings.optInfo.doseBatchBaseName ...
  129. % num2str((k-1)*optSettings.optInfo.Nperbatch) '.' ...
  130. % optSettings.optInfo.doseBatchExtension];
  131. % weightFileNames{k} = [optSettings.optInfo.outputFolder '/' ...
  132. % optSettings.optInfo.weightBatchBaseName ...
  133. % num2str((k-1)*optSettings.optInfo.Nperbatch) '.' ...
  134. % optSettings.optInfo.weightBatchExtension];
  135. fid = fopen(doseFileNames{k},'rb');
  136. if fid == -1 % batch file doesn't exist
  137. doseFileExist(k) = 0;
  138. else
  139. fclose(fid);
  140. doseFileExist(k) = 1; % mark that file exists
  141. end
  142. fid = fopen(weightFileNames{k},'rb');
  143. if fid == -1 % batch file doesn't exist
  144. weightFileExist(k) = 0;
  145. else
  146. fclose(fid);
  147. weightFileExist(k) = 1; % mark that file exists
  148. end
  149. end
  150. optResults = [];
  151. % assign names to optResults structure array
  152. optResults.batchNames = batchNames;
  153. if readAllOutputFiles == 1
  154. optResults.dose = {};
  155. optResults.weights = {};
  156. % read in the appropriate dose files
  157. for k=1:length(doseFileNames)
  158. if doseFileExist(k) == 1
  159. fid = fopen(doseFileNames{k},'rb');
  160. dose = fread(fid,'float=>float');
  161. fclose(fid);
  162. if numel(siz) == 3 % reshape dose if siz definedre
  163. dose = reshape(dose,siz);
  164. end
  165. optResults.dose{k} = dose;
  166. end
  167. end
  168. % read in the appropriate weights files
  169. for k=1:length(weightFileNames)
  170. if weightFileExist(k) == 1
  171. fid = fopen(weightFileNames{k},'rb');
  172. weights = fread(fid,'float=>float');
  173. fclose(fid);
  174. optResults.weights{k} = weights;
  175. end
  176. end
  177. else % read in only the last dose and weights files
  178. for k=length(doseFileNames):-1:0
  179. if doseFileExist(k) == 1 % found last dose file
  180. break;
  181. end
  182. end
  183. fid = fopen(doseFileNames{k},'rb');
  184. dose = fread(fid,'float=>float');
  185. fclose(fid);
  186. if numel(siz) == 3
  187. dose = reshape(dose,siz);
  188. end
  189. optResults.dose = dose;
  190. fid = fopen(weightFileNames{k},'rb');
  191. weights = fread(fid,'float=>float');
  192. fclose(fid);
  193. optResults.weights = weights;
  194. end
  195. % load DVH information
  196. if isfield(optSettings,'prescInfo') & isfield(optSettings.prescInfo,'presc') ...
  197. & isfield(optSettings.prescInfo.presc,'tissue') ...
  198. & isfield(optSettings.prescInfo.presc.tissue,'ind')
  199. optResults.presc = optSettings.prescInfo.presc;
  200. if readAllOutputFiles == 1
  201. % calculate a max dose vector
  202. for k=1:length(optResults.dose)
  203. dmax(k) = 1.1*max(optResults.dose{k}(:));
  204. end
  205. % calculate DVHs for each tissue, for each batch number
  206. for m=1:length(optResults.presc.tissue)
  207. tissMask = zeros(optResults.presc.siz,'int8');
  208. tissMask(optResults.presc.tissue(m).ind) = 1;
  209. for k=1:length(optResults.dose)
  210. dvhbins = [0:Ndvhbins-1]*dmax(k)/Ndvhbins;
  211. optResults.presc.tissue(m).dvhbins{k} = dvhbins;
  212. optResults.presc.tissue(m).dvh{k} ...
  213. = single(dvh(optResults.dose{k},tissMask,dvhbins))';
  214. end
  215. end
  216. else
  217. dmax = 1.1*max(optResults.dose(:));
  218. % calculate DVHs and EVHs for each tissue
  219. for m=1:length(optResults.presc.tissue)
  220. tissMask = zeros(optResults.presc.siz,'int8');
  221. tissMask(optResults.presc.tissue(m).ind) = 1;
  222. dvhbins = [0:Ndvhbins-1]*dmax/Ndvhbins;
  223. optResults.presc.tissue(m).dvhbins = dvhbins;
  224. optResults.presc.tissue(m).dvh ...
  225. = single(dvh(optResults.dose,tissMask,dvhbins))';
  226. % calculate the EDVH if presc is non-zero
  227. if sum(optResults.presc.tissue(m).dMinus) > 0
  228. % extract dose voxels
  229. doseRatio = zeros(optResults.presc.siz);
  230. doseRatio(optResults.presc.tissue(m).ind) = ...
  231. optResults.dose(optResults.presc.tissue(m).ind)./ ...
  232. optResults.presc.tissue(m).dMinus*100;
  233. % remove and Infs and NaNs
  234. doseRatio(isinf(doseRatio)) = 0;
  235. doseRatio(isnan(doseRatio)) = 0;
  236. % calculate edvh
  237. edvhbins = [0:Ndvhbins-1]*max(doseRatio(:))/Ndvhbins;
  238. optResults.presc.tissue(m).edvhbins = edvhbins;
  239. optResults.presc.tissue(m).edvh = ...
  240. single(dvh(doseRatio,tissMask,edvhbins))';
  241. end
  242. end
  243. end
  244. else
  245. fprintf('Unable to load DVH information: not enough prescription fields available.\n\n');
  246. end
  247. else
  248. fprintf('Unable to load beamlet weights and dose distributions: not enough fields available.\n\n');
  249. end
  250. % load the objective function
  251. if isfield(optSettings,'optInfo') & isfield(optSettings.optInfo,'outputFolder') ...
  252. & isfield(optSettings.optInfo,'objFuncFileName')
  253. fid = fopen([optSettings.optInfo.outputFolder '/' optSettings.optInfo.objFuncFileName],'rb');
  254. if fid ~= -1
  255. optResults.objFunc = fread(fid,'float=>float');
  256. fclose(fid);
  257. end
  258. else
  259. fprintf('Unable to load objective function file: not enough fields available.');
  260. end
  261. % sort all the loaded cell arrays based on batchNames
  262. batch_indices = cellfun(@str2num, optResults.batchNames);
  263. [delme perm_vector] = sort(batch_indices);
  264. optResults.batchNames = optResults.batchNames(perm_vector);
  265. optResults.dose = optResults.dose(perm_vector);
  266. optResults.weights = optResults.weights(perm_vector);
  267. % save file to matlab_files folder
  268. savefile = fullfile(optimizerFolder, 'matlab_files', 'optResults.mat');
  269. save(savefile,'optSettings','optResults');