singleRaytraceClean.cpp 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563
  1. /* singleRaytraceClean.cpp */
  2. /* Does a raytracing operation between point1 and point2 in a density grid
  3. defined by DENS, START, and INC. The visited voxels and their corresponding
  4. radiological depths are returned as the INDVISITED and DEFFVISITED vectors,
  5. respectively */
  6. #include "mex.h"
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include <malloc.h>
  10. #include <math.h>
  11. #include <string.h>
  12. #define SUCCESS 1
  13. #define FAILURE 0
  14. // macro for returning 3D grid index
  15. #define GRIDIND3D(x,y,z,dim) ((x) + dim[0]*(y) + dim[0]*dim[1]*(z))
  16. // 3D point
  17. typedef struct
  18. {
  19. float x;
  20. float y;
  21. float z;
  22. } POINT;
  23. // 3D grid
  24. typedef struct
  25. {
  26. POINT start;
  27. POINT inc;
  28. int x_count;
  29. int y_count;
  30. int z_count;
  31. float *matrix;
  32. } FLOAT_GRID;
  33. // macro to access grid values
  34. #define GRID_VALUE(GRID_ptr, i, j, k)\
  35. ((GRID_ptr)->matrix[(i) +\
  36. (GRID_ptr)->x_count *\
  37. ((j) + ((k) * (GRID_ptr)->y_count))])
  38. // User inputs
  39. #define DENS (prhs[0]) // electron density grid
  40. #define START (prhs[1]) // 3-element start vector
  41. #define INC (prhs[2]) // 3-element voxel size vector
  42. #define POINT1 (prhs[3]) // starting point for raytrace
  43. #define POINT2 (prhs[4]) // ending point for raytrace
  44. #define INDVISITED (plhs[0]) // indices of visited voxels
  45. #define DEFFVISITED (plhs[1]) // deff values of visited voxels
  46. // Prototypes
  47. int raytrace(FLOAT_GRID *,POINT,POINT,int *, float *, int *);
  48. double round(double);
  49. char errstr[200]; // error string that all routines have access to
  50. void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  51. {
  52. // scalars
  53. int k,M,N,Q;
  54. // pointers to vectors that do not need to be freed
  55. const int *densDim,*startDim,*incDim,*p1Dim,*p2Dim;
  56. float *densPtr,*startPtr,*incPtr,*p1Ptr,*p2Ptr;
  57. int *indVisitedPtr;
  58. float *deffVisitedPtr;
  59. int dims[2];
  60. int *indVisited;
  61. float *deffVisited; // indices and deffs visited along ray
  62. int maxVisited;
  63. int Nvisited;
  64. // point structures that will be used for raytrace operation
  65. POINT p1;
  66. POINT p2;
  67. // grid structures for raytracing operation
  68. FLOAT_GRID dens;
  69. // get dimensions of the input mxArrays
  70. densDim = mxGetDimensions(DENS);
  71. startDim = mxGetDimensions(START);
  72. incDim = mxGetDimensions(INC);
  73. p1Dim = mxGetDimensions(POINT1);
  74. p2Dim = mxGetDimensions(POINT2);
  75. // Check validity of input arguments.
  76. if ( ((mxGetNumberOfDimensions(DENS) != 2) && (mxGetNumberOfDimensions(DENS) != 3))
  77. || (mxGetClassID(DENS) != mxSINGLE_CLASS))
  78. mexErrMsgTxt("DENS must be a 3-D float array.");
  79. if ( (mxGetNumberOfDimensions(START) != 2)
  80. || (mxGetClassID(START) != mxSINGLE_CLASS)
  81. || (startDim[0]*startDim[1] != 3))
  82. mexErrMsgTxt("START must be a 3-element float vector.");
  83. if ( (mxGetNumberOfDimensions(INC) != 2)
  84. || (mxGetClassID(INC) != mxSINGLE_CLASS)
  85. || (incDim[0]*incDim[1] != 3))
  86. mexErrMsgTxt("INC must be a 3-element float vector.");
  87. if ( (mxGetNumberOfDimensions(POINT1) != 2)
  88. || (mxGetClassID(POINT1) != mxSINGLE_CLASS)
  89. || (p1Dim[0]*p1Dim[1] != 3))
  90. mexErrMsgTxt("POINT1 must be a 3-element float vector.");
  91. if ( (mxGetNumberOfDimensions(POINT2) != 2)
  92. || (mxGetClassID(POINT2) != mxSINGLE_CLASS)
  93. || (p2Dim[0]*p2Dim[1] != 3))
  94. mexErrMsgTxt("POINT2 must be a 3-element float vector.");
  95. // assign data pointers
  96. densPtr = (float *)mxGetData(DENS);
  97. startPtr = (float *)mxGetData(START);
  98. incPtr = (float *)mxGetData(INC);
  99. p1Ptr = (float *)mxGetData(POINT1);
  100. p2Ptr = (float *)mxGetData(POINT2);
  101. // get grid dimensions
  102. M = densDim[0];
  103. N = densDim[1];
  104. if (mxGetNumberOfDimensions(DENS) == 3)
  105. Q = densDim[2];
  106. else // handle 2D data as 3D data with the third dimension having 1 increment
  107. Q = 1;
  108. // printf("(M,N,Q) = (%d,%d,%d)\n",M,N,Q);
  109. // the maximum number of voxels that can be visited along a single line
  110. maxVisited = (int)(2*sqrt((double)(M*M+N*N+Q*Q)) + 1);
  111. indVisited = (int *)malloc(sizeof(int)*maxVisited);
  112. deffVisited = (float *)malloc(sizeof(int)*maxVisited);
  113. for (k=0;k<maxVisited;k++)
  114. {
  115. indVisited[k] = 0;
  116. deffVisited[k] = 0.0;
  117. }
  118. // Assign dens and deff parameters to structures for input to raytrace
  119. // Don't forget, dens.matrix and densPtr point to the same data, and
  120. // deff.matrix and deffPtr point to the same data.
  121. dens.matrix = densPtr;
  122. dens.start.x = startPtr[0];
  123. dens.start.y = startPtr[1];
  124. dens.start.z = startPtr[2];
  125. dens.inc.x = incPtr[0];
  126. dens.inc.y = incPtr[1];
  127. dens.inc.z = incPtr[2];
  128. dens.x_count = M;
  129. dens.y_count = N;
  130. dens.z_count = Q;
  131. // assign values to p1 and p2
  132. p1.x = p1Ptr[0];
  133. p1.y = p1Ptr[1];
  134. p1.z = p1Ptr[2];
  135. p2.x = p2Ptr[0];
  136. p2.y = p2Ptr[1];
  137. p2.z = p2Ptr[2];
  138. // do the raytrace, filling in voxels that are passed on the way
  139. if (raytrace(&dens,p1,p2,indVisited,deffVisited,&Nvisited) == FAILURE)
  140. mexErrMsgTxt(errstr);
  141. dims[0] = 1;
  142. dims[1] = Nvisited;
  143. // dims[1] = 1;
  144. // set up output vectors
  145. INDVISITED = mxCreateNumericArray(2,dims,mxINT32_CLASS,mxREAL);
  146. DEFFVISITED = mxCreateNumericArray(2,dims,mxSINGLE_CLASS,mxREAL);
  147. indVisitedPtr = (int *)mxGetData(INDVISITED);
  148. deffVisitedPtr = (float *)mxGetData(DEFFVISITED);
  149. // copy results array to Matlab outputs
  150. for (k=0;k<Nvisited;k++)
  151. {
  152. indVisitedPtr[k] = indVisited[k];
  153. deffVisitedPtr[k] = deffVisited[k];
  154. }
  155. }
  156. extern char errstr[200]; // error string that all routines have access to
  157. int raytrace(FLOAT_GRID *electron_density_grid,POINT point1,POINT point2,
  158. int *indVisited, float *deffVisited, int *Nvisited)
  159. /*
  160. --------------------------------------------------------------------------------
  161. NAME
  162. c_raytrace
  163. SYNOPSIS
  164. point1 and point2 are the end points of the ray.
  165. DESCRIPTION
  166. This function traces the ray from point x to point y (in real
  167. coords), assigning depth to any voxels which that ray crosses. The
  168. technique used is that described by Siddon R.L. (Med Phys 12 (2),
  169. 1985). This routine will not be understood without thorough reading
  170. of that paper! Point1 and point2 are the start and end points of the
  171. ray, respectively. External structures of type GRID are assumed to
  172. exist, where electron_density_grid are the electron densities, and
  173. radiological_depth_grid is the output grid for these calculations.
  174. Voxels in radiological_depth_grid are initially set -ve prior to
  175. calling this function.
  176. AUTHOR
  177. Written by David C. Murray
  178. University of Waikato
  179. Private Bag 3105
  180. Hamilton
  181. New Zealand
  182. and Copyright (1991) to
  183. David C. Murray and Peter W. Hoban,
  184. Cancer Society of New Zealand Inc., and
  185. University of Waikato.
  186. --------------------------------------------------------------------------------
  187. */
  188. {
  189. /* Variable descriptions:
  190. x1,x2,y1,y2,z1,z2 are the coordinates of the ray end points
  191. (i.e. (x1,y1,z1)=source, (x2,y2,z2)=point beyond phantom)
  192. xp1,yp1,zp1 are the real coords of voxel region origin (in cm)
  193. Nx,Ny,Nz are (no. of voxels + 1) in each direction
  194. dx,dy,dz are the widths in cm of the voxels
  195. */
  196. float x1,y1,z1,
  197. x2,y2,z2,
  198. xp1,yp1,zp1,
  199. dx,dy,dz;
  200. int Nx,Ny,Nz;
  201. /*General ray-trace algorithm variables*/
  202. float xpN, ypN, zpN; /*real coords in cm of region limits*/
  203. float alpha_x_min,alpha_y_min,alpha_z_min,alpha_x_max,alpha_y_max,alpha_z_max;
  204. /*limits of alpha x,y,z parameters*/
  205. float alpha_min, alpha_max; /*limits of alpha parameter*/
  206. int i_min,i_max,j_min,j_max,k_min,k_max;/*limits of indices for x,y,z dirns*/
  207. float *alpha_x,*alpha_y,*alpha_z; /*hold sets of x,y,z alpha values*/
  208. float *alpha; /*holds merged set of alpha values*/
  209. int i_index,j_index,k_index; /*loop indices for merging alphas*/
  210. int a; /*loop counter*/
  211. int max_index; /*max index of merged alpha array*/
  212. float d12; /*distance between ray end points*/
  213. float alpha_mid; /*mid-point of intersection length*/
  214. float length; /*intersection length*/
  215. int i,j,k; /*indices of voxel with int. length*/
  216. float rpl = 0.0; /*radiological path length in cm*/
  217. float voxel_density; /*temporary voxel density*/
  218. float lmax; // best possible penetration pathlength for a voxel
  219. float pnorm; // absolute difference between p1 and p2
  220. /* Assign variables */
  221. /******************************************************************************/
  222. x1 = point1.x;
  223. y1 = point1.y;
  224. z1 = point1.z;
  225. x2 = point2.x;
  226. y2 = point2.y;
  227. z2 = point2.z;
  228. Nx = electron_density_grid->x_count + 1;
  229. Ny = electron_density_grid->y_count + 1;
  230. Nz = electron_density_grid->z_count + 1;
  231. dx = electron_density_grid->inc.x;
  232. dy = electron_density_grid->inc.y;
  233. dz = electron_density_grid->inc.z;
  234. // (xp1,yp1,zp1) are the locations of the first grid planes for each dimension
  235. xp1 = electron_density_grid->start.x - (float)0.5*dx;
  236. yp1 = electron_density_grid->start.y - (float)0.5*dy;
  237. zp1 = electron_density_grid->start.z - (float)0.5*dz;
  238. pnorm = (float)sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
  239. // this is the largest pathlength possible for a ray through a voxel:
  240. lmax = (float)sqrt(pow((x2-x1)*dx,2.0f)+pow((y2-y1)*dy,2.0f)+pow((z2-z1)*dz,2.0f))/pnorm;
  241. /* Calculate xpN,ypN,zpN */
  242. /******************************************************************************/
  243. xpN = xp1 + (Nx-1)*dx;
  244. ypN = yp1 + (Ny-1)*dy;
  245. zpN = zp1 + (Nz-1)*dz;
  246. /*Calculate alpha_min and alpha_max*/
  247. /******************************************************************************/
  248. /*Avoid division by zero*/
  249. if (x1==x2)
  250. x2 += 0.00001;
  251. if (y1==y2)
  252. y2 += 0.00001;
  253. if (z1==z2)
  254. z2 += 0.00001;
  255. if ((fabs(x1-x2)<dx) && (fabs(y1-y2)<dy) && (fabs(z1-z2)<dz))
  256. {
  257. sprintf(errstr,"Error - ray trace region too small.");
  258. return(FAILURE);
  259. }
  260. alpha_x_min = (((xp1-x1)/(x2-x1))<((xpN-x1)/(x2-x1))) ? ((xp1-x1)/(x2-x1))
  261. : ((xpN-x1)/(x2-x1));
  262. alpha_y_min = (((yp1-y1)/(y2-y1))<((ypN-y1)/(y2-y1))) ? ((yp1-y1)/(y2-y1))
  263. : ((ypN-y1)/(y2-y1));
  264. alpha_z_min = (((zp1-z1)/(z2-z1))<((zpN-z1)/(z2-z1))) ? ((zp1-z1)/(z2-z1))
  265. : ((zpN-z1)/(z2-z1));
  266. alpha_x_max = (((xp1-x1)/(x2-x1))>((xpN-x1)/(x2-x1))) ? ((xp1-x1)/(x2-x1))
  267. : ((xpN-x1)/(x2-x1));
  268. alpha_y_max = (((yp1-y1)/(y2-y1))>((ypN-y1)/(y2-y1))) ? ((yp1-y1)/(y2-y1))
  269. : ((ypN-y1)/(y2-y1));
  270. alpha_z_max = (((zp1-z1)/(z2-z1))>((zpN-z1)/(z2-z1))) ? ((zp1-z1)/(z2-z1))
  271. : ((zpN-z1)/(z2-z1));
  272. alpha_min = (alpha_x_min>alpha_y_min) ? alpha_x_min : alpha_y_min;
  273. if (alpha_z_min>alpha_min)
  274. alpha_min = alpha_z_min;
  275. if (alpha_min<0)
  276. alpha_min = 0;
  277. alpha_max = (alpha_x_max<alpha_y_max) ? alpha_x_max : alpha_y_max;
  278. if (alpha_z_max<alpha_max)
  279. alpha_max = alpha_z_max;
  280. if (alpha_max>1)
  281. alpha_max = 1;
  282. // Monitor lines...
  283. /*
  284. printf(" alpha_x,y,z_min: %7.4f %7.4f %7.4f\n",
  285. alpha_x_min,alpha_y_min,alpha_z_min);
  286. printf(" alpha_x,y,z_max: %7.4f %7.4f %7.4f\n",
  287. alpha_x_max,alpha_y_max,alpha_z_max);
  288. printf(" alpha_min,alpha_max: %7.4f %7.4f\n",alpha_min,alpha_max);
  289. printf("Nx, xpN, x2,x1, dx = %d %5.2f %5.2f %5.2f %5.2f\n",Nx,xpN,x2,x1,dx); */
  290. /*Determine the ranges of i,j,k indices*/
  291. /******************************************************************************/
  292. /*The following assignments require conversion from float to integer types*/
  293. /*The value 0.001 is added/subtracted to ensure that the ceiling and floor*/
  294. /*functions convert to the correct value. Note that the range of these*/
  295. /*variables is from 1 to Nx,Ny,Nz, NOT 0 to Nx-1,Ny-1,Nz-1*/
  296. i_min = (x2>x1) ? (int) ceil((float) Nx - (xpN-alpha_min*(x2-x1)-x1)/dx-0.001)
  297. : (int) ceil((float) Nx - (xpN-alpha_max*(x2-x1)-x1)/dx-0.001);
  298. i_max = (x2>x1) ? (int) floor(1.0000 + (x1+alpha_max*(x2-x1)-xp1)/dx+0.001)
  299. : (int) floor(1.0000 + (x1+alpha_min*(x2-x1)-xp1)/dx+0.001);
  300. j_min = (y2>y1) ? (int) ceil((float) Ny - (ypN-alpha_min*(y2-y1)-y1)/dy-0.001)
  301. : (int) ceil((float) Ny - (ypN-alpha_max*(y2-y1)-y1)/dy-0.001);
  302. j_max = (y2>y1) ? (int) floor(1.0000 + (y1+alpha_max*(y2-y1)-yp1)/dy+0.001)
  303. : (int) floor(1.0000 + (y1+alpha_min*(y2-y1)-yp1)/dy+0.001);
  304. k_min = (z2>z1) ? (int) ceil((float) Nz - (zpN-alpha_min*(z2-z1)-z1)/dz-0.001)
  305. : (int) ceil((float) Nz - (zpN-alpha_max*(z2-z1)-z1)/dz-0.001);
  306. k_max = (z2>z1) ? (int) floor(1.0000 + (z1+alpha_max*(z2-z1)-zp1)/dz+0.001)
  307. : (int) floor(1.0000 + (z1+alpha_min*(z2-z1)-zp1)/dz+0.001);
  308. // Monitor lines...
  309. /* printf(" i,j,k_min: %3d %3d %3d\n",i_min,j_min,k_min);
  310. printf(" i,j,k_max: %3d %3d %3d\n",i_max,j_max,k_max); */
  311. // return if the indices do not make sense
  312. if ( i_max - i_min + 1 < 0 || i_max - i_min + 1 > Nx
  313. || j_max - j_min + 1 < 0 || j_max - j_min + 1 > Ny
  314. || k_max - k_min + 1 < 0 || k_max - k_min + 1 > Nz)
  315. {
  316. *Nvisited = 0; // no voxels were visited
  317. return(SUCCESS);
  318. }
  319. /*Generate sets of alpha values,reversing order if necessary*/
  320. /******************************************************************************/
  321. /*allocate array space on stack*/
  322. if ((alpha_x = (float*) calloc(Nx+1,sizeof(float))) == NULL)
  323. {
  324. sprintf(errstr,"Error - insufficient heap for alpha_x allocation.");
  325. return(FAILURE);
  326. }
  327. if ((alpha_y = (float*) calloc(Ny+1,sizeof(float))) == NULL)
  328. {
  329. sprintf(errstr,"Error - insufficient heap for alpha_y allocation.");
  330. return(FAILURE);
  331. }
  332. if ((alpha_z = (float*) calloc(Nz+1,sizeof(float))) == NULL)
  333. {
  334. sprintf(errstr,"Error - insufficient heap for alpha_z allocation.");
  335. return(FAILURE);
  336. }
  337. /*
  338. printf("Nx = %d, i_min = %d, i_max = %d\n",Nx,i_min,i_max);
  339. printf("Ny = %d, j_min = %d, j_max = %d\n",Ny,j_min,j_max);
  340. printf("Nz = %d, k_min = %d, k_max = %d\n",Nz,k_min,k_max); */
  341. if (i_min <= i_max)
  342. if (x2>x1)
  343. {
  344. alpha_x[0] = ((xp1+(i_min-1)*dx)-x1)/(x2-x1);
  345. for (a=1;a<=i_max-i_min;a++)
  346. alpha_x[a] = alpha_x[a-1]+dx/(x2-x1);
  347. }
  348. else
  349. {
  350. alpha_x[i_max-i_min] = ((xp1+(i_min-1)*dx)-x1)/(x2-x1);
  351. for (a=i_max-i_min-1;a>=0;a--)
  352. alpha_x[a] = alpha_x[a+1]+(dx/(x2-x1));
  353. }
  354. alpha_x[i_max-i_min+1] = 10000.0;
  355. if (j_min <= j_max)
  356. if (y2>y1)
  357. {
  358. alpha_y[0] = ((yp1+(j_min-1)*dy)-y1)/(y2-y1);
  359. for (a=1;a<=j_max-j_min;a++)
  360. alpha_y[a] = alpha_y[a-1]+dy/(y2-y1);
  361. }
  362. else
  363. {
  364. alpha_y[j_max-j_min] = ((yp1+(j_min-1)*dy)-y1)/(y2-y1);
  365. for (a=j_max-j_min-1;a>=0;a--)
  366. alpha_y[a] = alpha_y[a+1]+(dy/(y2-y1));
  367. }
  368. alpha_y[j_max-j_min+1] = 10001.0;
  369. if (k_min <= k_max)
  370. if (z2>z1)
  371. {
  372. alpha_z[0] = ((zp1+(k_min-1)*dz)-z1)/(z2-z1);
  373. for (a=1;a<=k_max-k_min;a++)
  374. alpha_z[a] = alpha_z[a-1]+(dz/(z2-z1));
  375. }
  376. else
  377. {
  378. alpha_z[k_max-k_min] = ((zp1+(k_min-1)*dz)-z1)/(z2-z1);
  379. for (a=k_max-k_min-1;a>=0;a--)
  380. alpha_z[a] = alpha_z[a+1]+(dz/(z2-z1));
  381. }
  382. alpha_z[k_max-k_min+1] = 10002.0;
  383. /* Monitor lines...
  384. if (i_max<i_min)
  385. fprintf(stdout," No alpha_x values\n");
  386. else
  387. fprintf(stdout," First & last alpha_x values: %7.4f %7.4f\n",
  388. alpha_x[0],alpha_x[i_max-i_min]);
  389. if (j_max<j_min)
  390. fprintf(stdout," No alpha_y values\n");
  391. else
  392. fprintf(stdout," First & last alpha_y values: %7.4f %7.4f\n",
  393. alpha_y[0],alpha_y[j_max-j_min]);
  394. if (k_max<k_min)
  395. fprintf(stdout," No alpha_z values\n");
  396. else
  397. fprintf(stdout," First & last alpha_z values: %7.4f %7.4f\n",
  398. alpha_z[0],alpha_z[k_max-k_min]);
  399. */
  400. /*Generate merged set of alpha values*/
  401. /******************************************************************************/
  402. if ((alpha = (float*) calloc(Nx+Ny+Nz+3,sizeof(float))) == NULL)
  403. {
  404. sprintf(errstr,"Error - insufficient heap for alpha allocation.");
  405. return(FAILURE);
  406. }
  407. max_index = (i_max-i_min+1)+(j_max-j_min+1)+(k_max-k_min+1)+1;
  408. alpha[0] = alpha_min;
  409. i_index = 0;
  410. j_index = 0;
  411. k_index = 0;
  412. for (a=1;a<=max_index-1;a++)
  413. if (alpha_x[i_index]<alpha_y[j_index])
  414. if (alpha_x[i_index]<alpha_z[k_index])
  415. {
  416. alpha[a] = alpha_x[i_index];
  417. i_index += 1;
  418. }
  419. else
  420. {
  421. alpha[a] = alpha_z[k_index];
  422. k_index += 1;
  423. }
  424. else
  425. if (alpha_y[j_index]<alpha_z[k_index])
  426. {
  427. alpha[a] = alpha_y[j_index];
  428. j_index += 1;
  429. }
  430. else
  431. {
  432. alpha[a] = alpha_z[k_index];
  433. k_index += 1;
  434. }
  435. alpha[max_index] = alpha_max;
  436. free(alpha_x); //deallocate temp array storage
  437. free(alpha_y);
  438. free(alpha_z);
  439. /*Monitor lines...
  440. fprintf(stdout," Number of elements in merged set = %4d\n",max_index+1);
  441. for (a=0;a<=max_index;a++)
  442. fprintf(stdout," Element %3d = %7.5f\n",a,alpha[a]);
  443. */
  444. /*Calculate voxel lengths and indices, and assign radiological depth*/
  445. /******************************************************************************/
  446. d12 = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
  447. //d12 is distance between ray end pts
  448. *Nvisited = 0; // start counter for number of voxel visited
  449. for (a=1;a<=max_index;a++)
  450. {
  451. length = d12*(alpha[a]-alpha[a-1]); //length is voxel intersection length
  452. if (fabs(length)>0.01) //do not process unless > 0.01 cm
  453. {
  454. alpha_mid = (alpha[a]+alpha[a-1])/2.0; //alpha_mid is middle of int. length
  455. //i,j,k are indices of voxel
  456. i = (int) round((x1 + alpha_mid*(x2-x1) - xp1)/dx - 0.5);
  457. j = (int) round((y1 + alpha_mid*(y2-y1) - yp1)/dy - 0.5);
  458. k = (int) round((z1 + alpha_mid*(z2-z1) - zp1)/dz - 0.5);
  459. // Remember that this function traces only a single ray.
  460. // rpl has been set to zero during initialisation.
  461. if (i>=0 && i<Nx-1 && j>=0 && j<Ny-1 && k>=0 && k<Nz-1)
  462. {
  463. voxel_density = GRID_VALUE(electron_density_grid,i,j,k);
  464. rpl += length * voxel_density/2.0; // add first half of int. length
  465. // store pathlength only if the voxel is intersected almost directly
  466. // by the ray
  467. // Include the visited indices in Matlab notation, since that's
  468. // how the data will eventually be returned.
  469. indVisited[*Nvisited] = i + j*electron_density_grid->x_count
  470. + k*electron_density_grid->x_count*electron_density_grid->y_count+1;
  471. deffVisited[*Nvisited] = rpl; // radiological depth at visited location
  472. (*Nvisited)++;
  473. rpl += length * voxel_density/2.0; //add second half of int. length
  474. }
  475. }
  476. }
  477. free(alpha); //deallocate remaining array storage
  478. return(SUCCESS);
  479. } /*End of s_raytrace routine*/
  480. double round(double x)
  481. // Rounds x to the nearest integer.
  482. {
  483. static double xfl;
  484. xfl = floor((double)x);
  485. if ( ((double)x - xfl) >= 0.5)
  486. return(xfl + 1.0); // round up
  487. else
  488. return(xfl); // round down
  489. }