powerDVH.m 8.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154
  1. function powerDVH(lin245cm,root245cm,quad245cm,power025,power4)
  2. addpath('C:\Documents and Settings\Steve\My Documents\Research\matlab')
  3. addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqUtil')
  4. addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input')
  5. Ndvhbins = 1000; %number of dvh bins
  6. prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower1\PTV70dosePlus.bin'); %open prescription in Matlab
  7. prescTumorRoot = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower2\PTV70dosePlus.bin');
  8. prescTumorQuad = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower3\PTV70dosePlus.bin');
  9. prescTumorPower025 = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower4\PTV70dosePlus.bin');
  10. prescTumorPower4 = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower5\PTV70dosePlus.bin');
  11. %dim = 1; %2D
  12. dim = 30; %3D
  13. fid = fopen(lin245cm,'rb'); %open plans in Matlab
  14. planLin245cm = reshape(fread(fid,'float'),72,72,dim);
  15. fclose(fid);
  16. fid = fopen(root245cm,'rb'); %open plans in Matlab
  17. planRoot245cm = reshape(fread(fid,'float'),72,72,dim);
  18. fclose(fid);
  19. fid = fopen(quad245cm,'rb'); %open plans in Matlab
  20. planQuad245cm = reshape(fread(fid,'float'),72,72,dim);
  21. fclose(fid);
  22. fid = fopen(power025,'rb'); %open plans in Matlab
  23. planPower025 = reshape(fread(fid,'float'),72,72,dim);
  24. fclose(fid);
  25. fid = fopen(power4,'rb'); %open plans in Matlab
  26. planPower4 = reshape(fread(fid,'float'),72,72,dim);
  27. fclose(fid);
  28. tumorMaskLin = single(~~(full3D(prescTumorLin))); %define tissue masks
  29. tumorMaskRoot = single(~~(full3D(prescTumorRoot)));
  30. tumorMaskQuad = single(~~(full3D(prescTumorQuad)));
  31. tumorMaskPower025 = single(~~(full3D(prescTumorPower025)));
  32. tumorMaskPower4 = single(~~(full3D(prescTumorPower4)));
  33. tumorDoseLin245cm = planLin245cm.*tumorMaskLin; %dose distributions in each ROI
  34. tumorDoseRoot245cm = planRoot245cm.*tumorMaskRoot;
  35. tumorDoseQuad245cm = planQuad245cm.*tumorMaskQuad;
  36. tumorDosePower025 = planPower025.*tumorMaskPower025;
  37. tumorDosePower4 = planPower4.*tumorMaskPower4;
  38. %convert 3D matrices to sparse matrices
  39. tumorDoseSparseLin245cm = sparse3D(tumorDoseLin245cm);
  40. tumorDoseSparseRoot245cm = sparse3D(tumorDoseRoot245cm);
  41. tumorDoseSparseQuad245cm = sparse3D(tumorDoseQuad245cm);
  42. tumorDoseSparsePower025 = sparse3D(tumorDosePower025);
  43. tumorDoseSparsePower4 = sparse3D(tumorDosePower4);
  44. tumorNormDoseLin245cm = ((tumorDoseSparseLin245cm.non_zero_values)./(prescTumorLin.non_zero_values)); %calculate equivalent uniform dose
  45. tumorAverageDoseLin245cm = mean(prescTumorLin.non_zero_values);
  46. tumorEdoseLin245cm = tumorNormDoseLin245cm.*tumorAverageDoseLin245cm;
  47. tumorEdoseLin245cm = tumorNormDoseLin245cm;
  48. indLin = find(((tumorNormDoseLin245cm-1.).^2).^(1/2) >= 0.05);
  49. volLin = 100.-double(length(indLin))./(double(tumorDoseSparseLin245cm.Nind))*100
  50. tumorNormDoseRoot245cm = ((tumorDoseSparseRoot245cm.non_zero_values)./(prescTumorRoot.non_zero_values)); %calculate equivalent uniform dose
  51. tumorAverageDoseRoot245cm = mean(prescTumorRoot.non_zero_values);
  52. tumorEdoseRoot245cm = tumorNormDoseRoot245cm.*tumorAverageDoseRoot245cm;
  53. tumorEdoseRoot245cm = tumorNormDoseRoot245cm;
  54. indRoot = find(((tumorNormDoseRoot245cm-1.).^2).^(1/2) >= 0.05);
  55. volRoot = 100.-double(length(indRoot))./(double(tumorDoseSparseRoot245cm.Nind))*100
  56. tumorNormDoseQuad245cm = ((tumorDoseSparseQuad245cm.non_zero_values)./(prescTumorQuad.non_zero_values)); %calculate equivalent uniform dose
  57. tumorAverageDoseQuad245cm = mean(prescTumorQuad.non_zero_values);
  58. tumorEdoseQuad245cm = tumorNormDoseQuad245cm.*tumorAverageDoseQuad245cm;
  59. tumorEdoseQuad245cm = tumorNormDoseQuad245cm;
  60. indQuad = find(((tumorNormDoseQuad245cm-1.).^2).^(1/2) >= 0.05);
  61. volQuad = 100.-double(length(indQuad))./(double(tumorDoseSparseQuad245cm.Nind))*100
  62. tumorNormDosePower025 = ((tumorDoseSparsePower025.non_zero_values)./(prescTumorPower025.non_zero_values)); %calculate equivalent uniform dose
  63. tumorAverageDosePower025 = mean(prescTumorPower025.non_zero_values);
  64. tumorEdosePower025 = tumorNormDosePower025.*tumorAverageDosePower025;
  65. tumorEdosePower025 = tumorNormDosePower025;
  66. indPower025 = find(((tumorNormDosePower025-1.).^2).^(1/2) >= 0.15);
  67. min(((tumorNormDosePower025-1.).^2).^(1/2))
  68. volPower025 = 100.-double(length(indPower025))./(double(tumorDoseSparsePower025.Nind))*100
  69. tumorNormDosePower4 = ((tumorDoseSparsePower4.non_zero_values)./(prescTumorPower4.non_zero_values)); %calculate equivalent uniform dose
  70. tumorAverageDosePower4 = mean(prescTumorPower4.non_zero_values);
  71. tumorEdosePower4 = tumorNormDosePower4.*tumorAverageDosePower4;
  72. tumorEdosePower4 = tumorNormDosePower4;
  73. indPower4 = find(((tumorNormDosePower4-1.).^2).^(1/2) >= 0.15);
  74. volPower4 = 100.-double(length(indPower4))./(double(tumorDoseSparsePower4.Nind))*100
  75. %
  76. dMax = 1.1*max(tumorEdoseLin245cm);
  77. dvhbins = [0:Ndvhbins-1]*dMax/Ndvhbins;
  78. dHistTumorLin245cm = histc(tumorEdoseLin245cm,dvhbins);
  79. dHistTumorRoot245cm = histc(tumorEdoseRoot245cm,dvhbins);
  80. dHistTumorQuad245cm = histc(tumorEdoseQuad245cm,dvhbins);
  81. dHistTumorPower025 = histc(tumorEdosePower025,dvhbins);
  82. dHistTumorPower4 = histc(tumorEdosePower4,dvhbins);
  83. % dHistTumorGomp245cm = histc(tumorEdoseGomp245cm,dvhbins);
  84. %
  85. %
  86. % dHistTumorUni245cm = histc(tumorDoseSparseUni245cm.non_zero_values,dvhbins);
  87. % dHistTumorLin = histc(tumorDoseSparseLin245cm.non_zero_values,dvhbins);
  88. % dHistTumorRoot = histc(tumorDoseSparseRoot245cm.non_zero_values,dvhbins);
  89. % dHistTumorQuad = histc(tumorDoseSparseQuad245cm.non_zero_values,dvhbins);
  90. % dHistTumorGomp = histc(tumorDoseSparseGomp245cm.non_zero_values,dvhbins);
  91. %
  92. %
  93. %
  94. % dHistTumorNormUni245cm = 100*dHistTumorUni245cm./(double(tumorDoseSparseUni245cm.Nind));
  95. dHistTumorNormLin245cm = 100*dHistTumorLin245cm./(double(tumorDoseSparseLin245cm.Nind));
  96. dHistTumorNormRoot245cm = 100*dHistTumorRoot245cm./(double(tumorDoseSparseRoot245cm.Nind));
  97. dHistTumorNormQuad245cm = 100*dHistTumorQuad245cm./(double(tumorDoseSparseQuad245cm.Nind));
  98. dHistTumorNormPower025 = 100*dHistTumorPower025./(double(tumorDoseSparsePower025.Nind));
  99. dHistTumorNormPower4 = 100*dHistTumorPower4./(double(tumorDoseSparsePower4.Nind));
  100. % dHistTumorNormGomp245cm = 100*dHistTumorGomp245cm./(double(tumorDoseSparseGomp245cm.Nind));
  101. %
  102. % dHistTumorNormLin = 100*dHistTumorLin./(double(tumorDoseSparseLin245cm.Nind));
  103. % dHistTumorNormRoot = 100*dHistTumorRoot./(double(tumorDoseSparseRoot245cm.Nind));
  104. % dHistTumorNormQuad = 100*dHistTumorQuad./(double(tumorDoseSparseQuad245cm.Nind));
  105. % dHistTumorNormGomp = 100*dHistTumorGomp/(prescTumorGomp.Nind);
  106. %
  107. %
  108. % DVHtumorUni245cm = flipud(cumsum(flipud(dHistTumorNormUni245cm)));
  109. DVHtumorLin245cm = flipud(cumsum(flipud(dHistTumorNormLin245cm)));
  110. DVHtumorRoot245cm = flipud(cumsum(flipud(dHistTumorNormRoot245cm)));
  111. DVHtumorQuad245cm = flipud(cumsum(flipud(dHistTumorNormQuad245cm)));
  112. DVHtumorPower025 = flipud(cumsum(flipud(dHistTumorNormPower025)));
  113. DVHtumorPower4 = flipud(cumsum(flipud(dHistTumorNormPower4)));
  114. % DVHtumorGomp245cm = flipud(cumsum(flipud(dHistTumorNormGomp245cm)));
  115. %
  116. % DVHtumorLin = flipud(cumsum(flipud(dHistTumorNormLin)));
  117. % DVHtumorRoot = flipud(cumsum(flipud(dHistTumorNormRoot)));
  118. % DVHtumorQuad = flipud(cumsum(flipud(dHistTumorNormQuad)));
  119. % DVHtumorGomp = flipud(cumsum(flipud(dHistTumorNormGomp)));
  120. normalfigure(dvhbins,DVHtumorLin245cm,DVHtumorRoot245cm,DVHtumorQuad245cm,DVHtumorPower025,DVHtumorPower4)
  121. % DVHrtLin245cm,DVHrtRoot245cm,DVHrtQuad245cm,...
  122. % DVHcordLin245cm,DVHcordRoot245cm,DVHcordQuad245cm,...
  123. % DVHlarynxLin245cm,DVHlarynxRoot245cm,DVHlarynxQuad245cm)
  124. end