1
0

terma_kerma.cpp 1.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960
  1. /* terma_kerma.cpp */
  2. /* Calculates the terma and the kerma of every voxel in the terma_mask. */
  3. #include "defs.h"
  4. extern char errstr[200]; // error string that all routines have access to
  5. int terma_kerma(FLOAT_GRID *deff, FLOAT_GRID *terma, FLOAT_GRID *kermac,
  6. MONO_KERNELS *mono, FLOAT_GRID *terma_mask)
  7. {
  8. // insideness is accounted for with terma_mask
  9. int i, j, k, e;
  10. float kermac0, terma0;
  11. //calculate T and Kc at zero depth for use in hardening correction
  12. //see Hoban et al 1994 (PMB)
  13. kermac0 = terma0 = 0.0;
  14. for (e=0;e<mono->nkernels;e++)
  15. {
  16. kermac0 = mono->fluence[e]*mono->energy[e]*mono->mu_en[e];
  17. terma0 = mono->fluence[e]*mono->energy[e]*mono->mu[e];
  18. }
  19. for (j=0;j<deff->y_count;j++)
  20. for (k=0;k<deff->z_count;k++)
  21. for (i=0;i<deff->x_count;i++)
  22. if (GRID_VALUE(terma_mask,i,j,k) > 0)
  23. {
  24. // The amount of each voxel that is inside the field (insideness) was
  25. // was accounted in the calculation of the terma_mask
  26. //sum terma and collision kerma over each energy in spectrum
  27. // (stored in mono kernel structure)
  28. for (e=0;e<mono->nkernels;e++)
  29. {
  30. GRID_VALUE(terma,i,j,k) += mono->fluence[e]*mono->energy[e]*mono->mu[e]
  31. * exp(-mono->mu[e]*GRID_VALUE(deff,i,j,k));
  32. GRID_VALUE(kermac,i,j,k) += mono->fluence[e]*mono->energy[e]*mono->mu_en[e]
  33. * exp(-mono->mu[e]*GRID_VALUE(deff,i,j,k));
  34. }
  35. //adjust terma and collision kerma according to insideness
  36. GRID_VALUE(terma,i,j,k) *= GRID_VALUE(terma_mask,i,j,k);
  37. GRID_VALUE(kermac,i,j,k) *= GRID_VALUE(terma_mask,i,j,k);
  38. // beam hardening correction
  39. if (terma0 <= 0.0 || kermac0 <= 0.0)
  40. nrerror("Input spectrum must not sum to zero.");
  41. else
  42. GRID_VALUE(terma,i,j,k) *= (GRID_VALUE(kermac,i,j,k)/GRID_VALUE(terma,i,j,k))
  43. /(kermac0/terma0);
  44. }
  45. else
  46. {
  47. GRID_VALUE(terma,i,j,k) = 0.0;
  48. GRID_VALUE(kermac,i,j,k) = 0.0;
  49. }
  50. return(SUCCESS);
  51. }