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/file11.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
11 #define _LARGEFILE_SOURCE
12 #define _FILE_OFFSET_BITS 64
33 void File11_class::CreaInPicture()
36 static const INTEGER pict = CHAR2INT('p', 'i', 'c', 't');
40 // READ ENERGY SPECTRAL DISTRIBUTION
41 ReadEnergySpectralDistribution();
43 // READ OBJECTS MAKING UP THE PHANTOM
46 // GET INFORMATION ABOUT CREATION OF PHANTOM
47 GetInformationABoutCreationOfPhantom(pict);
52 void File11_class::CreaInProj()
55 static const INTEGER proj = CHAR2INT('p', 'r', 'o', 'j');
59 // READ ENERGY SPECTRAL DISTRIBUTION
60 ReadEnergySpectralDistribution();
62 // READ OBJECTS MAKING UP THE PHANTOM
65 // READ RAYSUM INFORMATION
66 ReadRaysumInformation();
73 // READ IN AVERAGING INFORMATION FOR COMPUTING RAYSUMS
74 ReadIfoForComputionRaysums();
76 // INPUT ENTIRE GEOMETRY DATA
77 ReadGeometryData(proj);
79 // READ NUMBER OF RAYS AND RAY SPACING
80 ReadNumberOfRaysAndRaySpacing(proj);
82 // READ NUMBER OF PROJECTIONS
83 // READ LIST OF PROJECTION ANGLES INTO WORK AREA
84 ReadProjectionAngles(TRUE);
86 // MEASUREMENT STATISTICS -- NOISE AND SCATTER
87 GetMesurementStatistics();
89 // READ BACKGROUND ABSORPTION
90 ReadBackgroundAbsorption();
94 bool File11_class::ReadProjection(REAL* beta, REAL* ibase, INTEGER usrays)
98 if(fscanf(infile, "%lf", &val) != 1)
103 // skip angle in degrees
104 if(fscanf(infile, "%lf", &val) != 1)
107 for(int j = 0; j < usrays; j++){
108 if(fscanf(infile, "%lf", &val) != 1){
117 INTEGER File11_class::Open(BOOLEAN wr)
121 infile = fopen("file11", "wb");
125 infile = fopen("file11", "rb");
130 return -1; // error opening file
135 void File11_class::Close()
137 fclose(infile); // !!! core dump under cygwin
140 void File11_class::Flush()
145 void File11_class::Rewind()
147 fseek(infile, 0, SEEK_SET);
150 void File11_class::WriteRowOfPhanPixels(REAL* test_aray, INTEGER nelem)
153 for (int k = 0; k < nelem; k++)
157 fprintf(infile, "\n");
160 fprintf(infile, " %25.20f", test_aray[k]);
162 if (test_aray[k] != 0)
169 void File11_class::WriteEndOfPhan()
171 fprintf(infile, "\n");
174 void File11_class::WriteProj(REAL theta, REAL* nlo, INTEGER usrays)
176 REAL thedeg = theta * 180 / Consts.pi;
179 fprintf(infile, "\n %25.20f %25.20f", theta, thedeg);
181 for (int nr = 0; nr < usrays; nr++)
185 fprintf(infile, "\n");
187 fprintf(infile, " %25.20f", nlo[nr]);
191 void File11_class::WriteEndOfProj()
193 fprintf(infile, "\n");
196 void File11_class::ReadPicture(REAL* Pict, INTEGER nelem)
199 signed CHAR ch; // added by Lajos, Nov 18, 2004
201 //#pragma message( "to fix!" )
203 for(int i = 0; i < nelem; i++){
206 for(int j = 0; j < nelem; j++){
207 if(fscanf(infile, "%lf", &val) != 1){
208 fprintf(output,"\n **** end of file encountered on file file11");
209 fprintf(output, "\n **** program aborted\n");
216 while (((ch = fgetc(infile)) != '\n') && (ch != EOF))
217 { // added check for EOF. Lajos, Nov 18, 2004
223 void File11_class::wrhdr(INTEGER word)
226 static const INTEGER keywd[40] =
228 CHAR2INT('l', 'a', 's', 't'),
229 CHAR2INT('s', 'p', 'e', 'c'), CHAR2INT('m', 'o', 'n', 'o'),
230 CHAR2INT('p', 'o', 'l', 'y'), CHAR2INT('o', 'b', 'j', 'e'),
231 CHAR2INT('p', 'h', 'a', 'n'), CHAR2INT('a', 'v', 'e', 'r'),
232 CHAR2INT('r', 'a', 'y', 's'), CHAR2INT('g', 'e', 'o', 'm'),
233 CHAR2INT('p', 'a', 'r', 'a'), CHAR2INT('d', 'i', 'v', 'e'),
234 CHAR2INT('s', 't', 'r', 'i'), CHAR2INT('l', 'i', 'n', 'e'),
235 CHAR2INT('t', 'a', 'n', 'g'), CHAR2INT('a', 'r', 'c', ' '),
236 CHAR2INT('u', 'n', 'i', 'f'), CHAR2INT('v', 'a', 'r', 'i'),
237 CHAR2INT('a', 'n', 'g', 'l'), CHAR2INT('m', 'e', 'a', 's'),
238 CHAR2INT('p', 'e', 'r', 'f'), CHAR2INT('n', 'o', 'i', 's'),
239 CHAR2INT('q', 'u', 'a', 'n'), CHAR2INT('s', 'c', 'a', 't'),
240 CHAR2INT('m', 'u', 'l', 't'), CHAR2INT('a', 'd', 'd', 'i'),
241 CHAR2INT('s', 'e', 'e', 'd'), CHAR2INT('b', 'a', 'c', 'k'),
242 CHAR2INT('r', 'u', 'n', ' '), CHAR2INT('m', 'o', 'd', 'i'),
243 CHAR2INT('s', 'a', 'v', 'e'), CHAR2INT('p', 'i', 'x', 'e'),
244 CHAR2INT('e', 'q', 'u', 'a'), CHAR2INT('u', 's', 'e', 'r'),
245 CHAR2INT('a', 'n', 'd', ' '), CHAR2INT('a', 'f', 't', 'e'),
246 CHAR2INT('s', 't', 'e', 'p'), CHAR2INT('s', 'p', 'a', 'c'),
247 CHAR2INT('p', 'n', 'c', 'h'), CHAR2INT('l', 'i', 'n', 'o'), /// bug #135. Lajos, 08/02/2005
248 CHAR2INT(' ', ' ', ' ', ' ')
251 static const CHAR keywdStr[40][5] =
252 { "last", "spec", "mono", "poly", "obje", "phan", "aver", "rays", "geom",
253 "para", "dive", "stri", "line", "tang", "arc ", "unif", "vari",
254 "angl", "meas", "perf", "nois", "quan", "scat", "mult", "addi",
255 "seed", "back", "run ", "modi", "save", "pixe", "equa", "user",
256 "and ", "afte", "step", "spac", "pnch", "lino", /// bug #135. Lajos, 08/02/2005
259 static const CHAR kindStr[5][5] =
260 { "elip", "rect", "tria", "segm", "sect" };
262 static const INTEGER KWI_last = 0;
263 static const INTEGER KWI_spec = 1;
264 static const INTEGER KWI_mono = 2;
265 static const INTEGER KWI_poly = 3;
266 static const INTEGER KWI_obj = 4;
267 static const INTEGER KWI_phan = 5;
268 static const INTEGER KWI_aver = 6;
269 static const INTEGER KWI_ray = 7;
270 static const INTEGER KWI_geom = 8;
271 static const INTEGER KWI_angl = 17;
272 static const INTEGER KWI_meas = 18;
273 static const INTEGER KWI_perf = 19;
274 static const INTEGER KWI_nois = 20;
275 static const INTEGER KWI_quan = 21;
276 static const INTEGER KWI_scat = 22;
277 static const INTEGER KWI_mult = 23;
278 static const INTEGER KWI_add = 24;
279 static const INTEGER KWI_seed = 25;
280 static const INTEGER KWI_back = 26;
281 static const INTEGER KWI_pixel = 30;
282 static const INTEGER KWI_lino = 38; /// bug #135. Lajos, 08/02/2005
284 // DATA SPACE FOR OBJECTS IN PHANTOM AT TOP OF BLANK COMMON
288 equivalence (object, work(2) )
296 // VARIABLES INDEXING THE KEYWORDS IN THE ARRAY KEYWD
299 integer last, spec, mono, poly, obj, phan, aver, ray, geom, parl,
300 1 divr, stri, lin, itang, iarc, unif, vari, angl, meas, perf,
301 2 nois, quan, scat, mult, add, seed, back, run, modi, save,
304 data last, spec, mono, poly, obj, phan, aver, ray, geom, parl/
305 1 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 /
306 data divr, stri, lin, itang, iarc, unif, vari, angl, meas, perf/
307 1 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20 /
308 data nois, quan, scat, mult, add, seed, back, run, modi, save /
309 1 21 , 22 , 23 , 24 , 25 , 26 , 27 , 28 , 29 , 30 /
310 data pixel, equ, pnch / 31, 32, 38/
314 integer den1, denav, type, cx, cy, u, v, ang, sinang, cosang, rayl
315 data den1,denav,type,cx,cy, u, v,ang,sinang,cosang,rayl
316 + / 1, 8, 9, 10,11,12,13,14, 15, 16, 17/
319 //static const REAL one = 1.0;
324 quanfl = NoisePar.quanin > 0;
325 noisfl = quanfl || NoisePar.ultnfl || NoisePar.addnfl;
326 perfct = !(NoisePar.sctnfl || noisfl);
328 // SET THE OUTPUT FILE TO FILE11 OR PUNCH
332 if(word == keywd[KWI_pnch]) {
333 // ALSO SET VALUES PROPERLY FOR UNAVAILABLE QUANTITIES
337 Spctrm.energy[0] = -1;
341 // WRITE HEADER INFORMATION ON FILE11 FOR PHANTOM & RAYSUMS
342 // PUNCH OUT THE HEADER INFORMATION FOR THE PHANTOM AND PICTURES
344 fprintf(infile, "%s", Creacm.name);
347 if (Spctrm.nergy == 1)
349 fprintf(infile, "\n%s %s %4i", keywdStr[KWI_spec],
350 keywdStr[KWI_mono], Spctrm.energy[0]);
353 if (Spctrm.nergy > 1)
356 fprintf(infile, "\n%s %4s %4i", keywdStr[KWI_spec],
357 keywdStr[KWI_poly], Spctrm.nergy);
358 fprintf(infile, "\n");
360 for (i = 0; i < Spctrm.nergy; i++)
362 fprintf(infile, "%5i %5i ", Spctrm.energy[i], Creacm.percnt[i]);
368 fprintf(infile, "\n%s", keywdStr[KWI_obj]);
370 if (Creacm.nobj != 0)
373 for (i = 0; i < Creacm.nobj; i++)
375 itype = (int) objects[i].type;
377 fprintf(infile, "\n%4s %9.4f %9.4f %9.4f %9.4f %9.4f %25.18f",
378 kindStr[itype], objects[i].cx, objects[i].cy, objects[i].u,
379 objects[i].v, objects[i].ang, objects[i].den1[0]);
381 if (Spctrm.nergy > 1)
384 fprintf(infile, "\ndens ");
385 for (j = 1; j < Spctrm.nergy; j++)
387 fprintf(infile, " %15.10f", objects[i].den1[j]);
393 //jk 11/18/2007 introducing variability
394 //the additional values from LAST card have to be printed in file11
395 fprintf(infile, "\n%s %10.4f %ld %10.4f", keywdStr[KWI_last], mean,
396 Creacm.varseed, variability);
398 // SKIP PHANTOM AVERAGING AND OTHER INFO SUCH AS NELEM FOR RDPROJ
400 if (word != keywd[KWI_ray])
402 if (word == keywd[KWI_phan])
405 fprintf(infile, "\n%4s %4s %4i", keywdStr[KWI_phan],
406 keywdStr[KWI_aver], Creacm.nave1);
407 fprintf(infile, "\n%4s %5i size %10.4f",
408 keywdStr[KWI_pixel], GeoPar.nelem, GeoPar.pixsiz);
414 // HEADER INFO FOR RAYSUMS
417 fprintf(infile, "\n%4s %4s %4i\n", keywdStr[KWI_ray], keywdStr[KWI_aver],
419 for (i = 0; i < GeoPar.nave2; i++)
421 fprintf(infile, " %4i", GeoPar.naper[i]);
424 fprintf(infile, "\n%4s", keywdStr[KWI_geom]);
426 if (GeoPar.lino) /// bug #135. Lajos, 08/02/2005
428 fprintf(infile, "\n%s", keywdStr[KWI_lino]);
434 fprintf(infile, "\n%4s", int2str(Creacm.word1));
435 fprintf(infile, " %4s", int2str(Creacm.word2));
436 fprintf(infile, " %4s", int2str(Creacm.word3));
443 strcpy(wrd2, int2str(Creacm.word2));
445 "\n%4s %4s source at%10.4f det dist%10.4f",
446 int2str(Creacm.word1), wrd2, GeoPar.radius, GeoPar.stod);
449 fprintf(infile, "\n%s user %5i spacing %10.4f",
450 keywdStr[KWI_ray], GeoPar.usrays, GeoPar.pinc);
453 fprintf(infile, "\n%s %5i\n", keywdStr[KWI_angl], GeoPar.prjnum);
454 for (i = 0; i < GeoPar.prjnum; i++)
456 if (((i % 7) == 0) && (i != 0))
458 fprintf(infile, "\n");
460 fprintf(infile, " %9.4f", Creacm.pang[i]);
465 // MEASUREMENT STATISTICS
469 fprintf(infile, "\n%s %s", keywdStr[KWI_meas], keywdStr[KWI_perf]);
475 fprintf(infile, "\n%s %s", keywdStr[KWI_meas], keywdStr[KWI_nois]);
480 fprintf(infile, "\n%s %14.4f %14.4f cali %2i",
481 keywdStr[KWI_quan], NoisePar.quanmn, NoisePar.quancm,
487 fprintf(infile, "\n%s noise peak %20.10f width %20.10f",
488 keywdStr[KWI_scat], NoisePar.sctnpk, NoisePar.sctnwd);
493 fprintf(infile, "\n%s noise mean %20.10f std dev %20.10f",
494 keywdStr[KWI_mult], NoisePar.ultnmn, NoisePar.ultnsd);
499 fprintf(infile, "\n%s noise mean %20.10f std dev %20.10f",
500 keywdStr[KWI_add], NoisePar.addnmn, NoisePar.addnsd);
503 if (!perfct && noisfl)
505 fprintf(infile, "\n%s %4i", keywdStr[KWI_seed], Creacm.iseed);
508 fprintf(infile, "\n%4s ", keywdStr[KWI_back]);
509 for (i = 0; i < Spctrm.nergy; i++)
511 fprintf(infile, " %7.4f", Spctrm.backgr[i]);