loadResults.asv 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133
  1. % Loads the results of dose painting optimizations into the results
  2. % structure.
  3. resultsFile = 'resultsTomo.mat';
  4. % folder containing optimization results
  5. resultsFolder = 'F:\rtflynn\projects\dosePainting\dataAndAnalysis\protonsTomoContrastPhantom\tomo\optimization';
  6. inputFileBaseName = 'inputTomo';
  7. inputFileExtension = 'txt';
  8. % phantom radii
  9. Rp = [8 10 12 15];
  10. % radii for centers of boost regions
  11. Rb = [5 4 3 2];
  12. % contrasts
  13. contrasts = [-0.1 -0.2 -0.5 -1 0.1 0.2 0.5 1 2];
  14. % parameters for the dose painting phantom
  15. pProp.dx = 0.1; % pixel side length (cm)
  16. pProp.Rt = 6; % background radius (cm)
  17. pProp.Nboost = 6; % number of inserts, equi-angularly spaced
  18. pProp.rHi = 1; % highest insert radius
  19. pProp.rLo = 0.15; % lowest insert radius
  20. pProp.rhoPhant = 1; % density of the phantom (g/cc)
  21. results = []; % results structure
  22. for m=1 %1:length(Rp)
  23. for n=1:length(Rb)
  24. for q=1:length(contrasts)
  25. pProp.Rp = Rp(m); % phantom size (cm)
  26. pProp.Rb = Rb(n); % distance between phantom center and insert centers (cm)
  27. pProp.contrast = contrasts(q); % contrast between background and each insert
  28. pProp.Npix = ceil((2*pProp.Rp + 1)/pProp.dx); % number of pixels on each side of the grid
  29. results{m,n,q}.pProp = pProp;
  30. inputFileName = [resultsFolder '/' inputFileBaseName num2str((q-1) + n*length(contrasts) + m*length(contrasts)*length(Rb)) '.txt'];
  31. [optSettings,optResults] = loadOptResults(inputFileName,'last');
  32. % save the dose distribution, beamlet weights
  33. results{m,n,q}.dose = optResults.dose;
  34. results{m,n,q}.weights = optResults.weights;
  35. % save other useful information
  36. results{m,n,q}.objFunc = optResults.objFunc;
  37. results
  38. % loop through all tissues, extracting names and dvhs
  39. for k=1:length(optResults.presc.tissue)
  40. results{m,n,q}.tissue(k).name = optResults.presc.tissue(k).name;
  41. results{m,n,q}.tissue(k).dvhbins = optResults.presc.tissue(k).dvhbins;
  42. results{m,n,q}.tissue(k).dvh = optResults.presc.tissue(k).dvh;
  43. % create the differential DVH, calculate some dose statistics
  44. dvhDiff = [diff(results{m,n,q}.tissue(k).dvh) 0];
  45. dbins = results{m,n,q}.tissue(k).dvhbins;
  46. ED = sum(dvhDiff.*dbins)/sum(dvhDiff); % expectation value of dose
  47. ED2 = sum(dvhDiff.*dbins.^2)/sum(dvhDiff); % expectation of dose squared
  48. Dvar = ED2 - ED^2; % dose variance
  49. Dstd = sqrt(Dvar); % dose standard deviation
  50. results{m,n,q}.tissue(k).Davg = ED;
  51. results{m,n,q}.tissue(k).Dstd = Dstd;
  52. end
  53. end
  54. end
  55. end
  56. %
  57. %
  58. % for m=1:M
  59. % for n=1:N
  60. % for q=1:Q
  61. % siz = results{m,n,q}.optSettings.presc.siz;
  62. % Nbeamlets = numel(results{m,n,q}.optSettings.w0);
  63. % Niterations = results{m,n,q}.optSettings.Niterations;
  64. %
  65. % % load the dose results file
  66. % doseFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
  67. % doseFile = [doseFolder '/dosebatch' num2str(Niterations) '.img'];
  68. % fid = fopen(doseFile,'rb');
  69. % dose = fread(fid,'float=>float');
  70. % fclose(fid);
  71. % dose = reshape(dose,siz);
  72. %
  73. % % load the beamlet weights
  74. % weightsFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
  75. % weightsFile = [doseFolder '/weightbatch' num2str(Niterations)
  76. % '.img'];
  77. % fid = fopen(weightsFile,'rb');
  78. % weights = fread(fid,'float=>float');
  79. % fclose(fid);
  80. % weights = reshape(weights,[1 Nbeamlets]);
  81. %
  82. % % load the objective function
  83. % objFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
  84. % objFile = [doseFolder '/' results{1,1,1}.optSettings.objFuncFileName];
  85. % fid = fopen(objFile,'rb');
  86. % objFunc = fread(fid,'float=>float');
  87. % fclose(fid);
  88. % objFunc = reshape(objFunc,[1 Niterations+1]);
  89. %
  90. % % save results
  91. % results{m,n,q}.weights = weights;
  92. % results{m,n,q}.dose = dose;
  93. % results{m,n,q}.objFunc = objFunc;
  94. %
  95. % % calculate the dvh for each tissue
  96. % for k=1:length(results{m,n,q}.optSettings.presc.tissue)
  97. % mask = zeros(siz);
  98. % mask(results{m,n,q}.optSettings.presc.tissue(k).ind) = 1;
  99. % results{m,n,q}.optSettings.presc.tissue(k).dvhbins = dbins;
  100. % results{m,n,q}.optSettings.presc.tissue(k).dvh = dvh(results{m,n,q}.dose,mask,dbins);
  101. % dvhDiff = [diff(results{m,n,q}.optSettings.presc.tissue(k).dvh) 0]; % differential DVH
  102. % ED = sum(dvhDiff.*dbins)/sum(dvhDiff); % expectation value of dose
  103. % ED2 = sum(dvhDiff.*dbins.^2)/sum(dvhDiff); % expectation of dose squared
  104. % Dvar = ED2 - ED^2; % dose variance
  105. % Dstd = sqrt(Dvar); % dose standard deviation
  106. %
  107. % results{m,n,q}.optSettings.presc.tissue(k).Davg = ED;
  108. % results{m,n,q}.optSettings.presc.tissue(k).Dstd = Dstd;
  109. % end
  110. % end
  111. % end
  112. % fprintf('Finished m = %g of %g\n',m,M);
  113. % end
  114. %
  115. % results = results;
  116. %
  117. % save(resultsFile,'results');