parse_funcCondor.cpp 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825
  1. /* parse_func.cpp */
  2. /* Checks the validity of the arguments input to the convolution routine and
  3. transfers the Matlab-style arguments to the C structures defined in defs.h. */
  4. /* This parse function operates on file inputs rather than on Matlab inputs. */
  5. #include "defs.h"
  6. // Markers that are searched for inside the kernel_filenames file, which is
  7. // passed to the load_kernels routine. The line after each of these markers
  8. // in the kernel_filenames file is a filename corresponding to the marker.
  9. #define kernel_header_line "kernel_header"
  10. #define kernel_radii_line "kernel_radii"
  11. #define kernel_angles_line "kernel_angles"
  12. #define kernel_energies_line "kernel_energies"
  13. #define kernel_primary_line "kernel_primary"
  14. #define kernel_first_scatter_line "kernel_first_scatter"
  15. #define kernel_second_scatter_line "kernel_second_scatter"
  16. #define kernel_multiple_scatter_line "kernel_multiple_scatter"
  17. #define kernel_brem_annih_line "kernel_brem_annih"
  18. #define kernel_total_line "kernel_total"
  19. #define kernel_fluence_line "kernel_fluence"
  20. #define kernel_mu_line "kernel_mu"
  21. #define kernel_mu_en_line "kernel_mu_en"
  22. // Markers that are searched for inside the geometry_filenames, which is
  23. // passed to the load_geometry routine. These have the same meaning as
  24. // the kernel_filenames markers.
  25. #define geometry_header "./geometry_files/geometry_header.txt"
  26. #define geometry_density "./geometry_files/density.bin"
  27. #define beam_data "./geometry_files/beam_data.txt"
  28. #define geometry_header_line "geometry_header"
  29. #define geometry_density_line "geometry_density"
  30. #define beam_data_line "beam_data"
  31. extern char errstr[200]; // error string that all routines have access to
  32. int load_kernels(MONO_KERNELS *mono_kernels, char kernel_filenames[])
  33. /* Ensures that the kernel files have the following format:
  34. radii: [1xNradii float]
  35. angles: [1xNangles float]
  36. energies: [1xNenergies float]
  37. primary: [Nradii x Nangles x Nenergies float]
  38. first_scatter: [Nradii x Nangles x Nenergies float]
  39. second_scatter: [Nradii x Nangles x Nenergies float]
  40. multiple_scatter: [Nradii x Nangles x Nenergies float]
  41. brem_annih: [Nradii x Nangles x Nenergies float]
  42. total: [Nradii x Nangles x Nenergies float]
  43. mean_radii: [Nradii x Nangles x Nenergies float] (not used)
  44. mean_angles: [Nradii x Nangles x Nenergies float] (not used)
  45. helpfield: [any x any char]
  46. fluence: [1xNenergies float]
  47. mu: [1xNenergies float]
  48. mu_en: [1xNenergies float]
  49. mistakes or inconsistencies will result in errors.
  50. Results are then stored in mono_kernels.
  51. The names of the files containing all of these parameters are given
  52. by the kernel_filenames file.
  53. */
  54. {
  55. int j,k;
  56. int Nenergies, Nangles, Nradii;
  57. char str[200];
  58. // some strings to hold filenames
  59. char header_filename[200], radii_filename[200], angles_filename[200];
  60. char energies_filename[200], primary_filename[200], first_scatter_filename[200];
  61. char second_scatter_filename[200], multiple_scatter_filename[200];
  62. char brem_annih_filename[200], total_filename[200], fluence_filename[200];
  63. char mu_filename[200], mu_en_filename[200];
  64. // flags for file readin
  65. int header_flag = 0;
  66. int radii_flag = 0;
  67. int angles_flag = 0;
  68. int energies_flag = 0;
  69. int primary_flag = 0;
  70. int first_flag = 0;
  71. int second_flag = 0;
  72. int multiple_flag = 0;
  73. int brem_annih_flag = 0;
  74. int total_flag = 0;
  75. int fluence_flag = 0;
  76. int mu_flag = 0;
  77. int mu_en_flag = 0;
  78. FILE *fid; // generic filename
  79. FILE *primary, *first, *second, *multiple, *brem_annih, *total; // kernel files
  80. // read in the data filenames
  81. if ( (fid = fopen(kernel_filenames,"r")) == NULL) // open the kernel filenames
  82. {
  83. sprintf(errstr,"Could not open the kernel filenames file");
  84. return(FAILURE);
  85. }
  86. while (fscanf(fid,"%s\n",str) != EOF) // read until the end of the file
  87. {
  88. if (!strcmp(str,kernel_header_line))
  89. if (fscanf(fid,"%s\n",header_filename) != 1)
  90. {
  91. sprintf(errstr,"Failed to read-in the kernel_header filename from the kernel filenames file.");
  92. return(FAILURE);
  93. }
  94. else
  95. header_flag = 1;
  96. else if (!strcmp(str,kernel_radii_line))
  97. if (fscanf(fid,"%s\n",radii_filename) != 1)
  98. {
  99. sprintf(errstr,"Failed to read-in the kernel_radii filename from the kernel filenames file.");
  100. return(FAILURE);
  101. }
  102. else
  103. radii_flag = 1;
  104. else if (!strcmp(str,kernel_angles_line))
  105. if (fscanf(fid,"%s\n",angles_filename) != 1)
  106. {
  107. sprintf(errstr,"Failed to read-in the kernel_angles filename from the kernel filenames file.");
  108. return(FAILURE);
  109. }
  110. else
  111. angles_flag = 1;
  112. else if (!strcmp(str,kernel_energies_line))
  113. if (fscanf(fid,"%s\n",energies_filename) != 1)
  114. {
  115. sprintf(errstr,"Failed to read-in the kernel_energies filename from the kernel filenames file.");
  116. return(FAILURE);
  117. }
  118. else
  119. energies_flag = 1;
  120. else if (!strcmp(str,kernel_fluence_line))
  121. if (fscanf(fid,"%s\n",fluence_filename) != 1)
  122. {
  123. sprintf(errstr,"Failed to read-in the kernel_fluence filename from the kernel filenames file.");
  124. return(FAILURE);
  125. }
  126. else
  127. fluence_flag = 1;
  128. else if (!strcmp(str,kernel_primary_line))
  129. if (fscanf(fid,"%s\n",primary_filename) != 1)
  130. {
  131. sprintf(errstr,"Failed to read-in the kernel_primary filename from the kernel filenames file.");
  132. return(FAILURE);
  133. }
  134. else
  135. primary_flag = 1;
  136. else if (!strcmp(str,kernel_first_scatter_line))
  137. if (fscanf(fid,"%s\n",first_scatter_filename) != 1)
  138. {
  139. sprintf(errstr,"Failed to read-in the kernel_first_scatter filename from the kernel filenames file.");
  140. return(FAILURE);
  141. }
  142. else
  143. first_flag = 1;
  144. else if (!strcmp(str,kernel_second_scatter_line))
  145. if (fscanf(fid,"%s\n",second_scatter_filename) != 1)
  146. {
  147. sprintf(errstr,"Failed to read-in the kernel_second_scatter filename from the kernel filenames file.");
  148. return(FAILURE);
  149. }
  150. else
  151. second_flag = 1;
  152. else if (!strcmp(str,kernel_multiple_scatter_line))
  153. if (fscanf(fid,"%s\n",multiple_scatter_filename) != 1)
  154. {
  155. sprintf(errstr,"Failed to read-in the kernel_multiple_scatter filename from the kernel filenames file.");
  156. return(FAILURE);
  157. }
  158. else
  159. multiple_flag = 1;
  160. else if (!strcmp(str,kernel_brem_annih_line))
  161. if (fscanf(fid,"%s\n",brem_annih_filename) != 1)
  162. {
  163. sprintf(errstr,"Failed to read-in the kernel_brem_annih filename from the kernel filenames file.");
  164. return(FAILURE);
  165. }
  166. else
  167. brem_annih_flag = 1;
  168. else if (!strcmp(str,kernel_total_line))
  169. if (fscanf(fid,"%s\n",total_filename) != 1)
  170. {
  171. sprintf(errstr,"Failed to read-in the kernel_total filename from the kernel filenames file.");
  172. return(FAILURE);
  173. }
  174. else
  175. total_flag = 1;
  176. else if (!strcmp(str,kernel_mu_line))
  177. if (fscanf(fid,"%s\n",mu_filename) != 1)
  178. {
  179. sprintf(errstr,"Failed to read-in the kernel_mu filename from the kernel filenames file.");
  180. return(FAILURE);
  181. }
  182. else
  183. mu_flag = 1;
  184. else if (!strcmp(str,kernel_mu_en_line))
  185. if (fscanf(fid,"%s\n",mu_en_filename) != 1)
  186. {
  187. sprintf(errstr,"Failed to read-in the kernel_mu_en filename from the kernel filenames file.");
  188. return(FAILURE);
  189. }
  190. else
  191. mu_en_flag = 1;
  192. else
  193. {
  194. sprintf(errstr,"Unrecognized string in the kernel filenames file: %s\n",str);
  195. return(FAILURE);
  196. }
  197. }
  198. // confirm that all of the required filenames have been found
  199. if (header_flag == 0)
  200. {
  201. sprintf(errstr,"Unable to find a kernel header filename.");
  202. return(FAILURE);
  203. }
  204. else if (radii_flag == 0)
  205. {
  206. sprintf(errstr,"Unable to find a kernel radii filename.");
  207. return(FAILURE);
  208. }
  209. else if (angles_flag == 0)
  210. {
  211. sprintf(errstr,"Unable to find a kernel angles filename.");
  212. return(FAILURE);
  213. }
  214. else if (energies_flag == 0)
  215. {
  216. sprintf(errstr,"Unable to find a kernel energies filename.");
  217. return(FAILURE);
  218. }
  219. else if (primary_flag == 0)
  220. {
  221. sprintf(errstr,"Unable to find a primary kernel filename.");
  222. return(FAILURE);
  223. }
  224. else if (first_flag == 0)
  225. {
  226. sprintf(errstr,"Unable to find a first scatter kernel filename.");
  227. return(FAILURE);
  228. }
  229. else if (second_flag == 0)
  230. {
  231. sprintf(errstr,"Unable to find a second scatter kernel filename.");
  232. return(FAILURE);
  233. }
  234. else if (multiple_flag == 0)
  235. {
  236. sprintf(errstr,"Unable to find a multiple scatter kernel filename.");
  237. return(FAILURE);
  238. }
  239. else if (brem_annih_flag == 0)
  240. {
  241. sprintf(errstr,"Unable to find a brem_annih kernel filename.");
  242. return(FAILURE);
  243. }
  244. else if (total_flag == 0)
  245. {
  246. sprintf(errstr,"Unable to find a total kernel filename.");
  247. return(FAILURE);
  248. }
  249. else if (fluence_flag == 0)
  250. {
  251. sprintf(errstr,"Unable to find a kernel fluence filename.");
  252. return(FAILURE);
  253. }
  254. else if (mu_flag == 0)
  255. {
  256. sprintf(errstr,"Unable to find a kernel mu filename.");
  257. return(FAILURE);
  258. }
  259. else if (mu_en_flag == 0)
  260. {
  261. sprintf(errstr,"Unable to find a kernel mu_en filename.");
  262. return(FAILURE);
  263. }
  264. // read in the expected matrix sizes
  265. if ( (fid=fopen(header_filename,"r")) == NULL)
  266. {
  267. sprintf(errstr,"Missing kernel header file in kernel data folder.");
  268. return(FAILURE);
  269. }
  270. else
  271. {
  272. fgets(str,100,fid); // pop-off the first line of the header
  273. if (fscanf(fid,"%d %d %d\n",&Nradii,&Nangles,&Nenergies) != 3)
  274. {
  275. sprintf(errstr,"Incorrect amount of data in kernel header file.");
  276. return(FAILURE);
  277. }
  278. fclose(fid);
  279. }
  280. mono_kernels->nkernels = Nenergies;
  281. // read in the energies, fluences, mu, and mu_en values
  282. if ( (fid = fopen(energies_filename,"rb")) == NULL)
  283. {
  284. sprintf(errstr,"Missing kernel energies file.");
  285. return(FAILURE);
  286. }
  287. else
  288. {
  289. if ((mono_kernels->energy = (float *)malloc(sizeof(float)*Nenergies)) == NULL)
  290. {
  291. sprintf(errstr,"Failed in memory allocation for kernel energies.");
  292. return(FAILURE);
  293. }
  294. if( (int)fread(mono_kernels->energy,sizeof(float),Nenergies,fid) != Nenergies)
  295. {
  296. sprintf(errstr,"Incorrect number of energies in kernel_energies file.");
  297. return(FAILURE);
  298. }
  299. fclose(fid);
  300. }
  301. if ( (fid = fopen(fluence_filename,"rb")) == NULL)
  302. {
  303. sprintf(errstr,"Missing kernel fluences file.");
  304. return(FAILURE);
  305. }
  306. else
  307. {
  308. if ((mono_kernels->fluence = (float *)malloc(sizeof(float)*Nenergies)) == NULL)
  309. {
  310. sprintf(errstr,"Failed in memory allocation for kernel fluences.");
  311. return(FAILURE);
  312. }
  313. if( (int)fread(mono_kernels->fluence,sizeof(float),Nenergies,fid) != Nenergies)
  314. {
  315. sprintf(errstr,"Incorrect number of fluences in kernel_fluences file.");
  316. return(FAILURE);
  317. }
  318. fclose(fid);
  319. }
  320. if ( (fid = fopen(mu_filename,"rb")) == NULL)
  321. {
  322. sprintf(errstr,"Missing kernel mu file.");
  323. return(FAILURE);
  324. }
  325. else
  326. {
  327. if ((mono_kernels->mu = (float *)malloc(sizeof(float)*Nenergies)) == NULL)
  328. {
  329. sprintf(errstr,"Failed in memory allocation for kernel mu.");
  330. return(FAILURE);
  331. }
  332. if( (int)fread(mono_kernels->mu,sizeof(float),Nenergies,fid) != Nenergies)
  333. {
  334. sprintf(errstr,"Incorrect number of mu values in kernel_mu file.");
  335. return(FAILURE);
  336. }
  337. fclose(fid);
  338. }
  339. if ( (fid = fopen(mu_en_filename,"rb")) == NULL)
  340. {
  341. sprintf(errstr,"Missing kernel mu_en file.");
  342. return(FAILURE);
  343. }
  344. else
  345. {
  346. if ((mono_kernels->mu_en = (float *)malloc(sizeof(float)*Nenergies)) == NULL)
  347. {
  348. sprintf(errstr,"Failed in memory allocation for kernel mu_en.");
  349. return(FAILURE);
  350. }
  351. if( (int)fread(mono_kernels->mu_en,sizeof(float),Nenergies,fid) != Nenergies)
  352. {
  353. sprintf(errstr,"Incorrect number of mu_en values in kernel_mu_en file.");
  354. return(FAILURE);
  355. }
  356. fclose(fid);
  357. }
  358. // Only assign memory for bin boundaries for the first kernel. Just point the other
  359. // boundary pointers at the values in the first kernel.
  360. if ( (fid = fopen(radii_filename,"rb")) == NULL)
  361. {
  362. sprintf(errstr,"Missing kernel radii file.");
  363. return(FAILURE);
  364. }
  365. else
  366. {
  367. mono_kernels->kernel[0].nradii = Nradii;
  368. if( (mono_kernels->kernel[0].radial_boundary = (float *)malloc(sizeof(float)*Nradii)) == NULL)
  369. {
  370. sprintf(errstr,"Failed to allocate memory for radial boundaries.");
  371. return(FAILURE);
  372. }
  373. if( (int)fread(mono_kernels->kernel[0].radial_boundary,sizeof(float),Nradii,fid) != Nradii)
  374. {
  375. sprintf(errstr,"Incorrect number of radii values in kernel_radii file.");
  376. return(FAILURE);
  377. }
  378. fclose(fid);
  379. }
  380. if ( (fid = fopen(angles_filename,"rb")) == NULL)
  381. {
  382. sprintf(errstr,"Missing kernel angles file.");
  383. return(FAILURE);
  384. }
  385. else
  386. {
  387. mono_kernels->kernel[0].ntheta = Nangles;
  388. if( (mono_kernels->kernel[0].angular_boundary = (float *)malloc(sizeof(float)*Nangles)) == NULL)
  389. {
  390. sprintf(errstr,"Failed to allocate memory for angular boundaries.");
  391. return(FAILURE);
  392. }
  393. if( (int)fread(mono_kernels->kernel[0].angular_boundary,sizeof(float),Nangles,fid) != Nangles)
  394. {
  395. sprintf(errstr,"Incorrect number of angular values in kernel_angles file.");
  396. return(FAILURE);
  397. }
  398. fclose(fid);
  399. }
  400. // allocate space for kernel matrices and copy bin boundaries for other kernel categories
  401. for (j=0;j<Nenergies;j++) // loop through energies
  402. {
  403. if (j != 0)
  404. {
  405. mono_kernels->kernel[j].nradii = Nradii;
  406. mono_kernels->kernel[j].ntheta = Nangles;
  407. mono_kernels->kernel[j].radial_boundary = mono_kernels->kernel[0].radial_boundary;
  408. mono_kernels->kernel[j].angular_boundary = mono_kernels->kernel[0].angular_boundary;
  409. }
  410. // allocate space for kernels
  411. for (k=0;k<N_KERNEL_CATEGORIES;k++)
  412. if ((mono_kernels->kernel[j].matrix[k]
  413. = (float *)malloc(sizeof(float)*Nradii*Nangles*N_KERNEL_CATEGORIES)) == NULL)
  414. {
  415. sprintf(errstr,"Failed to allocate memory for kernel matrix.");
  416. return(FAILURE);
  417. }
  418. }
  419. // load-in the kernels one at a time and fill-up the monoenergetic kernels
  420. // open the kernel files
  421. if ( (primary = fopen(primary_filename,"rb")) == NULL)
  422. {
  423. sprintf(errstr,"Missing primary kernel.");
  424. return(FAILURE);
  425. }
  426. if ( (first = fopen(first_scatter_filename,"rb")) == NULL)
  427. {
  428. sprintf(errstr,"Missing first scatter kernel.");
  429. return(FAILURE);
  430. }
  431. if ( (second = fopen(second_scatter_filename,"rb")) == NULL)
  432. {
  433. sprintf(errstr,"Missing second scatter kernel.");
  434. return(FAILURE);
  435. }
  436. if ( (multiple = fopen(multiple_scatter_filename,"rb")) == NULL)
  437. {
  438. sprintf(errstr,"Missing multiple scatter kernel.");
  439. return(FAILURE);
  440. }
  441. if ( (brem_annih = fopen(brem_annih_filename,"rb")) == NULL)
  442. {
  443. sprintf(errstr,"Missing brem_annih kernel.");
  444. return(FAILURE);
  445. }
  446. if ( (total = fopen(total_filename,"rb")) == NULL)
  447. {
  448. sprintf(errstr,"Missing total kernel.");
  449. return(FAILURE);
  450. }
  451. // loop through all energies, reading the kernel files in to the kernel structures
  452. for (j=0;j<Nenergies;j++)
  453. {
  454. if ( (int)fread(mono_kernels->kernel[j].matrix[0],sizeof(float),Nradii*Nangles,primary)
  455. != Nradii*Nangles)
  456. {
  457. sprintf(errstr,"Could not read the correct number of kernel values.");
  458. return(FAILURE);
  459. }
  460. if ( (int)fread(mono_kernels->kernel[j].matrix[1],sizeof(float),Nradii*Nangles,first)
  461. != Nradii*Nangles)
  462. {
  463. sprintf(errstr,"Could not read the correct number of kernel values.");
  464. return(FAILURE);
  465. }
  466. if ( (int)fread(mono_kernels->kernel[j].matrix[2],sizeof(float),Nradii*Nangles,second)
  467. != Nradii*Nangles)
  468. {
  469. sprintf(errstr,"Could not read the correct number of kernel values.");
  470. return(FAILURE);
  471. }
  472. if ( (int)fread(mono_kernels->kernel[j].matrix[3],sizeof(float),Nradii*Nangles,multiple)
  473. != Nradii*Nangles)
  474. {
  475. sprintf(errstr,"Could not read the correct number of kernel values.");
  476. return(FAILURE);
  477. }
  478. if ( (int)fread(mono_kernels->kernel[j].matrix[4],sizeof(float),Nradii*Nangles,brem_annih)
  479. != Nradii*Nangles)
  480. {
  481. sprintf(errstr,"Could not read the correct number of kernel values.");
  482. return(FAILURE);
  483. }
  484. }
  485. // close the kernel files
  486. fclose(primary);
  487. fclose(first);
  488. fclose(second);
  489. fclose(multiple);
  490. fclose(brem_annih);
  491. return(SUCCESS);
  492. }
  493. int load_geometry(FLOAT_GRID *density, char geometry_filenames[])
  494. /* Ensure that the GEOMETRY structure has the following format:
  495. Geometry =
  496. start: [3 element float]
  497. voxel_size: [3 element float]
  498. data: [M x N x Q float]
  499. beam: [1x1 struct]
  500. Geometry.beam =
  501. ip: [3 element float]
  502. jp: [3 element float]
  503. kp: [3 element float]
  504. y_vec: [3 element float]
  505. xp: [float scalar]
  506. yp: [float scalar]
  507. del_xp: [float scalar]
  508. del_yp: [float scalar]
  509. SAD: [float scalar]
  510. and load the data into memory.
  511. The names of the files containing the above data are listed in two files:
  512. geometry_filenames and beam_filenames. The beam data are stored in a
  513. separate file since many dose calculations are done by changing the
  514. beam data and leaving the geometry data constant. This way one can use
  515. a separate beam file for each beam but the same geometry file.
  516. */
  517. {
  518. int Ntotal;
  519. char str[200]; // dummy string
  520. char header_filename[200], density_filename[200];
  521. // flags for filename readin
  522. int header_flag = 0;
  523. int density_flag = 0;
  524. FILE *fid; // dummy filename
  525. // read in the geometry filenames
  526. if ( (fid = fopen(geometry_filenames,"r")) == NULL) // open the geometry filenames
  527. {
  528. sprintf(errstr,"Could not open the geometry filenames file\n");
  529. return(FAILURE);
  530. }
  531. // read in the non-beam geometric data
  532. while (fscanf(fid,"%s\n",str) != EOF) // read until the end of the file
  533. {
  534. if (!strcmp(str,geometry_header_line))
  535. if (fscanf(fid,"%s\n",header_filename) != 1)
  536. {
  537. sprintf(errstr,"Failed to read-in the geometry_header filename from the geometry filenames file.");
  538. }
  539. else
  540. header_flag = 1;
  541. else if (!strcmp(str,geometry_density_line))
  542. if (fscanf(fid,"%s\n",density_filename) != 1)
  543. {
  544. sprintf(errstr,"Failed to read-in the geometry_density filename from the geometry filenames file.");
  545. return(FAILURE);
  546. }
  547. else
  548. density_flag = 1;
  549. else
  550. {
  551. sprintf("Unrecognized string in the geometry filenames file: %s\n",str);
  552. return(FAILURE);
  553. }
  554. }
  555. // confirm that all of the required filenames have been found
  556. if (header_flag == 0)
  557. {
  558. sprintf(errstr,"Unable to find a geometry header filename.\n");
  559. return(FAILURE);
  560. }
  561. else if (density_flag == 0)
  562. {
  563. sprintf(errstr,"Unable to find a geometry density filename.\n");
  564. return(FAILURE);
  565. }
  566. // read in geometric and beam data
  567. if( (fid = fopen(header_filename,"r")) == NULL)
  568. {
  569. sprintf(errstr,"Could not open geometry header file.");
  570. return(FAILURE);
  571. }
  572. else
  573. {
  574. // pop-off the first line of the header file
  575. if (fgets(str,100,fid) == NULL)
  576. {
  577. sprintf(errstr,"Could not read from geometry header file.");
  578. return(FAILURE);
  579. }
  580. // read-in the CT data grid size
  581. if (fscanf(fid,"%d %d %d\n",&density->x_count,&density->y_count,&density->z_count) != 3)
  582. {
  583. sprintf(errstr,"Could not read-in data grid size.");
  584. return(FAILURE);
  585. }
  586. // pop-off the next line of the header file
  587. if (fgets(str,100,fid) == NULL)
  588. {
  589. sprintf(errstr,"Could not read from geometry header file.");
  590. return(FAILURE);
  591. }
  592. // read-in the CT data grid size
  593. if (fscanf(fid,"%f %f %f\n",&density->start.x,&density->start.y,&density->start.z) != 3)
  594. {
  595. sprintf(errstr,"Could not read-in density grid start position.");
  596. return(FAILURE);
  597. }
  598. // pop-off the next line of the header file
  599. if (fgets(str,100,fid) == NULL)
  600. {
  601. sprintf(errstr,"Could not read from geometry header file.");
  602. return(FAILURE);
  603. }
  604. // read-in the CT data grid size
  605. if (fscanf(fid,"%f %f %f\n",&density->inc.x,&density->inc.y,&density->inc.z) != 3)
  606. {
  607. sprintf(errstr,"Could not read-in voxel size vector.");
  608. return(FAILURE);
  609. }
  610. Ntotal = density->x_count*density->y_count*density->z_count;
  611. }
  612. // read-in the CT density data
  613. if( (fid = fopen(density_filename,"rb")) == NULL)
  614. {
  615. sprintf(errstr,"Could not open CT density file.");
  616. return(FAILURE);
  617. }
  618. else
  619. {
  620. if( (density->matrix = (float *)malloc(sizeof(float)*Ntotal)) == NULL)
  621. {
  622. sprintf(errstr,"Unable to allocate space for CT density grid.");
  623. return(FAILURE);
  624. }
  625. else if ((int)fread(density->matrix,sizeof(float),Ntotal,fid) != Ntotal)
  626. {
  627. sprintf(errstr,"Unable to read-in CT density data.");
  628. return(FAILURE);
  629. }
  630. }
  631. return(SUCCESS);
  632. }
  633. int pop_beam(BEAM *bm, FILE *beamspec_batch_file)
  634. /* Pops a beam off of the beamspec_batch_file file and loads it into a BEAM structure.
  635. Ensure that the beam data file has the following format:
  636. Geometry.beam =
  637. ip: [3 element float]
  638. jp: [3 element float]
  639. kp: [3 element float]
  640. y_vec: [3 element float]
  641. xp: [float scalar]
  642. yp: [float scalar]
  643. del_xp: [float scalar]
  644. del_yp: [float scalar]
  645. SAD: [float scalar]
  646. num: [int scalar]
  647. and load the data into memory.
  648. */
  649. {
  650. char str[200]; // dummy string
  651. // read-in the beam data
  652. // pop-off the first line of the beam data file
  653. // read-in the CT data grid size
  654. if (fgets(str,100,beamspec_batch_file) == NULL)
  655. {
  656. sprintf(errstr,"Could not read from beam data file.");
  657. return(FAILURE);
  658. }
  659. if (fscanf(beamspec_batch_file,"%d\n",&bm->num) != 1)
  660. {
  661. sprintf(errstr,"Could not read-in beam data.");
  662. return(FAILURE);
  663. }
  664. if (fgets(str,100,beamspec_batch_file) == NULL)
  665. {
  666. sprintf(errstr,"Could not read from beam data file.");
  667. return(FAILURE);
  668. }
  669. // read-in the CT data grid size
  670. if (fscanf(beamspec_batch_file,"%f %f %f %f %f\n",&bm->SAD,&bm->xp,&bm->yp,&bm->del_xp,&bm->del_yp) != 5)
  671. {
  672. sprintf(errstr,"Could not read-in beam data.");
  673. return(FAILURE);
  674. }
  675. // pop-off the next line of the beam data file
  676. if (fgets(str,100,beamspec_batch_file) == NULL)
  677. {
  678. sprintf(errstr,"Could not read from from beam data file.");
  679. return(FAILURE);
  680. }
  681. // read-in the source position vector
  682. if (fscanf(beamspec_batch_file,"%f %f %f\n",&bm->y_vec[0],&bm->y_vec[1],&bm->y_vec[2]) != 3)
  683. {
  684. sprintf(errstr,"Could not read-in source location vector.");
  685. return(FAILURE);
  686. }
  687. // pop-off the next line of the beam data file
  688. if (fgets(str,100,beamspec_batch_file) == NULL)
  689. {
  690. sprintf(errstr,"Could not read from from beam data file.");
  691. return(FAILURE);
  692. }
  693. // read-in the beam's eye view vectors, ip first
  694. if (fscanf(beamspec_batch_file,"%f %f %f\n",&bm->ip[0],&bm->ip[1],&bm->ip[2]) != 3)
  695. {
  696. sprintf(errstr,"Could not read-in ip vector.");
  697. return(FAILURE);
  698. }
  699. // pop-off the next line of the beam data file
  700. if (fgets(str,100,beamspec_batch_file) == NULL)
  701. {
  702. sprintf(errstr,"Could not read from from beam data file.");
  703. return(FAILURE);
  704. }
  705. // read-in the second beam's eye view vector, jp
  706. if (fscanf(beamspec_batch_file,"%f %f %f\n",&bm->jp[0],&bm->jp[1],&bm->jp[2]) != 3)
  707. {
  708. sprintf(errstr,"Could not read-in jp vector.");
  709. return(FAILURE);
  710. }
  711. // pop-off the next line of the beam data file
  712. if (fgets(str,100,beamspec_batch_file) == NULL)
  713. {
  714. sprintf(errstr,"Could not read from from beam data file.");
  715. return(FAILURE);
  716. }
  717. // read-in the third beam's eye view vector, kp
  718. if (fscanf(beamspec_batch_file,"%f %f %f\n",&bm->kp[0],&bm->kp[1],&bm->kp[2]) != 3)
  719. {
  720. sprintf(errstr,"Could not read-in kp vector.");
  721. return(FAILURE);
  722. }
  723. // ensure that ip,jp,kp are orthonormal and form a right-handed coordinate system
  724. if ( (fabs(bm->ip[0]*bm->ip[0] + bm->ip[1]*bm->ip[1] + bm->ip[2]*bm->ip[2] - 1.0) > 1e-6)
  725. || (fabs(bm->jp[0]*bm->jp[0] + bm->jp[1]*bm->jp[1] + bm->jp[2]*bm->jp[2] - 1.0) > 1e-6)
  726. || (fabs(bm->kp[0]*bm->kp[0] + bm->kp[1]*bm->kp[1] + bm->kp[2]*bm->kp[2] - 1.0) > 1e-6)
  727. || (fabs(bm->ip[0]*bm->jp[0] + bm->ip[1]*bm->jp[1] + bm->ip[2]*bm->jp[2]) > 1e-6)
  728. || (fabs(bm->ip[0]*bm->kp[0] + bm->ip[1]*bm->kp[1] + bm->ip[2]*bm->kp[2]) > 1e-6)
  729. || (fabs(bm->jp[0]*bm->kp[0] + bm->jp[1]*bm->kp[1] + bm->jp[2]*bm->kp[2]) > 1e-6))
  730. {
  731. sprintf(errstr,"Beam's eye view unit vectors ip, jp, and kp must be orthonormal.");
  732. return(FAILURE);
  733. }
  734. // check that cross(ip,jp) = kp
  735. if ( (fabs(bm->ip[1]*bm->jp[2] - bm->ip[2]*bm->jp[1] - bm->kp[0]) > 1e-6)
  736. || (fabs(bm->ip[2]*bm->jp[0] - bm->ip[0]*bm->jp[2] - bm->kp[1]) > 1e-6)
  737. || (fabs(bm->ip[0]*bm->jp[1] - bm->ip[1]*bm->jp[0] - bm->kp[2]) > 1e-6))
  738. {
  739. sprintf(errstr,"Must have cross(ip,jp) = kp for beam's eye view");
  740. return(FAILURE);
  741. }
  742. return(SUCCESS);
  743. }