1
0

terma_dose_masks.cpp 10.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301
  1. /* terma_dose_masks.cpp */
  2. /* All voxels that lie within a cylinder of radius RAPERTURE + RSAFE, centered
  3. on the central beam axis, are marked in the dose_mask grid. These voxels will
  4. be used in the dose calculation. The voxels that are marked in the terma_mask
  5. grid are a subset of those that are marked in the dose_mask grid. The terma_mask
  6. grid contains values that specify the fraction of each voxel that lies inside
  7. of the beam aperture in the beam's eye view plane, which exists at a distance
  8. SAD from the x-ray source location. */
  9. /* Determination of voxel "insideness" for the terma_mask:
  10. If it is discovered that a voxel is intersected by an aperture boundary in the
  11. beam's eye view plane, then that voxel is broken-up into Mus x Nus x Qus sub-
  12. voxels (upsampled). The fraction of the sub-voxels that lie inside the beam
  13. aperture in the beam's eye view then becomes the insideness for that voxel. */
  14. #include "defs.h"
  15. extern char errstr[200]; // error string that all routines have access to
  16. int terma_dose_masks(FLOAT_GRID *terma_mask, FLOAT_GRID *dose_mask, BEAM *bm)
  17. // The user must allocate space for terma_mask and dose_mask before calling
  18. // this function.
  19. // The validity of all of the arguments to this function are assumed to be
  20. // checked with the parse_func, which acts as a gatekeeper for the convolution
  21. // program.
  22. {
  23. // scalars
  24. int i,j,k,m; // indices
  25. int ius,jus,kus; // upsample voxel indices
  26. int M,N,Q; // grid dimensions
  27. float dx,dy,dz; // voxel dimensions
  28. float xp,yp,SAD; // aperture center in ip,jp,kp coordinate system
  29. float xsub,ysub,zsub; // coordinates of sub-voxels
  30. float del_xp,del_yp; // width of beamlet field in plane at SAD
  31. float Rmax,Rcurr,Rcyl,Rvox,Rcyl_sqr,Rproj;
  32. float eta,gamma,delta; // dummy variables
  33. // vectors
  34. float *x,*y,*z; // x,y,z coordinate vectors
  35. float *dxus,*dyus,*dzus; // upsample vectors
  36. float *y_vec,*ip,*jp,*kp,*d_hat; // beam vectors
  37. // matrices
  38. float **c; // matrix of the 8 corners of voxel grid
  39. // Ensure that upsample factors are all greater than zero
  40. if (Mus < 1 || Nus < 1 || Qus < 1)
  41. {
  42. sprintf(errstr,"Upsample factors must all be greater than zero.");
  43. return(FAILURE);
  44. }
  45. // record the sizes of terma_mask and dose_mask
  46. M = terma_mask->x_count;
  47. N = terma_mask->y_count;
  48. Q = terma_mask->z_count;
  49. if (M != dose_mask->x_count || N != dose_mask->y_count
  50. || Q != dose_mask->z_count)
  51. {
  52. sprintf(errstr,"dose_mask and terma_mask dimensions are incompatible.");
  53. return(FAILURE);
  54. }
  55. dx = terma_mask->inc.x;
  56. dy = terma_mask->inc.y;
  57. dz = terma_mask->inc.z;
  58. // initialize vectors
  59. x = fvector(0,M-1);
  60. y = fvector(0,N-1);
  61. z = fvector(0,Q-1);
  62. dxus = fvector(0,Mus-1);
  63. dyus = fvector(0,Nus-1);
  64. dzus = fvector(0,Qus-1);
  65. y_vec = fvector(0,2);
  66. ip = fvector(0,2);
  67. jp = fvector(0,2);
  68. kp = fvector(0,2);
  69. d_hat = fvector(0,2);
  70. // initialize matrices
  71. c = fmatrix(0,7,0,2);
  72. // fill-in the voxel coordinate vectors
  73. for (i=0;i<M;i++) x[i] = terma_mask->start.x + (float)i*dx;
  74. for (j=0;j<N;j++) y[j] = terma_mask->start.y + (float)j*dy;
  75. for (k=0;k<Q;k++) z[k] = terma_mask->start.z + (float)k*dz;
  76. // fill-in the upsample vectors
  77. // For a voxel at (i,j,k), the xyz coordinates of a sub-voxel (ius,jus,kus)
  78. // are (x[i] + dxus[ius], y[j] + dyus[jus], z[k] + dzus[k])
  79. dxus[0] = -0.5*dx*(1-1/Mus);
  80. for (ius=1;ius<Mus;ius++) dxus[ius] = dxus[ius-1] + dx/Mus;
  81. dyus[0] = -0.5*dy*(1-1/Nus);
  82. for (jus=1;jus<Nus;jus++) dyus[jus] = dyus[jus-1] + dy/Nus;
  83. dzus[0] = -0.5*dz*(1-1/Qus);
  84. for (kus=1;kus<Qus;kus++) dzus[kus] = dzus[kus-1] + dz/Qus;
  85. // rewrite beam vectors for quick reference later
  86. y_vec[0] = bm->y_vec[0];
  87. y_vec[1] = bm->y_vec[1];
  88. y_vec[2] = bm->y_vec[2];
  89. // coordinate system in beam's eye view
  90. for (j=0;j<3;j++) ip[j] = bm->ip[j];
  91. for (j=0;j<3;j++) jp[j] = bm->jp[j];
  92. for (j=0;j<3;j++) kp[j] = bm->kp[j];
  93. // aperture center in beam's eye view
  94. xp = bm->xp;
  95. yp = bm->yp;
  96. // aperture size in beam's eye view
  97. del_xp = bm->del_xp;
  98. del_yp = bm->del_yp;
  99. // source-axis distance of the beam
  100. SAD = bm->SAD;
  101. // calculate the max distance between the source vector and each grid corner
  102. Rmax = 0.0;
  103. // the the matrix c with the coordinates of the eight corners of the grid.
  104. for (i=0;i<=1;i++)
  105. for (j=0;j<=1;j++)
  106. for (k=0;k<=1;k++)
  107. {
  108. // 0 => lower corners 1 => upper corners
  109. c[i+j*2+k*4][0] = (1-(float)i)*(x[0]-dx/2.0) + (float)i*(x[M-1]+dx/2.0);
  110. c[i+j*2+k*4][1] = (1-(float)j)*(y[0]-dy/2.0) + (float)j*(y[N-1]+dy/2.0);
  111. c[i+j*2+k*4][2] = (1-(float)k)*(z[0]-dz/2.0) + (float)k*(z[Q-1]+dz/2.0);
  112. }
  113. // find which corner is the farthest from the source
  114. for (m=0;m<=7;m++)
  115. {
  116. Rcurr = sqrt(pow(y_vec[0] - c[m][0],2.0f)
  117. + pow(y_vec[1] - c[m][1],2.0f)
  118. + pow(y_vec[2] - c[m][2],2.0f));
  119. if (Rcurr > Rmax)
  120. Rmax = Rcurr;
  121. }
  122. // Fill the dose_mask
  123. // Radius of cylinder about y_vec + kp that will contain all voxels
  124. // to be used in the dose calculation.
  125. Rcyl = sqrt(del_xp*del_xp + del_yp*del_yp)*Rmax/SAD; // project Rmax to aperture plane
  126. Rcyl_sqr = Rcyl*Rcyl + RSAFE*RSAFE; // add an empirical safety margin to the radius squared
  127. // calculate the true source direction
  128. // d_hat points at the center of the aperture from the source location
  129. d_hat[0] = xp*ip[0] + yp*jp[0] + SAD*kp[0];
  130. d_hat[1] = xp*ip[1] + yp*jp[1] + SAD*kp[1];
  131. d_hat[2] = xp*ip[2] + yp*jp[2] + SAD*kp[2];
  132. // normalize d_hat so it's a true "hat" vector
  133. delta = sqrt(d_hat[0]*d_hat[0] + d_hat[1]*d_hat[1] + d_hat[2]*d_hat[2]);
  134. d_hat[0] = d_hat[0]/delta;
  135. d_hat[1] = d_hat[1]/delta;
  136. d_hat[2] = d_hat[2]/delta;
  137. for (k=0;k<Q;k++)
  138. for (j=0;j<N;j++)
  139. for (i=0;i<M;i++)
  140. {
  141. // squared distance between the voxel and the source:
  142. eta = (x[i] - y_vec[0])*(x[i] - y_vec[0])
  143. + (y[j] - y_vec[1])*(y[j] - y_vec[1])
  144. + (z[k] - y_vec[2])*(z[k] - y_vec[2]);
  145. // squared distance between the voxel and the source along source direction,
  146. // y_vec + kp:
  147. gamma = pow((x[i] - y_vec[0])*d_hat[0]
  148. + (y[j] - y_vec[1])*d_hat[1]
  149. + (z[k] - y_vec[2])*d_hat[2],2.0f);
  150. // printf("%lf %lf\n",eta, gamma);
  151. // squared difference between the voxel and the axis of the cylinder
  152. delta = eta - gamma;
  153. // If the voxel is inside the cylinder about y_vec + ip of radius Rcyl plus
  154. // a safety margin then mark it in dose_mask.
  155. if (delta <= Rcyl_sqr)
  156. GRID_VALUE(dose_mask,i,j,k) = 1.0;
  157. else
  158. GRID_VALUE(dose_mask,i,j,k) = 0.0;
  159. }
  160. // Fill the terma_mask, including the insideness of each voxel
  161. // Each voxel is enclosed with in a sphere of the following radius:
  162. Rvox = 0.5*sqrt(dx*dx + dy*dy + dz*dz);
  163. for (i=0;i<M;i++)
  164. for (j=0;j<N;j++)
  165. for (k=0;k<Q;k++)
  166. // deal only with voxels that are marked in the dose_mask
  167. if (GRID_VALUE(dose_mask,i,j,k) > 0.0)
  168. {
  169. // distance between source and voxel in the kp direction
  170. delta = kp[0]*(x[i] - y_vec[0])
  171. + kp[1]*(y[j] - y_vec[1])
  172. + kp[2]*(z[k] - y_vec[2]);
  173. // voxel's projected offset on the aperture plane in the ip direction:
  174. eta = (SAD/delta)*( ip[0]*(x[i] - y_vec[0])
  175. + ip[1]*(y[j] - y_vec[1])
  176. + ip[2]*(z[k] - y_vec[2])) - xp;
  177. // voxel's projected offset on the aperture plane in the jp direction:
  178. gamma = (SAD/delta)*( jp[0]*(x[i] - y_vec[0])
  179. + jp[1]*(y[j] - y_vec[1])
  180. + jp[2]*(z[k] - y_vec[2])) - yp;
  181. // take absolute value of offsets
  182. eta = fabs(eta);
  183. gamma = fabs(gamma);
  184. Rproj = Rvox*SAD/delta; // voxel radius projected to plane at SAD
  185. // Determine where the voxel lies with respect to the aperture:
  186. if (eta <= 0.5*del_xp+Rproj && gamma <= 0.5*del_yp+Rproj)
  187. // voxel is inside aperture plus a half-voxel margin:
  188. if (eta >= 0.5*del_xp-Rproj || gamma >= 0.5*del_yp-Rproj)
  189. // voxel is between the aperture plus/minus a half-voxel margin:
  190. // (true at this point if the voxel size is larger than the aperture size)
  191. {
  192. // Determine insideness of the voxel by breaking it up
  193. // into subvoxels and projecting the parts to the aperture plane.
  194. m = 0; // number of subvoxels inside aperture
  195. // project each subvoxel onto the aperture at SAD
  196. for (ius=0;ius<Mus;ius++)
  197. for (jus=0;jus<Nus;jus++)
  198. for (kus=0;kus<Qus;kus++)
  199. {
  200. // find the center of the subvoxel
  201. xsub = x[i] + dxus[ius];
  202. ysub = y[j] + dyus[jus];
  203. zsub = z[k] + dzus[kus];
  204. // project the subvoxel onto the aperture
  205. // distance between source and subvoxel in the kp direction
  206. delta = kp[0]*(xsub - y_vec[0])
  207. + kp[1]*(ysub - y_vec[1])
  208. + kp[2]*(zsub - y_vec[2]);
  209. // projected offset on aperture plane in the ip direction:
  210. eta = (SAD/delta)*( ip[0]*(xsub - y_vec[0])
  211. + ip[1]*(ysub - y_vec[1])
  212. + ip[2]*(zsub - y_vec[2])) - xp;
  213. // projected offset on aperture plane in the jp direction:
  214. gamma = (SAD/delta)*( jp[0]*(xsub - y_vec[0])
  215. + jp[1]*(ysub - y_vec[1])
  216. + jp[2]*(zsub - y_vec[2])) - yp;
  217. eta = fabs(eta);
  218. gamma = fabs(gamma);
  219. // check if the subvoxel is inside the aperture at SAD
  220. if (eta <= 0.5*del_xp && gamma <= 0.5*del_yp)
  221. m++;
  222. }
  223. // the fraction of subvoxels inside the aperture becomes the insidness
  224. GRID_VALUE(terma_mask,i,j,k) = (float)m/(float)(Mus*Nus*Qus);
  225. }
  226. else
  227. // voxel is inside the aperture minus a half-voxel margin
  228. GRID_VALUE(terma_mask,i,j,k) = 1.0;
  229. else
  230. // voxel outside the aperture plus the half-voxel margin
  231. GRID_VALUE(terma_mask,i,j,k) = 0.0;
  232. }
  233. // free vectors
  234. free_fvector(x,0,M-1);
  235. free_fvector(y,0,N-1);
  236. free_fvector(z,0,Q-1);
  237. free_fvector(dxus,0,Mus-1);
  238. free_fvector(dyus,0,Nus-1);
  239. free_fvector(dzus,0,Qus-1);
  240. free_fvector(y_vec,0,2);
  241. free_fvector(ip,0,2);
  242. free_fvector(jp,0,2);
  243. free_fvector(kp,0,2);
  244. free_fvector(d_hat,0,2);
  245. // free matrices
  246. free_fmatrix(c,0,7,0,2);
  247. return(SUCCESS);
  248. }