raytrace.cpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345
  1. /* C_RAYTRACE.C ***************************************************************/
  2. #include "defs.h"
  3. extern char errstr[200]; // error string that all routines have access to
  4. int raytrace(FLOAT_GRID *electron_density_grid,FLOAT_GRID *radiological_depth_grid,
  5. POINT point1,POINT point2)
  6. /*
  7. --------------------------------------------------------------------------------
  8. NAME
  9. c_raytrace
  10. SYNOPSIS
  11. point1 and point2 are the end points of the ray.
  12. DESCRIPTION
  13. This function traces the ray from point x to point y (in real
  14. coords), assigning depth to any voxels which that ray crosses. The
  15. technique used is that described by Siddon R.L. (Med Phys 12 (2),
  16. 1985). This routine will not be understood without thorough reading
  17. of that paper! Point1 and point2 are the start and end points of the
  18. ray, respectively. External structures of type GRID are assumed to
  19. exist, where electron_density_grid are the electron densities, and
  20. radiological_depth_grid is the output grid for these calculations.
  21. Voxels in radiological_depth_grid are initially set -ve prior to
  22. calling this function.
  23. AUTHOR
  24. Written by David C. Murray
  25. University of Waikato
  26. Private Bag 3105
  27. Hamilton
  28. New Zealand
  29. and Copyright (1991) to
  30. David C. Murray and Peter W. Hoban,
  31. Cancer Society of New Zealand Inc., and
  32. University of Waikato.
  33. --------------------------------------------------------------------------------
  34. */
  35. {
  36. /* Variable descriptions:
  37. x1,x2,y1,y2,z1,z2 are the coordinates of the ray end points
  38. (i.e. (x1,y1,z1)=source, (x2,y2,z2)=point beyond phantom)
  39. xp1,yp1,zp1 are the real coords of voxel region origin (in cm)
  40. Nx,Ny,Nz are (no. of voxels + 1) in each direction
  41. dx,dy,dz are the widths in cm of the voxels
  42. */
  43. float x1,y1,z1,
  44. x2,y2,z2,
  45. xp1,yp1,zp1,
  46. dx,dy,dz;
  47. int Nx,Ny,Nz;
  48. /*General ray-trace algorithm variables*/
  49. float xpN, ypN, zpN; /*real coords in cm of region limits*/
  50. float alpha_x_min,alpha_y_min,alpha_z_min,alpha_x_max,alpha_y_max,alpha_z_max;
  51. /*limits of alpha x,y,z parameters*/
  52. float alpha_min, alpha_max; /*limits of alpha parameter*/
  53. int i_min,i_max,j_min,j_max,k_min,k_max;/*limits of indices for x,y,z dirns*/
  54. float *alpha_x,*alpha_y,*alpha_z; /*hold sets of x,y,z alpha values*/
  55. float *alpha; /*holds merged set of alpha values*/
  56. int i_index,j_index,k_index; /*loop indices for merging alphas*/
  57. int a; /*loop counter*/
  58. int max_index; /*max index of merged alpha array*/
  59. float d12; /*distance between ray end points*/
  60. float alpha_mid; /*mid-point of intersection length*/
  61. float length; /*intersection length*/
  62. int i,j,k; /*indices of voxel with int. length*/
  63. float rpl = 0.0; /*radiological path length in cm*/
  64. float voxel_density; /*temporary voxel density*/
  65. float lmax; // best possible penetration pathlength for a voxel
  66. float pnorm; // absolute difference between p1 and p2
  67. /* Assign variables */
  68. /******************************************************************************/
  69. x1 = point1.x;
  70. y1 = point1.y;
  71. z1 = point1.z;
  72. x2 = point2.x;
  73. y2 = point2.y;
  74. z2 = point2.z;
  75. Nx = electron_density_grid->x_count + 1;
  76. Ny = electron_density_grid->y_count + 1;
  77. Nz = electron_density_grid->z_count + 1;
  78. dx = electron_density_grid->inc.x;
  79. dy = electron_density_grid->inc.y;
  80. dz = electron_density_grid->inc.z;
  81. // (xp1,yp1,zp1) are the locations of the first grid planes for each dimension
  82. xp1 = electron_density_grid->start.x - 0.5*dx;
  83. yp1 = electron_density_grid->start.y - 0.5*dy;
  84. zp1 = electron_density_grid->start.z - 0.5*dz;
  85. pnorm = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
  86. // this is the largest pathlength possible for a ray through a voxel:
  87. lmax = sqrt(pow((x2-x1)*dx,2.0f)+pow((y2-y1)*dy,2.0f)+pow((z2-z1)*dz,2.0f))/pnorm;
  88. /* Calculate xpN,ypN,zpN */
  89. /******************************************************************************/
  90. xpN = xp1 + (Nx-1)*dx;
  91. ypN = yp1 + (Ny-1)*dy;
  92. zpN = zp1 + (Nz-1)*dz;
  93. /*Calculate alpha_min and alpha_max*/
  94. /******************************************************************************/
  95. /*Avoid division by zero*/
  96. if (x1==x2)
  97. x2 += 0.00001;
  98. if (y1==y2)
  99. y2 += 0.00001;
  100. if (z1==z2)
  101. z2 += 0.00001;
  102. if ((fabs(x1-x2)<dx) && (fabs(y1-y2)<dy) && (fabs(z1-z2)<dz))
  103. {
  104. sprintf(errstr,"Error - ray trace region too small.");
  105. return(FAILURE);
  106. }
  107. alpha_x_min = (((xp1-x1)/(x2-x1))<((xpN-x1)/(x2-x1))) ? ((xp1-x1)/(x2-x1))
  108. : ((xpN-x1)/(x2-x1));
  109. alpha_y_min = (((yp1-y1)/(y2-y1))<((ypN-y1)/(y2-y1))) ? ((yp1-y1)/(y2-y1))
  110. : ((ypN-y1)/(y2-y1));
  111. alpha_z_min = (((zp1-z1)/(z2-z1))<((zpN-z1)/(z2-z1))) ? ((zp1-z1)/(z2-z1))
  112. : ((zpN-z1)/(z2-z1));
  113. alpha_x_max = (((xp1-x1)/(x2-x1))>((xpN-x1)/(x2-x1))) ? ((xp1-x1)/(x2-x1))
  114. : ((xpN-x1)/(x2-x1));
  115. alpha_y_max = (((yp1-y1)/(y2-y1))>((ypN-y1)/(y2-y1))) ? ((yp1-y1)/(y2-y1))
  116. : ((ypN-y1)/(y2-y1));
  117. alpha_z_max = (((zp1-z1)/(z2-z1))>((zpN-z1)/(z2-z1))) ? ((zp1-z1)/(z2-z1))
  118. : ((zpN-z1)/(z2-z1));
  119. alpha_min = (alpha_x_min>alpha_y_min) ? alpha_x_min : alpha_y_min;
  120. if (alpha_z_min>alpha_min)
  121. alpha_min = alpha_z_min;
  122. if (alpha_min<0)
  123. alpha_min = 0;
  124. alpha_max = (alpha_x_max<alpha_y_max) ? alpha_x_max : alpha_y_max;
  125. if (alpha_z_max<alpha_max)
  126. alpha_max = alpha_z_max;
  127. if (alpha_max>1)
  128. alpha_max = 1;
  129. // Monitor lines...
  130. /*
  131. printf(" alpha_x,y,z_min: %7.4f %7.4f %7.4f\n",
  132. alpha_x_min,alpha_y_min,alpha_z_min);
  133. printf(" alpha_x,y,z_max: %7.4f %7.4f %7.4f\n",
  134. alpha_x_max,alpha_y_max,alpha_z_max);
  135. printf(" alpha_min,alpha_max: %7.4f %7.4f\n",alpha_min,alpha_max);
  136. printf("Nx, xpN, x2,x1, dx = %d %5.2f %5.2f %5.2f %5.2f\n",Nx,xpN,x2,x1,dx); */
  137. /*Determine the ranges of i,j,k indices*/
  138. /******************************************************************************/
  139. /*The following assignments require conversion from float to integer types*/
  140. /*The value 0.001 is added/subtracted to ensure that the ceiling and floor*/
  141. /*functions convert to the correct value. Note that the range of these*/
  142. /*variables is from 1 to Nx,Ny,Nz, NOT 0 to Nx-1,Ny-1,Nz-1*/
  143. i_min = (x2>x1) ? (int) ceil((float) Nx - (xpN-alpha_min*(x2-x1)-x1)/dx-0.001)
  144. : (int) ceil((float) Nx - (xpN-alpha_max*(x2-x1)-x1)/dx-0.001);
  145. i_max = (x2>x1) ? (int) floor(1.0000 + (x1+alpha_max*(x2-x1)-xp1)/dx+0.001)
  146. : (int) floor(1.0000 + (x1+alpha_min*(x2-x1)-xp1)/dx+0.001);
  147. j_min = (y2>y1) ? (int) ceil((float) Ny - (ypN-alpha_min*(y2-y1)-y1)/dy-0.001)
  148. : (int) ceil((float) Ny - (ypN-alpha_max*(y2-y1)-y1)/dy-0.001);
  149. j_max = (y2>y1) ? (int) floor(1.0000 + (y1+alpha_max*(y2-y1)-yp1)/dy+0.001)
  150. : (int) floor(1.0000 + (y1+alpha_min*(y2-y1)-yp1)/dy+0.001);
  151. k_min = (z2>z1) ? (int) ceil((float) Nz - (zpN-alpha_min*(z2-z1)-z1)/dz-0.001)
  152. : (int) ceil((float) Nz - (zpN-alpha_max*(z2-z1)-z1)/dz-0.001);
  153. k_max = (z2>z1) ? (int) floor(1.0000 + (z1+alpha_max*(z2-z1)-zp1)/dz+0.001)
  154. : (int) floor(1.0000 + (z1+alpha_min*(z2-z1)-zp1)/dz+0.001);
  155. /*Monitor lines...
  156. fprintf(stdout," i,j,k_min: %3d %3d %3d\n",i_min,j_min,k_min);
  157. fprintf(stdout," i,j,k_max: %3d %3d %3d\n",i_max,j_max,k_max);
  158. */
  159. /*Generate sets of alpha values,reversing order if necessary*/
  160. /******************************************************************************/
  161. /*allocate array space on stack*/
  162. if ((alpha_x = (float*) calloc(Nx+1,sizeof(float))) == NULL)
  163. {
  164. sprintf(errstr,"Error - insufficient heap for alpha_x allocation.");
  165. return(FAILURE);
  166. }
  167. if ((alpha_y = (float*) calloc(Ny+1,sizeof(float))) == NULL)
  168. {
  169. sprintf(errstr,"Error - insufficient heap for alpha_y allocation.");
  170. return(FAILURE);
  171. }
  172. if ((alpha_z = (float*) calloc(Nz+1,sizeof(float))) == NULL)
  173. {
  174. sprintf(errstr,"Error - insufficient heap for alpha_z allocation.");
  175. return(FAILURE);
  176. }
  177. /*
  178. printf("Nx = %d, i_min = %d, i_max = %d\n",Nx,i_min,i_max);
  179. printf("Ny = %d, j_min = %d, j_max = %d\n",Ny,j_min,j_max);
  180. printf("Nz = %d, k_min = %d, k_max = %d\n",Nz,k_min,k_max); */
  181. if (i_min <= i_max)
  182. if (x2>x1)
  183. {
  184. alpha_x[0] = ((xp1+(i_min-1)*dx)-x1)/(x2-x1);
  185. for (a=1;a<=i_max-i_min;a++)
  186. alpha_x[a] = alpha_x[a-1]+dx/(x2-x1);
  187. }
  188. else
  189. {
  190. alpha_x[i_max-i_min] = ((xp1+(i_min-1)*dx)-x1)/(x2-x1);
  191. for (a=i_max-i_min-1;a>=0;a--)
  192. alpha_x[a] = alpha_x[a+1]+(dx/(x2-x1));
  193. }
  194. alpha_x[i_max-i_min+1] = 10000.0;
  195. if (j_min <= j_max)
  196. if (y2>y1)
  197. {
  198. alpha_y[0] = ((yp1+(j_min-1)*dy)-y1)/(y2-y1);
  199. for (a=1;a<=j_max-j_min;a++)
  200. alpha_y[a] = alpha_y[a-1]+dy/(y2-y1);
  201. }
  202. else
  203. {
  204. alpha_y[j_max-j_min] = ((yp1+(j_min-1)*dy)-y1)/(y2-y1);
  205. for (a=j_max-j_min-1;a>=0;a--)
  206. alpha_y[a] = alpha_y[a+1]+(dy/(y2-y1));
  207. }
  208. alpha_y[j_max-j_min+1] = 10001.0;
  209. if (k_min <= k_max)
  210. if (z2>z1)
  211. {
  212. alpha_z[0] = ((zp1+(k_min-1)*dz)-z1)/(z2-z1);
  213. for (a=1;a<=k_max-k_min;a++)
  214. alpha_z[a] = alpha_z[a-1]+(dz/(z2-z1));
  215. }
  216. else
  217. {
  218. alpha_z[k_max-k_min] = ((zp1+(k_min-1)*dz)-z1)/(z2-z1);
  219. for (a=k_max-k_min-1;a>=0;a--)
  220. alpha_z[a] = alpha_z[a+1]+(dz/(z2-z1));
  221. }
  222. alpha_z[k_max-k_min+1] = 10002.0;
  223. /*Monitor lines...
  224. if (i_max<i_min)
  225. fprintf(stdout," No alpha_x values\n");
  226. else
  227. fprintf(stdout," First & last alpha_x values: %7.4f %7.4f\n",
  228. alpha_x[0],alpha_x[i_max-i_min]);
  229. if (j_max<j_min)
  230. fprintf(stdout," No alpha_y values\n");
  231. else
  232. fprintf(stdout," First & last alpha_y values: %7.4f %7.4f\n",
  233. alpha_y[0],alpha_y[j_max-j_min]);
  234. if (k_max<k_min)
  235. fprintf(stdout," No alpha_z values\n");
  236. else
  237. fprintf(stdout," First & last alpha_z values: %7.4f %7.4f\n",
  238. alpha_z[0],alpha_z[k_max-k_min]);
  239. */
  240. /*Generate merged set of alpha values*/
  241. /******************************************************************************/
  242. if ((alpha = (float*) calloc(Nx+Ny+Nz+3,sizeof(float))) == NULL)
  243. {
  244. sprintf(errstr,"Error - insufficient heap for alpha allocation.");
  245. return(FAILURE);
  246. }
  247. max_index = (i_max-i_min+1)+(j_max-j_min+1)+(k_max-k_min+1)+1;
  248. alpha[0] = alpha_min;
  249. i_index = 0;
  250. j_index = 0;
  251. k_index = 0;
  252. for (a=1;a<=max_index-1;a++)
  253. if (alpha_x[i_index]<alpha_y[j_index])
  254. if (alpha_x[i_index]<alpha_z[k_index])
  255. {
  256. alpha[a] = alpha_x[i_index];
  257. i_index += 1;
  258. }
  259. else
  260. {
  261. alpha[a] = alpha_z[k_index];
  262. k_index += 1;
  263. }
  264. else
  265. if (alpha_y[j_index]<alpha_z[k_index])
  266. {
  267. alpha[a] = alpha_y[j_index];
  268. j_index += 1;
  269. }
  270. else
  271. {
  272. alpha[a] = alpha_z[k_index];
  273. k_index += 1;
  274. }
  275. alpha[max_index] = alpha_max;
  276. free(alpha_x); //deallocate temp array storage
  277. free(alpha_y);
  278. free(alpha_z);
  279. /*Monitor lines...
  280. fprintf(stdout," Number of elements in merged set = %4d\n",max_index+1);
  281. for (a=0;a<=max_index;a++)
  282. fprintf(stdout," Element %3d = %7.5f\n",a,alpha[a]);
  283. */
  284. /*Calculate voxel lengths and indices, and assign radiological depth*/
  285. /******************************************************************************/
  286. d12 = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
  287. //d12 is distance between ray end pts
  288. // printf("made it this far in raytrace.\n");
  289. for (a=1;a<=max_index;a++)
  290. {
  291. length = d12*(alpha[a]-alpha[a-1]); //length is voxel intersection length
  292. if (fabs(length)>0.01) //do not process unless > 0.01 cm
  293. {
  294. alpha_mid = (alpha[a]+alpha[a-1])/2.0;
  295. //alpha_mid is middle of int. length
  296. i = (int) floor((x1 + alpha_mid*(x2-x1) - xp1)/dx);
  297. j = (int) floor((y1 + alpha_mid*(y2-y1) - yp1)/dy);
  298. k = (int) floor((z1 + alpha_mid*(z2-z1) - zp1)/dz);
  299. //i,j,k are indices of voxel
  300. // Remember that this function traces only a single ray.
  301. // rpl has been set to zero during initialisation.
  302. voxel_density = GRID_VALUE(electron_density_grid,i,j,k);
  303. rpl += length * voxel_density/2.0; // add first half of int. length
  304. // store pathlength only if the voxel is intersected almost directly
  305. // by the ray
  306. if (length>=0.75/2*lmax && GRID_VALUE(radiological_depth_grid,i,j,k)<0.0)
  307. GRID_VALUE(radiological_depth_grid, i, j, k) = rpl;
  308. rpl += length * voxel_density/2.0; //add second half of int. length
  309. }
  310. }
  311. free(alpha); //deallocate remaining array storage
  312. return(SUCCESS);
  313. } /*End of s_raytrace routine*/