dens2mstp.m 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166
  1. function [Smw Fmw2] = dens2mstp(dens_ratio)
  2. % DENS2MSTP Convert physical density to mass stopping power and Fmw2
  3. %
  4. % Usage:
  5. % [Smw Fmw2] = dens2mstp ( density_matrix )
  6. % INPUT:
  7. % density_matrix = density matrix (g/cm^3) with any dimensionality
  8. % OUTPUT:
  9. % Smw = stopping power ratio
  10. % Fmw2 = squared F value ratio to water
  11. %
  12. % Given an array of densities relative to water (assumed to be from a CT
  13. % scan), the mass stopping power ratio (relative to water) array is
  14. % calculated. The conversion is done by assuming that every voxel in the CT
  15. % consists of either air, adipose tissue, muscle, or bone. A lookup table
  16. % of mass stopping power ratios between each material and water is then
  17. % constructed by computing the average mass stopping power ratio over some
  18. % proton energy region, usually between 50 MeV and 250 MeV. The mass
  19. % stopping powers are looked-up in tables obtained from the NIST website.
  20. %
  21. % This code was intended to be used as a way of converting CT
  22. % data to mass stopping power data for proton dose calculations. The energy
  23. % of the proton beam is assumed to be between 50 MeV and 250 MeV, a regime
  24. % in which the approximation that stopping power ratios are energy
  25. % independent holds.
  26. %
  27. % Note that the value returned is the stopping power ratio. The stopping power
  28. % has units of MeV/cm. This is an important distinction!
  29. %
  30. % Copyleft: Ryan Flynn, Xiaohu Mo
  31. % Version: 1.2
  32. % Set up the boundary points for the density-to-material lookup table.
  33. air_low = 0;
  34. air_high = 0.05; fat_low = air_high;
  35. fat_high = 0.99; wat_low = fat_high;
  36. wat_high = 1.01; mus_low = wat_high;
  37. mus_high = 1.25; bon_low = mus_high;
  38. bon_high = 20.0;
  39. % check the input argument
  40. if ~isempty(find(dens_ratio < air_low, 1))
  41. error('Cannot have negative densities. Redefine input argument.');
  42. end
  43. if ~isempty(find(dens_ratio > bon_high, 1))
  44. error('Extremely high densities are present in the input array, data may not be scaled properly');
  45. end
  46. % Product of Fmw and rho for tabluated materials
  47. fmw_air_to_water = 990.8 * 1.205E-3;
  48. fmw_fat_to_water = 0.972 * 0.92;
  49. fmw_mus_to_water = 0.972 * 1.04;
  50. fmw_bon_to_water = 0.656 * 1.85;
  51. fmw_wat_to_water = 1.0;
  52. % Stopping power ranges that will be used to calculate the average stopping
  53. % power ratios between each material and water.
  54. Tlo = 1; % lower energy boundary in MeV
  55. Thi = 250; % upper energy boundary in MeV
  56. % load the stopping powers
  57. load('MSTPliquid_water.mat');
  58. water = MSTP;
  59. load('MSTPair.mat');
  60. air = MSTP;
  61. load('MSTPadipose.mat');
  62. fat = MSTP;
  63. load('MSTPskeletal_muscle.mat');
  64. mus = MSTP;
  65. load('MSTPcompact_bone.mat');
  66. bon = MSTP;
  67. % calculate the spacing between energy bins
  68. dT = [water.T(1) diff(water.T)];
  69. % find the indices of integration that will be used
  70. indices_low = find(water.T <= Tlo);
  71. ind_low = indices_low(end);
  72. indices_high = find(water.T >= Thi);
  73. ind_high = indices_high(1);
  74. % Calculate the average stopping power ratio for the material and water
  75. % between the energy range specified above.
  76. smw_air_to_water = sum(air.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
  77. smw_fat_to_water = sum(fat.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
  78. smw_wat_to_water = 1.0;
  79. smw_mus_to_water = sum(mus.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
  80. smw_bon_to_water = sum(bon.tot(ind_low:ind_high)./water.tot(ind_low:ind_high).*dT(ind_low:ind_high))/(Thi-Tlo);
  81. % material mask
  82. region_air = find (dens_ratio >= air_low & dens_ratio < air_high);
  83. region_fat = find (dens_ratio >= fat_low & dens_ratio < fat_high);
  84. region_wat = find (dens_ratio >= wat_low & dens_ratio < wat_high);
  85. region_mus = find (dens_ratio >= mus_low & dens_ratio < mus_high);
  86. region_bon = find (dens_ratio >= bon_low);
  87. % fill the mass stopping power ratio
  88. Smw = zeros ( size ( dens_ratio ) );
  89. Smw ( region_air ) = smw_air_to_water * dens_ratio ( region_air );
  90. Smw ( region_fat ) = smw_fat_to_water * dens_ratio ( region_fat );
  91. Smw ( region_wat ) = smw_wat_to_water * dens_ratio ( region_wat );
  92. Smw ( region_mus ) = smw_mus_to_water * dens_ratio ( region_mus );
  93. Smw ( region_bon ) = smw_bon_to_water * dens_ratio ( region_bon );
  94. % fill the Fmw2
  95. % divide by density to appropriately scale Fmw
  96. Fmw2 = zeros ( size ( dens_ratio ) );
  97. Fmw2 ( region_air ) = fmw_air_to_water ./ dens_ratio ( region_air );
  98. Fmw2 ( region_fat ) = fmw_fat_to_water ./ dens_ratio ( region_fat );
  99. Fmw2 ( region_wat ) = fmw_wat_to_water ./ dens_ratio ( region_wat );
  100. Fmw2 ( region_mus ) = fmw_mus_to_water ./ dens_ratio ( region_mus );
  101. Fmw2 ( region_bon ) = fmw_bon_to_water ./ dens_ratio ( region_bon );
  102. Fmw2 = Fmw2.^2;
  103. % % plot the stopping power ratio as a function of energy for bone and water
  104. % f1 = figure; set(gcf,'color',[1 1 1]);
  105. % p = plot(water.T(ind_low:ind_high),bone.tot(ind_low:ind_high)./water.tot(ind_low:ind_high));
  106. % set(p,'LineWidth',2);
  107. % set(gca,'YLim',[0.9 0.95]);
  108. % set(gca,'FontSize',14);
  109. % title('Mass stopping power ratio between bone and water','FontSize',14);
  110. % xlabel('Energy (MeV)','FontSize',14);
  111. % ylabel('Ratio','FontSize',14);
  112. % print -dbitmap MSTP_ratio_with_energy.png;
  113. %
  114. % % plot the density ratio to mass stopping power ratio lookup table:
  115. % xmax = 2.5;
  116. % N = 50;
  117. % dx = xmax/N;
  118. % x = [0:N-1]*dx;
  119. % y = zeros(size(x));
  120. % y(x >= air_low & x <= air_high) = air_to_water;
  121. % y(x > adipose_low & x <= adipose_high) = adipose_to_water;
  122. % y(x > muscle_low & x <= muscle_high) = muscle_to_water;
  123. % y(x > bone_low) = bone_to_water;
  124. %
  125. % z = y.*x;
  126. %
  127. % f2 = figure; set(gcf,'color',[1 1 1]);
  128. % plot(x(x >= air_low & x <= air_high),y(x >= air_low & x <= air_high),'.b',...
  129. % x(x > adipose_low & x <= adipose_high),y(x > adipose_low & x <= adipose_high),'+b',...
  130. % x(x > muscle_low & x <= muscle_high),y(x > muscle_low & x <= muscle_high),'xb',...
  131. % x(x > bone_low),y(x > bone_low),'*b');
  132. % l = legend('air','adipose','muscle','bone');
  133. % set(l,'FontSize',14);
  134. % set(gca,'FontSize',14);
  135. % xlabel('Density ratio','FontSize',14);
  136. % ylabel('Mass stopping power ratio','FontSize',14);
  137. % title('Change in mass stopping power ratio with density','FontSize',14);
  138. % print -dbitmap MSTP_ratio_with_density.png;
  139. %
  140. % f3 = figure; set(gcf,'color',[1 1 1]);
  141. % plot(x(x >= air_low & x <= air_high),z(x >= air_low & x <= air_high),'.b',...
  142. % x(x > adipose_low & x <= adipose_high),z(x > adipose_low & x <= adipose_high),'+b',...
  143. % x(x > muscle_low & x <= muscle_high),z(x > muscle_low & x <= muscle_high),'xb',...
  144. % x(x > bone_low),z(x > bone_low),'*b');
  145. % l = legend('air','adipose','muscle','bone');
  146. % set(l,'FontSize',14,'position',[0.6258 0.2433 0.2250 0.2460]);
  147. % set(gca,'FontSize',14);
  148. % xlabel('Density ratio','FontSize',14);
  149. % ylabel('Stopping power ratio','FontSize',14);
  150. % title('Change in stopping power ratio with density','FontSize',14);
  151. % print -dbitmap STP_ratio_with_density.png;