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/inputfile.cpp $
5 $LastChangedRevision: 86 $
6 $Date: 2014-07-02 16:45:09 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 THE FIRST 7 ROWS OF 'OBJECT' STORE THE DENSITY FOR THE 7 ENERGY
11 LEVELS, 8TH STORES THE AVERAGE DENSITY;
12 9TH ROW SPECIFIES TYPE OF OBJECT, 1 FOR ELLIPSE/CIRCLE
13 2 FOR RECTANGLE/SQUARE
15 4 FOR SEGMENT OF CIRCLE
16 5 FOR SECTOR OF CIRCLE
17 10TH ROW CX X-COORDINATE OF CENTER
18 11TH ROW CY Y-COORDINATE OF CENTER
19 12TH ROW U SEMI X-AXIS LENGTH
20 13TH ROW V SEMI Y-AXIS LENGTH
21 14TH ROW ANG ANGLE OF ROTATION FROM X-AXIS
22 15TH ROW SINANG SIN OF ANGLE OF ROTATION
23 16TH ROW COSANG COSINE OF ANGLE OF ROTATION
24 17TH ROW RAYL LENGTH OF RAY INTERSECTED BY OBJECT
25 THIS VALUE OF SUBSCRIPT REFERRED TO IN CREAPR ONLY
47 #include "inputfile.h"
49 //jk 11/18/2007 introducing variability
52 //long newseed; //the variable Creacm.varseed is used instead
53 //the array stores values of variability added to each pixel of the phantom
54 //for each energy level
55 REAL * variab_array[7];
56 //the flag is set if the value of variability in the input file is other than zero
58 //saves the value of the origical scale which is the mean value for variability calculations
61 REAL * phantom_poly_array[7];
63 void InputFile_class::Open() //open input file (stdin) and
68 fprintf(output, "\n **** unable to open input file");
69 fprintf(output, "\n **** program aborted\n");
74 void InputFile_class::GetName() // get name
77 signed CHAR ch; // changed "CHAR" to "signed CHAR". lajos, Nov 18, 2004
79 fprintf(output, "\n");
81 //assign the name of reconstruction to the 'name' member of
82 //Creacm (defined in global scope in creacm.cpp
83 for (i = 0; i < 80; i++)
84 { // changed "i < 81" to "i < 80". Lajos, Nov 18, 2004
85 //if there are fewer than 80 characters, stop earlier
86 if (((ch = fgetc(infile)) == '\n') || (ch == EOF))
92 //skip over remaining characters on the line (only first 80 characters
93 //are used for the name)
94 while (ch != '\n' && ch != EOF)
95 ch = fgetc(infile); // bug117 - swr - 6/29/05
96 // insufficient records in input; EOF encountered
99 creaer(25); //generate an error
103 Creacm.name[i] = 0; //set the last character to null (it is
104 //a null terminated string)
106 fprintf(output, "\n %s\n", Creacm.name);
109 void InputFile_class::ReadEnergySpectralDistribution() //READ ENERGY
110 //SPECTRAL DISTRIBUTION
114 INTEGER ndred, tword;
117 static const INTEGER spect_codes[3] =
118 { CHAR2INT('s', 'p', 'e', 'c'), CHAR2INT('m', 'o', 'n', 'o'), CHAR2INT('p',
121 tword = getwrd(TRUE, &eol, spect_codes, 1); /*getwrd returns as its value
122 the first four letters of the next string of alpha characters
123 in the input record calls function 'getnxt'.
124 returns eol = TRUE when end of line is reached.*/
127 creaer(spect_codes[0]); //generate error if there are
128 //fewer than 4 characters on the line
130 Creacm.word2 = getwrd(FALSE, &eol, &(spect_codes[1]), 2); //gets the code
131 //for either 'mono' or 'poly' from the input file
133 if (Creacm.word2 == spect_codes[2]) //if 'poly'
135 // PROCESS POLYCHROMATIC DATA HERE; FOR MONOCHROMATIC GO ELESEWHERE
137 Spctrm.nergy = getint(FALSE, &eol);
138 if ((Spctrm.nergy <= 0) || (Spctrm.nergy > 7))
143 for (i = 0; i < Spctrm.nergy; i++)
145 Spctrm.energy[i] = getint(FALSE, &eol);
148 // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
153 Creacm.percnt[i] = getint(FALSE, &eol);
156 // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
163 "\n energy spectrum is polychromatic with %3i energy levels",
165 fprintf(output, "\n with the following photon distribution");
166 fprintf(output, "\n energy level percent");
167 for (i = 0; i < Spctrm.nergy; i++)
169 fprintf(output, "\n%17i %11i", Spctrm.energy[i], Creacm.percnt[i]);
172 // COMPUTE THE WEIGHTS FOR EACH ENERGY LEVEL FOR AVERAGING THE
173 // ABSORPTION OVER THE ENERGY LEVELS
176 for (i = 0; i < Spctrm.nergy; i++)
178 Spctrm.engwt[i] = Creacm.percnt[i] * (REAL) 0.01;
179 ndred += Creacm.percnt[i];
180 if (Creacm.percnt[i] < 0)
189 if (Creacm.word2 == spect_codes[1])
191 // MONOCHROMATIC DATA
194 Creacm.percnt[0] = 100;
195 Spctrm.engwt[0] = 1.0;
196 Spctrm.energy[0] = getint(FALSE, &eol);
199 "\n energy spectrum is monochromatic at energy level %5i\n",
207 void InputFile_class::GetInformationABoutCreationOfPhantom(INTEGER invoke)
208 // GET INFORMATION ABOUT CREATION OF PHANTOM
211 static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
214 static const INTEGER phan_codes[2] =
215 { CHAR2INT('p', 'h', 'a', 'n'), CHAR2INT('a', 'v', 'e', 'r') };
217 tword = getwrd(TRUE, &eol, phan_codes, 1);
219 creaer(phan_codes[0]);
221 fprintf(output, "\n");
223 Creacm.word2 = getwrd(FALSE, &eol, &(phan_codes[1]), 1);
224 Creacm.picflg = (Creacm.word2 == phan_codes[1]);
226 // IF PHANTOM NOT CREATED DONT READ NAVE1
232 fprintf(output, "\n this run will generate a phantom");
234 Creacm.nave1 = getint(FALSE, &eol);
236 "\n density in each pixel is obtained as the average of %i x %i points\n",
237 Creacm.nave1, Creacm.nave1);
239 if ((Creacm.nave1 <= 0) || (Creacm.nave1 == Creacm.nave1 / 2 * 2))
242 // READ IN GRID SIZE AND PIXEL SIZE .
244 GeoPar.nelem = getint(TRUE, &eol);
245 GeoPar.pixsiz = getnum(FALSE, &eol);
246 fprintf(output, "\n picture size %i x %i, pixel size %10.4f\n",
247 GeoPar.nelem, GeoPar.nelem, GeoPar.pixsiz);
249 if ((GeoPar.nelem <= 0) || (GeoPar.nelem == GeoPar.nelem / 2 * 2))
251 if (GeoPar.pixsiz < Consts.zero)
253 GeoPar.midpix = (GeoPar.nelem + 1) / 2;
254 GeoPar.area = GeoPar.nelem * GeoPar.nelem;
259 // READ BACKGROUND ABSORPTION
260 void InputFile_class::ReadBackgroundAbsorption()
264 static const INTEGER hback = CHAR2INT('b', 'a', 'c', 'k');
268 for (i = 0; i < Spctrm.nergy; i++)
270 Spctrm.backgr[i] = getnum(FALSE, &eol);
273 fprintf(output, "\n at levels");
274 fprintf(output, "\n ");
275 for (i = 0; i < Spctrm.nergy; i++)
277 fprintf(output, "%9i ", Spctrm.energy[i]);
280 fprintf(output, "\n background absorption");
281 for (i = 0; i < Spctrm.nergy; i++)
283 fprintf(output, "%9.4f ", Spctrm.backgr[i]);
287 // READ RAYSUM INFORMATION
288 void InputFile_class::ReadRaysumInformation()
291 // static const INTEGER haver = CHAR2INT('a','v','e','r'); // commented out since it is unused. Lajos, Jan 13, 2005
294 static const INTEGER rays_codes[2] =
295 { CHAR2INT('r', 'a', 'y', 's'), CHAR2INT('a', 'v', 'e', 'r') };
297 tword = getwrd(TRUE, &eol, rays_codes, 1);
299 creaer(rays_codes[0]);
301 fprintf(output, "\n");
303 Creacm.word2 = getwrd(FALSE, &eol, &(rays_codes[1]), 1);
304 Creacm.prjflg = (Creacm.word2 == rays_codes[1]);
307 // INPUT ENTIRE GEOMETRY DATA
308 void InputFile_class::ReadGeometryData(INTEGER invoke)
313 static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
315 static const INTEGER geom_codes[10] = /// bug #135. Lajos, 08/01/2005
317 CHAR2INT('g', 'e', 'o', 'm'), CHAR2INT('p', 'a', 'r', 'a'),
318 CHAR2INT('d', 'i', 'v', 'e'), CHAR2INT('l', 'i', 'n', 'o'), /// bug #135. Lajos, 08/01/2005
319 CHAR2INT('u', 'n', 'i', 'f'), CHAR2INT('v', 'a', 'r', 'i'),
320 CHAR2INT('l', 'i', 'n', 'e'), CHAR2INT('s', 't', 'r', 'i'),
321 CHAR2INT('a', 'r', 'c', ' '), CHAR2INT('t', 'a', 'n', 'g') };
323 fprintf(output, "\n");
325 tword = getwrd(TRUE, &eol, geom_codes, 1);
327 creaer(geom_codes[0]);
329 fprintf(output, "\n");
331 // GET RAY TYPE, PARALLEL, DIVERGENT, OR LINOGRAM, UNIF, VARI, TANG OR ARC /// bug #135. Lajos, 08/01/2005
332 tword = getwrd(TRUE, &eol, &(geom_codes[1]), 3); /// bug #135. Lajos, 08/01/2005
336 "\n **** keyword parallel, divergent, or linogram is missing"); /// bug #135. Lajos, 08/01/2005
337 fprintf(output, "\n **** program aborted\n");
340 Creacm.word1 = tword;
342 GeoPar.lino = FALSE; /// bug #135. Lajos, 08/02/2005
344 if (tword == geom_codes[1])
346 tword = getwrd(FALSE, &eol, &(geom_codes[4]), 2); /// bug #135. Lajos, 08/01/2005
349 fprintf(output, "\n **** keyword uniform or variable is missing");
350 fprintf(output, "\n **** program aborted\n");
354 else if (tword == geom_codes[2]) /// bug #135. Lajos, 08/01/2005
356 tword = getwrd(FALSE, &eol, &(geom_codes[8]), 2); /// bug #135. Lajos, 08/01/2005
359 fprintf(output, "\n **** keyword arc or tangent is missing");
360 fprintf(output, "\n **** program aborted\n");
364 else /// bug #135. Lajos, 08/01/2005
366 // assuming PARALLEL VARIABLE LINE when LINOGRAM has been specified
368 Creacm.word1 = geom_codes[1];
369 tword = geom_codes[5];
372 Creacm.word2 = tword;
374 GeoPar.par = (Creacm.word1 == geom_codes[1]);
375 GeoPar.div = (Creacm.word1 == geom_codes[2]); // moved out of the "if(!GeoPar.par) {...}" block. Lajos, 08/02/2005
380 // DIVERGENT GEOMETRY READ SPECS
383 GeoPar.arc = (Creacm.word2 == geom_codes[8]); /// bug #135. Lajos, 08/01/2005
384 GeoPar.tang = (Creacm.word2 == geom_codes[9]); /// bug #135. Lajos, 08/01/2005
385 if (!(GeoPar.div && (GeoPar.arc || GeoPar.tang)))
388 fprintf(output, "\n rays are divergent from point sources");
390 GeoPar.radius = getnum(FALSE, &eol);
391 GeoPar.stod = getnum(FALSE, &eol);
392 fprintf(output, "\n source to origin distance %10.4f",
398 "\n the detectors lie on an arc with source to detector distance = %10.4f",
405 "\n the detectors lie on a st. line with source to middle detector distance = %10.4f",
411 if (GeoPar.radius <= GeoPar.nelem * GeoPar.pixsiz / 1.414)
417 // PARALLEL GEOMETRY READ SPECS; IS SPACING UNIFORM OR VARIABLE
419 GeoPar.uni = (Creacm.word2 == geom_codes[4]); /// bug #135. Lajos, 08/01/2005
420 GeoPar.vri = (Creacm.word2 == geom_codes[5]); /// bug #135. Lajos, 08/01/2005
423 tword = geom_codes[6]; /// bug #135. Lajos, 08/01/2005
426 tword = getwrd(FALSE, &eol, &(geom_codes[6]), 2); /// bug #135. Lajos, 08/01/2005
429 fprintf(output, "\n **** keyword line or strip is missing");
430 fprintf(output, "\n **** program aborted\n");
434 Creacm.word3 = tword;
436 GeoPar.line = (Creacm.word3 == geom_codes[6]); /// bug #135. Lajos, 08/01/2005
437 GeoPar.strip = (Creacm.word3 == geom_codes[7]); /// bug #135. Lajos, 08/01/2005
439 if (!(GeoPar.par && (GeoPar.strip || GeoPar.line)))
441 if (!(GeoPar.uni || GeoPar.vri))
444 fprintf(output, "\n rays are parallel with");
447 fprintf(output, " uniform spacing between rays");
451 fprintf(output, " variable spacing between rays");
455 fprintf(output, "\n data collected along strips\n");
459 fprintf(output, "\n data collected along lines\n");
462 fprintf(output, "\n");
465 void InputFile_class::ReadNumberOfRaysAndRaySpacing(INTEGER invoke)
468 INTEGER melen, tword;
471 static const INTEGER proj = CHAR2INT('p', 'r', 'o', 'j');
472 static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
474 static const INTEGER rays_codes[3] =
475 { CHAR2INT('r', 'a', 'y', 's'), CHAR2INT('u', 's', 'e', 'r'), CHAR2INT('p',
478 if (GeoPar.lino) /// bug #135. Lajos, 08/01/2005
480 // assuming RAYS USER if LINOGRAM has been specified
481 if (GeoPar.nelem == 0)
484 "\n **** nelem must be specified using phantom average if linogram geometry has been chosen");
485 fprintf(output, "\n **** program aborted\n");
488 tword = rays_codes[1];
489 GeoPar.usrays = 2 * GeoPar.nelem + 1;
490 GeoPar.pinc = GeoPar.pixsiz;
494 tword = getwrd(TRUE, &eol, rays_codes, 1);
496 creaer(rays_codes[0]);
500 tword = getwrd(FALSE, &eol, &(rays_codes[1]), 2);
504 if (tword == rays_codes[1])
507 GeoPar.usrays = getint(FALSE, &eol);
509 if (!((GeoPar.usrays == 0) && (invoke == crea)))
511 if ((GeoPar.usrays % 2) <= 0)
517 melen = getint(FALSE, &eol);
518 zisxip = getnum(FALSE, &eol);
521 GeoPar.pinc = getnum(FALSE, &eol);
522 if (GeoPar.pinc < Consts.zero)
526 if ((invoke == proj) && (GeoPar.nelem == 0))
530 if ((GeoPar.pinc * GeoPar.usrays - 1.0) / GeoPar.stod > Consts.pi)
534 if (tword == rays_codes[2])
535 GeoPar.usrays = GeoPar.numray(melen, zisxip);
537 if (GeoPar.nelem > 0)
538 GeoPar.snrays = GeoPar.numray(GeoPar.nelem, GeoPar.pixsiz);
539 if (GeoPar.nelem == 0)
540 GeoPar.snrays = GeoPar.usrays;
541 if (GeoPar.usrays == 0)
542 GeoPar.usrays = GeoPar.snrays;
544 fprintf(output, "\n number of rays per projection %5i",
549 fprintf(output, "\n snark computed number of rays %5i",
552 fprintf(output, "\n at detector spacing %10.4f", GeoPar.pinc);
554 // SET THE APPROPRIATE VALUES FOR THE FIRST AND LAST RAY NUMBERS
555 // FOR THE USER RAY SET AND THE SNARK RAY SET
556 // THE VALUE OF NRAYS IN /USER/ COMMON NOW IS NO OF USER SPECIFIED
557 // RAYS IF GIVEN OR SNRAYS COMPUTED BY SNARK IF DEFAULT
559 GeoPar.nrays = MAX0(GeoPar.snrays, GeoPar.usrays);
560 GeoPar.midray = GeoPar.nrays / 2;
562 GeoPar.fsnray = (GeoPar.nrays - GeoPar.snrays) / 2;
563 GeoPar.lsnray = GeoPar.fsnray + GeoPar.snrays - 1;
565 GeoPar.fusray = (GeoPar.nrays - GeoPar.usrays) / 2;
566 GeoPar.lusray = GeoPar.fusray + GeoPar.usrays - 1;
569 // MEASUREMENT STATISTICS -- NOISE AND SCATTER
571 void InputFile_class::GetMesurementStatistics()
574 BOOLEAN perfct, noisfl;
578 static const INTEGER hquan = CHAR2INT('q', 'u', 'a', 'n');
579 static const INTEGER hscat = CHAR2INT('s', 'c', 'a', 't');
580 static const INTEGER hmult = CHAR2INT('m', 'u', 'l', 't');
581 static const INTEGER haddi = CHAR2INT('a', 'd', 'd', 'i');
583 static const INTEGER meas_codes[10] =
585 CHAR2INT('m', 'e', 'a', 's'), CHAR2INT('p', 'e', 'r', 'f'),
586 CHAR2INT('n', 'o', 'i', 's'), CHAR2INT('s', 'e', 'e', 'd'),
587 CHAR2INT('b', 'a', 'c', 'k'), CHAR2INT('q', 'u', 'a', 'n'),
588 CHAR2INT('s', 'c', 'a', 't'), CHAR2INT('m', 'u', 'l', 't'),
589 CHAR2INT('a', 'd', 'd', 'i'), CHAR2INT('c', 'a', 'l', 'i')
592 NoisePar.setQuantumNoise(0, 1.0, 1.0); //bug221 - swr - 03/03/07
594 word = getwrd(TRUE, &eol, meas_codes, 1);
596 creaer(meas_codes[0]);
598 word = getwrd(FALSE, &eol, &(meas_codes[1]), 2);
601 fprintf(output, "\n **** keyword perfect or noisy is missing");
602 fprintf(output, "\n **** program aborted\n");
606 perfct = (word == meas_codes[1]);
607 noisfl = (word == meas_codes[2]);
611 // IF PERFECT MEASUREMENTS THEN SKIP NOISE INFORMATION
613 creaer(meas_codes[1]);
615 fprintf(output, "\n projection data are noiseless\n");
617 word = getwrd(TRUE, &eol, &(meas_codes[4]), 1);
624 "\n noise characteristics of projection data follow");
625 fprintf(output, "\n nature characteristics");
629 // CHECK TYPE OF NOISE AND READ RELEVANT DATA
630 word = getwrd(TRUE, &eol, &(meas_codes[3]), 6);
632 if ((word == meas_codes[3]) || (word == meas_codes[4]))
640 NoisePar.quanmn = getnum(FALSE, &eol);
641 NoisePar.quancm = getnum(FALSE, &eol);
643 // GET CALIBRATION TYPE AND ECHO ON OUTPUT
645 if (getwrd(FALSE, &eol, &(meas_codes[9]), 1) == meas_codes[9])
646 NoisePar.quanin = getint(FALSE, &eol);
647 if (NoisePar.quanin == 4)
649 fprintf(output, "\n Emission tomography");
654 if ((NoisePar.quanin < 1) || (NoisePar.quanin > 4))
657 "\n quantum mean number of photons in actual measurement %12.4e",
660 "\n mean number of photons in calibration measurement is %10.4f times above",
663 if (NoisePar.quanin == 1)
666 "\n calibration is constant for all rays in any single projection");
669 if (NoisePar.quanin == 2)
672 "\n calibration is constant over all projections for any single ray");
675 if (NoisePar.quanin == 3)
678 "\n calibration is variable for every ray in every projection");
681 if ((NoisePar.quanmn < Consts.zero)
682 || (NoisePar.quancm < Consts.zero))
685 // DISTRIBUTE PHOTONS ACCORDING TO APERTURE WEIGHTS COMPUTED EARLIER
687 for (n = 0; n < GeoPar.nave2; n++)
689 Creacm.aperwt[n] *= NoisePar.quanmn;
693 "\n mean numbers of photons entering each substrip are");
694 fprintf(output, "\n ");
695 for (n = 0; n < GeoPar.nave2; n++)
697 fprintf(output, "%12.4e ", Creacm.aperwt[n]);
700 NoisePar.setQuantumNoise(); //bug221 - swr - 03/03/07
705 ///if(word == hscat) {
707 //bug221 - swr - 03/03/07
709 REAL peak = getnum(FALSE, &eol);
710 REAL width = getnum(FALSE, &eol);
711 NoisePar.setScatter(peak, width);
713 "\n scatter the fraction scattered in line is %10.4f",
716 "\n scatter extends to a distance of %10.4f",
719 "\n dying linearly from the middle");
721 if ((peak < Consts.zero) || (width < Consts.zero))
728 // MULTIPLICATIVE NOISE
729 //bug221 - swr - 03/03/07
731 REAL mean = getnum(FALSE, &eol);
732 REAL std = getnum(FALSE, &eol);
733 NoisePar.setMultiplicativeNoise(mean, std);
735 "\n multiplicative mean %10.4f std. dev. %10.4f",
737 if (fabs(NoisePar.ultnmn) < Consts.zero)
744 //bug221 - swr - 03/03/07
746 REAL mean = getnum(FALSE, &eol);
747 REAL std = getnum(FALSE, &eol);
748 NoisePar.setAdditiveNoise(mean, std);
750 " additive %10.4f std. dev. %10.4f",
759 noisfl = ((NoisePar.quanin > 0) || NoisePar.ultnfl || NoisePar.addnfl);
763 if (word != meas_codes[3])
764 creaer(meas_codes[3]);
766 // INPUT SEED FOR RANDOM NUMBER GENERATOR
768 Creacm.iseed = getint(FALSE, &eol);
770 "\n seed for random number generator is %10i\n",
773 word = getwrd(TRUE, &eol, &(meas_codes[4]), 1);
778 // READ NUMBER OF PROJECTIONS
779 // READ LIST OF PROJECTION ANGLES INTO WORK AREA
781 void InputFile_class::ReadProjectionAngles(BOOLEAN f11)
784 INTEGER JJ, IANG, NUMRAYS, II, LIMIT, NN;
794 static const INTEGER angl_codes[2] = /// bug #135. Lajos, 08/01/2005
795 { CHAR2INT('a', 'n', 'g', 'l'), CHAR2INT('e', 'q', 'u', 'a') /// bug #135. Lajos, 08/01/2005
796 // CHAR2INT('l','i','n','o') /// bug #135. Lajos, 08/01/2005
799 fprintf(output, "\n");
801 if (GeoPar.lino) /// bug #135. Lajos, 08/01/2005
803 GeoPar.prjnum = 2 * (2 * GeoPar.nelem + 1);
807 tword = getwrd(TRUE, &eol, angl_codes, 1);
809 creaer(angl_codes[0]);
811 GeoPar.prjnum = getint(FALSE, &eol);
813 fprintf(output, "\n total number of projections %5i",
816 if (GeoPar.prjnum <= 0)
819 // FIRST ALLOCATE STORAGE AREA
820 if (Creacm.pang != NULL)
821 delete[] Creacm.pang;
822 Creacm.pang = new REAL[GeoPar.prjnum];
823 Anglst.Init(GeoPar.prjnum);
825 // WANT TO GENERATE EQUALLY-SPACED ANGLES ONLY IF SPECIFIED
832 linosp = TRUE; /// bug #135. Lajos, 08/01/2005
835 tword = getwrd(FALSE, &eol, &(angl_codes[1]), 1); /// bug #135. Lajos, 08/01/2005
836 if (tword == angl_codes[1])
850 *Creacm.pang = getnum(FALSE, &eol);
852 creaer(25); // added error checking. Lajos, 03/29/2005
853 *(Creacm.pang + 1) = getnum(FALSE, &eol);
855 creaer(25); // added error checking. Lajos, 03/29/2005
856 //C intang = intang + 1
858 // CHANGED CODE FOR READING EQUALLY SPACED ANGLES TO MAKE USERS
859 // LIFE BIT MORE COMFORTABLE FOR VERSION 1.02 SNARK87 10th AUG
860 // FOR flag DESCRIPTION SEE CODE AT LABEL 40
861 // -----S.SADANANDA AUG 10 1987
863 if (GeoPar.prjnum <= 1)
866 delta = (Creacm.pang[1] - Creacm.pang[0]) / (GeoPar.prjnum - 1);
868 fprintf(output,"\n %36.30f", (double) delta);
870 for (np = 1; np < GeoPar.prjnum; np++)
872 Creacm.pang[np] = np * delta + Creacm.pang[0];
874 fprintf(output,"\n %5i %36.30f", np+1, Creacm.pang[np]);
879 // Modified code to enable generation of angles for linogram
881 // ----- J. Browne, May 1993
885 if (!(GeoPar.par && GeoPar.vri))
888 "\n **** for linogram reconstruction use PARALLEL geometry with VARIABLE detector spacing");
889 fprintf(output, "\n **** program aborted\n");
892 NN = (GeoPar.prjnum / 2 - 3) / 4;
894 NUMRAYS = 4 * NN + 3;
895 RAD_DEG = (REAL) 180.0 / Consts.pi;
897 for (II = 1; II <= 2; II++)
899 for (IANG = -LIMIT; IANG <= LIMIT; IANG++)
904 Creacm.pang[JJ] = RAD_DEG
905 * (REAL) atan((REAL) 2.0 * IANG / NUMRAYS) + 90;
910 Creacm.pang[JJ] = RAD_DEG
911 * (REAL) atan((REAL) 2.0 * IANG / NUMRAYS) + 180;
925 fprintf(output, "\n");
926 //bug 276, jklukowska, allow free format in file11
927 for (i = 0; i < GeoPar.prjnum; i++)
929 Creacm.pang[i] = getnum(FALSE, &eol); // try to get number
932 getnxt(TRUE, FALSE); // get another line
933 Creacm.pang[i] = getnum(FALSE, &eol); // try again
943 for (i = 0; i < GeoPar.prjnum; i++)
945 Creacm.pang[i] = getnum(FALSE, &eol); // try to get number
948 getnxt(TRUE); // get another line
949 Creacm.pang[i] = getnum(FALSE, &eol); // try again
955 fprintf(output, "\n projection angles");
957 for (j = 0; j < GeoPar.prjnum; j++)
959 if (((j % 10) == 0) && (j > 0))
960 fprintf(output, "\n ");
961 fprintf(output, " %9.4f", Creacm.pang[j]);
964 // CONVERT TO RADIANS AND COMPUTE SIN AND COS TABLES
966 Anglst.InitDiv(Creacm.pang, GeoPar.prjnum);
968 fprintf(output, "\n");
971 void InputFile_class::ReadIfoForComputionRaysums()
976 GeoPar.nave2 = getint(FALSE, &eol);
978 "\n projection data are calculated by dividing each ray interval into %i substrips",
981 if ((GeoPar.nave2 <= 0) || (GeoPar.nave2 > 13)
982 || (GeoPar.nave2 == GeoPar.nave2 / 2 * 2))
986 /////////////////////////////////////////////////////////////
987 for (i = 0; i < GeoPar.nave2; i++)
988 { /// Is this what is needed MK ???
989 GeoPar.naper[i] = getint(FALSE, &eol);
992 // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
997 fprintf(output, "\n with aperture (substrip) weights");
998 for (n = 0; n < GeoPar.nave2; n++)
1000 fprintf(output, " %5i", GeoPar.naper[n]);
1003 // GET SUM OF THE APERTURE WEIGHTS FOR LATER USE
1005 for (n = 0; n < GeoPar.nave2; n++)
1007 if (GeoPar.naper[n] < 0)
1009 ns += GeoPar.naper[n];
1014 // COMPUTE NORMALIZED APERTURE WEIGHTS FOR USE IN CREATING
1016 for (n = 0; n < GeoPar.nave2; n++)
1018 Creacm.aperwt[n] = ((REAL) GeoPar.naper[n]) / ((REAL) (ns));
1022 // READ OBJECTS MAKING UP THE PHANTOM
1024 void InputFile_class::ReadObjects(INTEGER invoke)
1031 INTEGER itype, tword;
1043 SNARK_Object_Node *f_object, *t_object;
1045 static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
1047 static const INTEGER object_codes[8] =
1049 CHAR2INT('o', 'b', 'j', 'e'), CHAR2INT('e', 'l', 'i', 'p'),
1050 CHAR2INT('r', 'e', 'c', 't'), CHAR2INT('t', 'r', 'i', 'a'),
1051 CHAR2INT('s', 'e', 'g', 'm'), CHAR2INT('s', 'e', 'c', 't'),
1052 CHAR2INT('l', 'a', 's', 't'), CHAR2INT('d', 'e', 'n', 's') };
1054 tword = getwrd(TRUE, &eol, object_codes, 1);
1055 if (tword != object_codes[0])
1056 creaer(object_codes[0]);
1058 fprintf(output, "\n description of objects");
1059 fprintf(output, "\n density at levels");
1060 fprintf(output, "\n numb type x-coord y-coord x-length y-length angle av dens");
1061 for (i = 0; i < Spctrm.nergy; i++)
1062 fprintf(output, " %9i", Spctrm.energy[i]);
1064 // READ IN OBJECT PARAMETERS IN TEMPORARY LOCATIONS SO AS NOT TO
1065 // CLOBBER WORK AREA WHILE DOING INPUT FOR RDPICT AND RDPROJ
1069 f_object = new SNARK_Object_Node;
1070 t_object = f_object;
1075 getnxt(TRUE); // read line of input
1076 itype = getwrd(FALSE, &eol, &(object_codes[1]), 6);
1080 fprintf(output, "\n **** keyword elip, rect, tria, segm, sect or last is missing");
1081 fprintf(output, "\n **** program aborted\n");
1085 // CHECK TO SEE IF LAST OBJECT HAS BEEN READ
1087 if (itype == object_codes[6])
1090 cxt = getnum(FALSE, &eol);
1093 fprintf(output, "\n **** parameter cx of object %s is missing", int2str(itype));
1094 fprintf(output, "\n **** program aborted\n");
1098 cyt = getnum(FALSE, &eol);
1101 fprintf(output, "\n **** parameter cy of object %s is missing", int2str(itype));
1102 fprintf(output, "\n **** program aborted\n");
1106 ut = getnum(FALSE, &eol);
1109 fprintf(output, "\n **** parameter u of object %s is missing", int2str(itype));
1110 fprintf(output, "\n **** program aborted\n");
1114 vt = getnum(FALSE, &eol);
1117 fprintf(output, "\n **** parameter v of object %s is missing", int2str(itype));
1118 fprintf(output, "\n **** program aborted\n");
1122 angt = getnum(FALSE, &eol);
1125 fprintf(output, "\n **** parameter ang of object %s is missing", int2str(itype));
1126 fprintf(output, "\n **** program aborted\n");
1130 dent[0] = getnum(FALSE, &eol);
1133 fprintf(output, "\n **** parameter den(0) of object %s is missing", int2str(itype));
1134 fprintf(output, "\n **** program aborted\n");
1138 // IF POLYCHROMATIC READ IN ABSORPTION FOR OTHER LEVELS
1140 if (Spctrm.nergy > 1)
1142 getnxt(TRUE); // read line of input
1143 getwrd(FALSE, &eol, &(object_codes[7]), 1);
1147 fprintf(output, "\n **** keyword dens is missing");
1148 fprintf(output, "\n **** program aborted\n");
1152 for (j = 1; j < Spctrm.nergy; j++)
1154 dent[j] = getnum(FALSE, &eol);
1157 fprintf(output, "\n **** parameter den(%d) of object %s is missing", j, int2str(itype));
1158 fprintf(output, "\n **** program aborted\n");
1164 // ALSO COMPUTE THE WEIGHTED AVERAGE ABSORPTION
1167 for (j = 0; j < Spctrm.nergy; j++)
1168 densit += Spctrm.engwt[j] * dent[j];
1170 fprintf(output, "\n %4i %s %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f", i + 1, int2str(itype), cxt, cyt, ut, vt, angt, densit);
1171 for (j = 0; j < Spctrm.nergy; j++)
1172 fprintf(output, " %9.4f", dent[j]);
1174 // IDENTIFY OBJECT TYPE
1176 for (otype = SOT_elip; otype <= SOT_sect; otype++)
1178 if (itype == object_codes[otype + 1])
1182 // IF INVOKED BY RDPICT OR RDPROJ NO NEED TO REMEMBER THE ABOVE VALUE
1185 if ((ut < Consts.zero) || (vt < Consts.zero))
1188 // TRANSFER VALUES TO A TEMPORARY LIST
1190 t_object->data.type = (SNARK_Object_type) otype;
1191 for (j = 0; j < Spctrm.nergy; j++)
1192 t_object->data.den1[j] = dent[j];
1193 t_object->data.denav = densit;
1194 t_object->data.cx = cxt;
1195 t_object->data.cy = cyt;
1196 t_object->data.u = ut;
1197 t_object->data.v = vt;
1198 t_object->data.ang = angt;
1199 phi = (REAL) (angt * Consts.pi / 180.0);
1200 t_object->data.sinang = (REAL) sin(phi);
1201 t_object->data.cosang = (REAL) cos(phi);
1203 t_object->next = new SNARK_Object_Node;
1204 t_object = t_object->next;
1210 MAX_NUMBER_OF_OBJECTS = i;
1211 objects = new SNARK_Object[MAX_NUMBER_OF_OBJECTS];
1213 t_object = f_object;
1214 for (i = 0; i < MAX_NUMBER_OF_OBJECTS; i++)
1216 objects[i].type = t_object->data.type;
1217 for (j = 0; j < Spctrm.nergy; j++)
1218 objects[i].den1[j] = t_object->data.den1[j];
1219 objects[i].denav = t_object->data.denav;
1220 objects[i].cx = t_object->data.cx;
1221 objects[i].cy = t_object->data.cy;
1222 objects[i].u = t_object->data.u;
1223 objects[i].v = t_object->data.v;
1224 objects[i].ang = t_object->data.ang;
1225 objects[i].sinang = t_object->data.sinang;
1226 objects[i].cosang = t_object->data.cosang;
1228 t_object = t_object->next;
1230 f_object = t_object;
1234 // GRAB THE SCALE FACTOR FROM THE LAST CARD
1236 scale = getnum(FALSE, &eol);
1239 fprintf(output, "\n **** scale factor is missing");
1240 fprintf(output, "\n **** program aborted\n");
1245 "\n scale factor multiplying object densities %10.4f",
1249 if (fabs(scale) < Consts.zero)
1252 // IF PICT OR PROJ NO NEED TO SCALE
1256 // GET THE NUMBER OF OBJECTS READ IN
1259 if (Creacm.nobj != 0)
1261 // SCALE DENSITIES FIRST
1262 for (n = 0; n < Creacm.nobj; n++)
1264 objects[n].denav *= scale;
1265 for (j = 0; j < Spctrm.nergy; j++)
1266 objects[n].den1[j] *= scale;
1270 fprintf(output, "\n");
1272 //jk 11/18/2007 introducing variability
1273 //GRAB THE VARIABILITY VALUE FROM THE LAST CARD
1274 //jklukowska 12/9/8, bug 274, varseed should be read as an integer, not real number
1275 Creacm.varseed = getint(FALSE, &eol);
1281 fprintf(output, "\n seed set to %d", Creacm.varseed);
1284 variability = getnum(FALSE, &eol);
1289 else if (variability < 0)
1291 fprintf(output, "\n **** inhomogeneity should be >= 0");
1292 fprintf(output, "\n **** program aborted\n");
1296 if (variability != 0)
1298 else if (variability == 0)
1301 fprintf(output, "\n inhomogeneity set to %10.4f", variability);
1305 // added by Lajos, Jan 25, 2005
1307 ECHOLINE PRINTS THE CONTENTS OF THE CURRENT INPUT LINE (I.E. THE ARRAY 'DATA') TO THE OUTPUT.
1308 IF flag == true, THE OUTPUT WILL START WITH '<*>'. OTHERWISE, '<#>' WILL BE USED.
1310 void InputFile_class::echoline(BOOLEAN flag)
1313 fprintf(output, "\n <*> %s", data);
1315 fprintf(output, "\n <#> %s", data);