convolutionCondor.cpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390
  1. /* Cconvolution.cpp */
  2. #include "defs.h"
  3. //function prototypes
  4. int load_kernels(MONO_KERNELS *, char []);
  5. int load_geometry(FLOAT_GRID *, char []);
  6. int pop_beam(BEAM *, FILE *);
  7. int calc_deff(FLOAT_GRID *,FLOAT_GRID *, FLOAT_GRID *, BEAM *);
  8. int terma_kerma(FLOAT_GRID *,FLOAT_GRID *,FLOAT_GRID *,MONO_KERNELS *,FLOAT_GRID *);
  9. int make_poly_kernel(MONO_KERNELS *, KERNEL *);
  10. int calc_dose(FLOAT_GRID *,FLOAT_GRID *,FLOAT_GRID *,KERNEL *,BEAM *, FLOAT_GRID *);
  11. int terma_dose_masks(FLOAT_GRID *, FLOAT_GRID *, BEAM *);
  12. int convolution(MONO_KERNELS *, BEAM *, FLOAT_GRID *, FILE *);
  13. char errstr[200]; // error string that all routines have access to
  14. int main(int argc, char *argv[])
  15. // Expects four input arguments:
  16. // 1) kernel_filenames -- contains a list of the kernel files
  17. // 2) geometry_filenames -- contains a list of the geometry files
  18. // 3) beamspec batch file -- locations of beamspec files in batch
  19. // 4) beamlet batch file -- locations of resulting beamlets in batch
  20. {
  21. int i,j,b,B;
  22. char tmpstr[200];
  23. FLOAT_GRID density;
  24. MONO_KERNELS mono_kernels;
  25. BEAM beam;
  26. FILE *beamspec_batch_file;
  27. FILE *beamlet_batch_file;
  28. /* // print the arguments
  29. printf("Input arguments:\n");
  30. for (j=1;j<argc;j++)
  31. printf("%s\n",argv[j]); */
  32. if (argc != 5)
  33. {
  34. printf("Expecting four input command line arguments, received %d.\n",argc);
  35. return(FAILURE);
  36. }
  37. if (load_kernels(&mono_kernels,argv[1]) == FAILURE)
  38. {
  39. sprintf(tmpstr,"Failed at loading the kernels.\n");
  40. strcat(tmpstr,errstr);
  41. strcpy(errstr,tmpstr);
  42. printf("%s\n",errstr);
  43. return(FAILURE);
  44. }
  45. if (load_geometry(&density,argv[2]) == FAILURE)
  46. {
  47. sprintf(tmpstr,"Failed at loading the geometry.\n");
  48. strcat(tmpstr,errstr);
  49. strcpy(errstr,tmpstr);
  50. printf("%s\n",errstr);
  51. return(FAILURE);
  52. }
  53. /*
  54. // diagnostic lines
  55. printf("SAD = %lf \n",beam.SAD);
  56. printf("xp = %lf \n",beam.xp);
  57. printf("yp = %lf \n",beam.yp);
  58. printf("del_xp = %lf \n",beam.del_xp);
  59. printf("del_yp = %lf \n",beam.del_yp);
  60. printf("y_vec = (%lf,%lf,%lf) \n",beam.y_vec[0],beam.y_vec[1],beam.y_vec[2]);
  61. printf("ip = (%lf,%lf,%lf) \n",beam.ip[0],beam.ip[1],beam.ip[2]);
  62. printf("jp = (%lf,%lf,%lf) \n",beam.jp[0],beam.jp[1],beam.jp[2]);
  63. printf("kp = (%lf,%lf,%lf) \n",beam.kp[0],beam.kp[1],beam.kp[2]);
  64. printf("Xcount = %d, Ycount = %d, Zcount = %d \n",density.x_count,density.y_count,density.z_count);
  65. printf("start = (%lf,%lf,%lf) \n",density.start.x,density.start.y,density.start.z);
  66. printf("inc = (%lf,%lf,%lf) \n",density.inc.x,density.inc.y,density.inc.z); */
  67. // open the beam specification batch file
  68. if ((beamspec_batch_file = fopen(argv[3],"r")) == NULL)
  69. {
  70. printf("Failed to open beamspec batch file %s\n",argv[3]);
  71. return(FAILURE);
  72. }
  73. // open the dose batch file
  74. if ((beamlet_batch_file = fopen(argv[4],"wb")) == NULL)
  75. {
  76. printf("Failed to open beamlet batch file %s\n",argv[4]);
  77. return(FAILURE);
  78. }
  79. // get the number of beams from the beamspec batch file
  80. if (fgets(tmpstr,100,beamspec_batch_file) == NULL) // pop off the first line
  81. {
  82. sprintf(errstr,"Could not read from beam data file.");
  83. printf("%s\n",errstr);
  84. return(FAILURE);
  85. }
  86. if (fscanf(beamspec_batch_file,"%d\n",&B) != 1)
  87. {
  88. sprintf(errstr,"Could not read-in number of beams from beamspec file.");
  89. printf("%s\n",errstr);
  90. return(FAILURE);
  91. }
  92. // write the number of beamlets in this batch as the first entry
  93. fwrite(&B,sizeof(int),1,beamlet_batch_file);
  94. // Do convolution calculations for all consistent beamspec and dose
  95. // filenames. If a calculation for a beamlet fails, print an error
  96. // and move on to the next beamspec file.
  97. for (b=0;b<B;b++)
  98. {
  99. // pop off a beam
  100. if(pop_beam(&beam,beamspec_batch_file) == FAILURE)
  101. {
  102. sprintf(tmpstr,"Failed to load beamspec number %d:\n",b);
  103. strcat(tmpstr,errstr);
  104. strcpy(errstr,tmpstr);
  105. printf("%s\n",errstr);
  106. }
  107. else if (convolution(&mono_kernels,&beam,&density,beamlet_batch_file) == FAILURE)
  108. // An error occurred, so print the error string to standard out
  109. // but do not terminate the remaining beam batches.
  110. {
  111. j = 0;
  112. fwrite(&beam.num,sizeof(int),1,beamlet_batch_file);
  113. fwrite(&density.x_count,sizeof(int),1,beamlet_batch_file);
  114. fwrite(&density.y_count,sizeof(int),1,beamlet_batch_file);
  115. fwrite(&density.z_count,sizeof(int),1,beamlet_batch_file);
  116. fwrite(&j,sizeof(int),1,beamlet_batch_file);
  117. printf("Error in the calculation for beamlet number %d,\n so all zeros were saved for the resulting beamlet file.\n",b);
  118. printf("%s\n",errstr);
  119. }
  120. else
  121. printf("Successfully calculated beamlet %d of %s.\n",b,argv[3]);
  122. }
  123. // close the beamlet file
  124. fclose(beamlet_batch_file);
  125. // free the density grid
  126. free(density.matrix);
  127. // only need to free angular and radial boundaries for the first
  128. // kernel, since other kernel boundaries just point to the same place
  129. free(mono_kernels.kernel[0].angular_boundary);
  130. free(mono_kernels.kernel[0].radial_boundary);
  131. // free the kernels
  132. free(mono_kernels.energy);
  133. free(mono_kernels.fluence);
  134. free(mono_kernels.mu);
  135. free(mono_kernels.mu_en);
  136. for (i=0;i<mono_kernels.nkernels;i++)
  137. for (j=0;j<N_KERNEL_CATEGORIES;j++)
  138. free(mono_kernels.kernel[i].matrix[j]);
  139. return(SUCCESS);
  140. }
  141. int convolution(MONO_KERNELS *mono_kernels, BEAM *beam, FLOAT_GRID *density, FILE *beamlet_batch_file)
  142. // routine that actually performs the convolution for a given kernel, beam, and geometry
  143. {
  144. int i,j,M,N,Q,Nind,Ntotal; // dimensions of the CT density grid
  145. int *dose_ind;
  146. float *dose_data, doseMax;
  147. char tmpstr[200]; // temporary string
  148. FLOAT_GRID terma_mask;
  149. FLOAT_GRID dose_mask;
  150. FLOAT_GRID deff;
  151. FLOAT_GRID terma;
  152. FLOAT_GRID kermac;
  153. FLOAT_GRID dose;
  154. KERNEL poly_kernel;
  155. // copy the density grid dimensions to the calculation grids
  156. copy_grid_geometry(density,&terma_mask);
  157. copy_grid_geometry(density,&dose_mask);
  158. copy_grid_geometry(density,&deff);
  159. copy_grid_geometry(density,&terma);
  160. copy_grid_geometry(density,&kermac);
  161. copy_grid_geometry(density,&dose);
  162. // dimensions of all the grids
  163. M = density->x_count;
  164. N = density->y_count;
  165. Q = density->z_count;
  166. Ntotal = M*N*Q;
  167. // Allocate memory for all of the grids and fill them all with zeros
  168. if ((terma_mask.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  169. {
  170. sprintf(errstr,"Could not allocate memory for terma_mask.");
  171. return(FAILURE);
  172. }
  173. if ((dose_mask.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  174. {
  175. sprintf(errstr,"Could not allocate memory for dose_mask.");
  176. return(FAILURE);
  177. }
  178. if ((deff.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  179. {
  180. sprintf(errstr,"Could not allocate memory for deff.");
  181. return(FAILURE);
  182. }
  183. if ((terma.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  184. {
  185. sprintf(errstr,"Could not allocate memory for terma.");
  186. return(FAILURE);
  187. }
  188. if ((kermac.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  189. {
  190. sprintf(errstr,"Could not allocate memory for kermac.");
  191. return(FAILURE);
  192. }
  193. if ((dose.matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  194. {
  195. sprintf(errstr,"Could not allocate memory for dose.");
  196. return(FAILURE);
  197. }
  198. for (i=0;i<Ntotal;i++)
  199. {
  200. terma_mask.matrix[i] = 0.0;
  201. dose_mask.matrix[i] = 0.0;
  202. deff.matrix[i] = 0.0;
  203. terma.matrix[i] = 0.0;
  204. kermac.matrix[i] = 0.0;
  205. dose.matrix[i] = 0.0;
  206. }
  207. /* start calculations */
  208. // If a failure occurs in any calculation, append the error
  209. // onto the error string and then return a FAILURE.
  210. // create terma and dose masks
  211. if (SUCCESS != terma_dose_masks(&terma_mask,&dose_mask,beam))
  212. {
  213. sprintf(tmpstr,"Failed in terma_dose_masks!\n");
  214. strcat(tmpstr,errstr);
  215. strcpy(errstr,tmpstr);
  216. return(FAILURE);
  217. }
  218. //create polyenergetic kernel from mono kernels and fluence,mu data
  219. if (SUCCESS != make_poly_kernel(mono_kernels,&poly_kernel) )
  220. {
  221. sprintf(tmpstr,"Failed making polyenergetic kernel!\n");
  222. strcat(tmpstr,errstr);
  223. strcpy(errstr,tmpstr);
  224. return(FAILURE);
  225. }
  226. //create effective depth array from density array
  227. if (SUCCESS != calc_deff(density,&deff,&terma_mask,beam))
  228. {
  229. sprintf(tmpstr,"Failed in calc_deff!\n");
  230. strcat(tmpstr,errstr);
  231. strcpy(errstr,tmpstr);
  232. return(FAILURE);
  233. }
  234. //create kerma and terma arrays
  235. //note kerma is collision kerma and is used for a kernel hardening correction
  236. if (SUCCESS != terma_kerma(&deff,&terma,&kermac,mono_kernels,&terma_mask))
  237. {
  238. sprintf(tmpstr,"Failed in terma_kerma calculation!\n");
  239. strcat(tmpstr,errstr);
  240. strcpy(errstr,tmpstr);
  241. return(FAILURE);
  242. }
  243. //use all this stuff to calculate dose
  244. if ( (SUCCESS != calc_dose(density,&terma,&dose,&poly_kernel,beam,&dose_mask)) )
  245. {
  246. sprintf(tmpstr,"Failed calculating dose!\n");
  247. strcat(tmpstr,errstr);
  248. strcpy(errstr,tmpstr);
  249. return(FAILURE);
  250. }
  251. /* //diagnostic lines:
  252. FILE *fid;
  253. fid = fopen("dose.bin","wb");
  254. fwrite(dose.matrix,sizeof(float),Ntotal,fid);
  255. fclose(fid);
  256. fid = fopen("terma.bin","wb");
  257. fwrite(terma.matrix,sizeof(float),Ntotal,fid);
  258. fclose(fid);
  259. fid = fopen("kermac.bin","wb");
  260. fwrite(kermac.matrix,sizeof(float),Ntotal,fid);
  261. fclose(fid);
  262. fid = fopen("deff.bin","wb");
  263. fwrite(deff.matrix,sizeof(float),Ntotal,fid);
  264. fclose(fid);
  265. fid = fopen("terma_mask.bin","wb");
  266. fwrite(terma_mask.matrix,sizeof(float),Ntotal,fid);
  267. fclose(fid);
  268. fid = fopen("dose_mask.bin","wb");
  269. fwrite(dose_mask.matrix,sizeof(float),Ntotal,fid);
  270. fclose(fid);
  271. fid = fopen("density.bin","wb");
  272. fwrite(density->matrix,sizeof(float),Ntotal,fid);
  273. fclose(fid); */
  274. // find maximum dose
  275. doseMax = 0.0;
  276. for (i=0;i<Ntotal;i++)
  277. if (dose.matrix[i] > doseMax)
  278. doseMax = dose.matrix[i];
  279. // count the number of non-zero dose values
  280. Nind = 0;
  281. for (i=0;i<Ntotal;i++)
  282. if (dose.matrix[i] >= doseMax*doseCutoffThreshold
  283. && dose.matrix[i] > 0.0)
  284. Nind++;
  285. else
  286. dose.matrix[i] = 0.0; // turn off doses below threshold
  287. // allocate memory for sparse dose data
  288. dose_ind = (int *)malloc(sizeof(int)*Nind);
  289. dose_data = (float *)malloc(sizeof(float)*Nind);
  290. // store the sparse data
  291. j = 0; // index just for the sparse data
  292. for (i=0;i<Ntotal;i++)
  293. if (dose.matrix[i] >= doseMax*doseCutoffThreshold
  294. && dose.matrix[i] > 0.0)
  295. {
  296. dose_ind[j] = i;
  297. dose_data[j] = dose.matrix[i];
  298. j++;
  299. }
  300. // save dose to a file
  301. // save the total file size first, then the number of non-zero elements
  302. fwrite(&beam->num,sizeof(int),1,beamlet_batch_file);
  303. fwrite(&M,sizeof(int),1,beamlet_batch_file);
  304. fwrite(&N,sizeof(int),1,beamlet_batch_file);
  305. fwrite(&Q,sizeof(int),1,beamlet_batch_file);
  306. fwrite(&Nind,sizeof(int),1,beamlet_batch_file);
  307. fwrite(dose_ind,sizeof(int),Nind,beamlet_batch_file);
  308. fwrite(dose_data,sizeof(float),Nind,beamlet_batch_file);
  309. free(dose_ind);
  310. free(dose_data);
  311. // free the calculation grids
  312. free(terma_mask.matrix);
  313. free(dose_mask.matrix);
  314. free(deff.matrix);
  315. free(terma.matrix);
  316. free(kermac.matrix);
  317. free(dose.matrix);
  318. free(poly_kernel.angular_boundary);
  319. free(poly_kernel.radial_boundary);
  320. // free(poly_kernel.total_matrix);
  321. for (j=0;j<N_KERNEL_CATEGORIES;j++)
  322. free(poly_kernel.matrix[j]);
  323. return(SUCCESS);
  324. }