2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/rdproj.cpp $
5 $LastChangedRevision: 91 $
6 $Date: 2014-07-02 17:34:27 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 RDPROJ GENERATES THE RANDOM ACCESS VERSION OF THE PROJECTION DATA
12 THE MODIFIER ON PROJECTION CARD IS ONE OF THE FOLLOWING -
13 REAL THIS IS THE DEFAULT VALUE. PROJECTION DATA ARE READ FROM
14 FILE11 AND PUT ON PRJFIL
15 PSEUDO CREATES PSEUDO PROJECTION DATA FROM INFORMATION ON FILE11
16 AND INPUT AND PUTS ON PRJFIL. THIS REQUIRES THAT
17 PICTURE TEST BE ALREADY SPECIFIED.
18 OLD THIS ASSUMES THAT PRJFIL ALREADY EXISTS FROM A PREVIOUS
19 SNARK RUN. THE GEOMETRY AND NOISE INFORMATION ARE IN THE
20 FIRST RECORD OF PRJFIL.
22 SEE SNARK MANUAL FOR MORE DETAILS.
25 *Modification date: July 13, 2008
26 *Author: Joanna Klukowska (jk)
27 *Changes: Added functionality for beam hardening correction.
28 * If the user sets the options for correction in the input
29 * file, rdproj() createas nergy temporary prjfil's with to
30 * pass to bh_correction::correction() for beam hardening correction
31 * which is responsible for creation of the final prjfil
32 * that is used in the remainder of the program.
33 * The option of beam hardening correction is only appropriate
34 * when polychromatic data is provided.
43 #include <sys/types.h>
60 #include "bh_correction.h"
68 #include "bh_correction.h"
75 static const INTEGER proj_codes[4] =
76 { CHAR2INT('o', 'l', 'd', ' '), CHAR2INT('r', 'e', 'a', 'l'), CHAR2INT('p',
77 's', 'e', 'u'), CHAR2INT('b', 'h', 'p', 's') //jk 9/19/2008
78 //used only for recursive runs of snark during the
79 //beam hardening correction process
80 //should not be listed as an option in the manual
83 // jk 7/13/2008 added for beam hardening correction
84 static const INTEGER bhc_codes[1] =
85 { CHAR2INT('b', 'e', 'a', 'm') };
87 static BOOLEAN called = FALSE;
107 INTEGER* lbase = NULL;
108 REAL* wbase = NULL; //wei 3/2005
115 // jk 7/13/2008 added for beam hardening correction
117 INTEGER bhcword; //BHC keyword read from the input file
118 SnarkPrjFile * BHC_PrjFile = NULL; //temporary prjfil's
121 bh_correction bh_cor;
123 ///////////////////////////////////////////////////////////
125 // CHECK FOR MULTIPLE CALLS
131 dtor = (REAL) (Consts.pi / 180.0);
133 // CHECK THAT RDPICT HAS BEEN PERFORMED BEFORE THIS
135 if (GeoPar.nelem <= 0)
138 // DETERMINE THE PROJECTION MODIFIER AND THE INPUT FILE
140 // changed last parameter from "4" to "3". lajos, Nov 18, 2004
141 nword = InFile.getwrd(FALSE, &eol, proj_codes, 4);
143 // IF OLD IS OPTED JUST READ HEADER & QUIT
145 if (nword == proj_codes[0]) // OLD
147 // OPEN PROJECTION FILE FOR PROJECTION OLD CASE
149 if (ProjFile.Open((char *) "prjfil") != 0)
151 fprintf(output, "\n **** unable to open prjfil");
152 fprintf(output, "\n **** program aborted\n");
156 ProjFile.ReadHeader(Creacm.name, &NoisePar, &Spctrm, &Anglst);
158 // compute average density - fixes bug 232 of snark05
159 double *projdata = new double[GeoPar.nrays];
162 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
163 for (np = 0; np < GeoPar.prjnum; ++np)
165 ProjFile.ReadProj(np, projdata, GeoPar.nrays);
166 for (nr = 0; nr < GeoPar.nrays; nr++)
175 * MAX0(fabs(sinth), fabs(costh))));
179 totden += projdata[nr] / GeoPar.pinc;
184 totden += projdata[nr];
186 posit(np, nr, &ax, &ay, &mx, &my);
187 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
191 GeoPar.aveden = totden / totlen;
192 fprintf(output, "\n estimate of totlen = %16.6f", totlen);
193 fprintf(output, "\n estimate of totden = %16.6f", totden);
194 fprintf(output, "\n estimate of average density = %10.4f",
201 // READ PREFACE TO PROJECTION DATA
203 if (nword == proj_codes[2] || nword == proj_codes[3]) // PSEUDO
205 NoisePar.resetNoise(); // jklukowska, bug 263, in case the noise was set
206 // previously, we need to reset it to initial values
211 if (nword == proj_codes[1]) // REAL
215 if (File11.Open(FALSE) != 0)
218 "\n **** unable to open file11 for reading");
219 fprintf(output, "\n **** program aborted\n");
222 File11.isOpen = TRUE;
230 "\n **** OLD, REAL or PSEUDO must be specified");
231 fprintf(output, "\n **** program aborted\n");
239 // jk 7/13/2008 added for beam hardening correction
240 //check if bh_correction flag is set
241 bhcword = InFile.getwrd(FALSE, &eol, bhc_codes, 1);
242 if (bhcword == bhc_codes[0])
243 { //beam hardening correction flag set
245 //check if data is polychromatic
246 if (Spctrm.nergy == 1) //bhc cannot be performed
249 "\n **** Beam hardening correction is not appropriate\n"
250 "for monochromatic data (nergy = 1).\n"
251 "Program will continue ignoring BHC command.\n\n");
257 //read in the parameters for beam hardening correction
258 bh_cor.readInputFilePoly();
260 //allocate temporary prjfil's
261 BHC_PrjFile = new SnarkPrjFile[1];
264 time(&t); //get the current time
265 //create temporary directory
266 sprintf(bh_cor.tempDir, "bhc%i", (int) t);
267 sprintf(command, "mkdir %s", bh_cor.tempDir);
268 if (system(command) != 0)
270 fprintf(stdout, "\nERROR creating temporary directory %s\n",
277 if (ProjFile.Open((char*) "prjfil", GeoPar.usrays, GeoPar.prjnum) != 0)
279 fprintf(output, "\n **** unable to open prjfil");
280 fprintf(output, "\n **** program aborted\n");
284 // jk 7/13/2008 added for beam hardening correction
285 //open temporary prjfil's
288 for (i = 0; i < 1; i++)
290 sprintf(filename, "%s/q_p0", bh_cor.tempDir);
291 if (BHC_PrjFile[i].Open(filename, GeoPar.usrays, GeoPar.prjnum)
294 fprintf(output, "\n **** unable to open temporary prjfil");
295 fprintf(output, "\n **** program aborted\n");
301 ProjFile.WriteHeader(Creacm.name, &GeoPar, &NoisePar, &Spctrm);
303 ProjFile.WriteAngles(Anglst.bth);
305 // jk 7/13/2008 added for beam hardening correction
306 // write header info to temporary prjfil's
309 for (i = 0; i < 1; i++)
311 BHC_PrjFile[i].WriteHeader(Creacm.name, &GeoPar, &NoisePar,
313 BHC_PrjFile[i].WriteAngles(Anglst.bth);
317 // ALLOCATE SPACE FOR ARRAYS LIST,WEIGHT,SPHI,CPHI,PBASE
319 bldlst(&lbase, &wbase);
326 Anglst.pbase = new REAL[GeoPar.nrays];
328 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
330 Anglst.pbase[nr] = 0.0;
333 // IF PROJECTION DATA ARE NOISY WITH BIAS DUE TO SPECIFICATION OF
334 // SKEWED DISTRIBUTIONS PREPARE HERE TO REMOVE BIAS
338 if (NoisePar.ultnmn != 0.0)
340 umean = (REAL) 1.0 / NoisePar.ultnmn;
347 amean = NoisePar.addnmn;
350 // INITIALIZE FOR COMPUTATION OF AVERAGE DENSITY
354 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
356 // CHECK TO SEE IF PHANTOM IS PRESENT WHEN PSEUDO IS OPTED
358 if ((nword == proj_codes[2] || nword == proj_codes[3])
359 && (Creacm.test == NULL))
364 // SELECT PROJECTION BY PROJECTION AND PREPARE PRJFIL
366 for (np = 0; np < GeoPar.prjnum; np++)
368 if (nword == proj_codes[1])
371 if (!File11.ReadProjection(&theta, Anglst.pbase + GeoPar.fusray,
372 GeoPar.lusray - GeoPar.fusray + 1))
375 "\n **** end of file encountered on file11");
376 fprintf(output, "\n **** program aborted\n");
380 delete[] wbase; //wei 3/2005
384 sinth = (REAL) sin(theta);
385 costh = (REAL) cos(theta);
388 the = *(Anglst.bth + np);
391 if ((REAL) fabs(theta - the) > (REAL) 0.001)
393 fprintf(output, "\n **** angles specified for projection "
394 "%4i do not match", np);
395 fprintf(output, "\n angle given in header is %10.4f "
396 "degrees = %10.4f radians", angle, the);
397 fprintf(output, "\n angle given in data is %10.4f "
398 "degrees = %10.4f radians", ang, theta);
399 fprintf(output, "\n **** program aborted\n");
403 delete[] wbase; //wei 3/2005
409 // PROJECTION PSEUDO or BHPSEUDO OPTION
411 theta = Creacm.pang[np];
413 sinth = (REAL) sin(theta);
414 costh = (REAL) cos(theta);
416 // bug221 start - swr - 03/03/07
417 for (nr = 0; nr < GeoPar.nrays; nr++)
418 Anglst.pbase[nr] = 0.0;
420 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
423 if (nword == proj_codes[3])
426 //jk 07/20/2008 modyfying how pseudo projections are
427 // computed - they need to be polychromatic for
428 //beam hardening correction
430 for (j = 0; j < Spctrm.nergy; j++)
432 temp += Spctrm.engwt[j]
434 pseudo(phantom_poly_array[j], np,
435 nr, lbase, wbase, &numb,
443 temp = NoisePar.raysum(
444 pseudo(Creacm.test, np, nr, lbase, wbase, &numb,
445 &snorm, GeoPar.line, FALSE));
447 if (NoisePar.quanin > 0 && NoisePar.quanin != 4)
448 temp *= NoisePar.quanmn;
449 Anglst.pbase[nr] = NoisePar.applyQuantumNoise(np, nr, temp);
452 // TAKE CARE OF NOISE
454 NoisePar.scatterNoise(Anglst.pbase);
456 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
458 REAL temp = NoisePar.computeAttenuation(np, nr,
460 temp = NoisePar.multiplicativeNoise(temp);
461 Anglst.pbase[nr] = NoisePar.additiveNoise(temp);
463 // bug221 end - swr - 03/03/07
465 // REMOVE BIAS IN PROJECTION DATA POSSIBLY INTRODUCED BY A SKEWED
466 // NOISE DISTRIBUTION
468 if (NoisePar.ultnfl || NoisePar.addnfl)
470 for (i = GeoPar.fusray; i <= GeoPar.lusray; i++)
472 Anglst.pbase[i] = (Anglst.pbase[i] - amean) * umean;
476 // ADD UP FOR ESTIMATION OF AVERAGE DENSITY FROM RAYSUMS
478 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
480 // Changed here for SNARK89 Version 1.01 7/13/90 to conform to manual
486 (REAL) (Anglst.pbase[nr]
488 * MAX0(fabs(sinth), fabs(costh))));
492 totden += Anglst.pbase[nr] / GeoPar.pinc;
497 totden += Anglst.pbase[nr];
500 posit(np, nr, &ax, &ay, &mx, &my);
502 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
507 // jk 7/13/2008 added for beam hardening correction
508 // write data to temporary prjfil's if bhc is set
511 for (i = 0; i < 1; i++)
514 //iterate over values of Anglst.pbase and multiply by
515 //coefficients of the correction polynomial
516 //for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
517 INTEGER NoOfRays = GeoPar.lusray - GeoPar.fusray + 1;
519 REAL * temp = new REAL[NoOfRays];
520 for (k = 0; k < NoOfRays; k++)
522 for (nr = 0; nr < NoOfRays; nr++)
524 t = Anglst.pbase[nr + GeoPar.fusray];
526 for (j = 0; j <= bh_cor.noDegree; j++)
528 temp[nr] += bh_cor.a[j] * tt;
533 BHC_PrjFile[i].WriteProj(np, temp); //&(Anglst.pbase[GeoPar.fusray]));
537 else //write data to the final prjfil
541 ProjFile.WriteProj(np, &(Anglst.pbase[GeoPar.fusray]));
546 GeoPar.aveden = totden / totlen;
547 fprintf(output, "\n\n estimate of totlen = %16.6f", totlen);
549 fprintf(output, "\n estimate of totden = %16.6f", totden);
551 fprintf(output, "\n estimate of average density = %10.4f",
554 // WRITE UPDATED HEADER ON PRJFIL
556 // jk 7/13/2008 added for beam hardening correction
557 // write header to temporary prjfil's if bhc is set
560 for (i = 0; i < 1; i++)
562 BHC_PrjFile[i].WriteHeader(Creacm.name, &GeoPar, &NoisePar,
566 else //write header to the final prjfil
569 ProjFile.WriteHeader(Creacm.name, &GeoPar, &NoisePar, &Spctrm);
572 // PREPARE TO EXIT FROM ROUTINE
574 //free(Anglst.bth); deleted now in the destructor
579 // jk 7/13/2008 added for beam hardening correction
582 for (i = 0; i < 1; i++)
584 BHC_PrjFile[i].Close();
593 // jk 7/13/2008 added for beam hardening correction
594 // CALL BEAM HARDENING CORRECTION ROUTINES
598 bh_cor.correction(BHC_PrjFile);
600 // OPEN PROJECTION FILE JUST CORRECTED FOR BEAM HARDENING
602 if (ProjFile.Open((char*) "prjfil") != 0)
605 "\n **** unable to open prjfil after beam hardening");
606 fprintf(output, "\n **** program aborted\n");
610 ProjFile.ReadHeader(Creacm.name, &NoisePar, &Spctrm, &Anglst);
612 // compute average density
613 double *projdata = new double[GeoPar.nrays];
616 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
617 for (np = 0; np < GeoPar.prjnum; ++np)
619 ProjFile.ReadProj(np, projdata, GeoPar.nrays);
620 for (nr = 0; nr < GeoPar.nrays; nr++)
629 * MAX0(fabs(sinth), fabs(costh))));
633 totden += projdata[nr] / GeoPar.pinc;
638 totden += projdata[nr];
640 posit(np, nr, &ax, &ay, &mx, &my);
641 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
645 GeoPar.aveden = totden / totlen;
646 fprintf(output, "\n estimate of totlen = %16.6f", totlen);
647 fprintf(output, "\n estimate of totden = %16.6f", totden);
648 fprintf(output, "\n estimate of average density = %10.4f",
656 fprintf(output, "\n projection data read");
657 fprintf(output, "\n %s", Creacm.name);
659 // WRITE PROJECTION HEADER INTO RECORD nelem + 3 OF
660 // RECONSTRUCTION FILE
662 RecFile.SetProjName(Creacm.name);
664 //delete the created directory bhc
667 sprintf(command, "rm -rf %s", bh_cor.tempDir);
668 if (system(command) != 0)
670 fprintf(stdout, "\n\n****ERROR deleting temporary directory %s\n\n",
678 delete[] wbase; //wei 3/2005
679 if (BHC_PrjFile != NULL)
680 delete[] BHC_PrjFile;
682 for (i = 6; i >= 0; i--)
684 if (phantom_poly_array[i] != NULL)
685 delete[] phantom_poly_array[i];