convolutionCondor.cpp 11 KB

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