src/snark/file11.cpp

   1 /*
   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) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #define _LARGEFILE_SOURCE
  12 #define _FILE_OFFSET_BITS 64
  13
  14 #include <cstdlib>
  15 #include <cstdio>
  16 #include <cstring>
  17
  18 #include "blkdta.h"
  19 #include "creacm.h"
  20 #include "geom.h"
  21 #include "spctrm.h"
  22 #include "noise.h"
  23 #include "objects.h"
  24 #include "consts.h"
  25 #include "uiod.h"
  26
  27 #include "int2str.h"
  28
  29 #include "file11.h"
  30
  31 File11_class File11;
  32
  33 void File11_class::CreaInPicture()
  34 {
  35
  36         static const INTEGER pict = CHAR2INT('p', 'i', 'c', 't');
  37
  38         GetName();
  39
  40         // READ ENERGY SPECTRAL DISTRIBUTION
  41         ReadEnergySpectralDistribution();
  42
  43         // READ OBJECTS MAKING UP THE PHANTOM
  44         ReadObjects(pict);
  45
  46         // GET INFORMATION ABOUT CREATION OF PHANTOM
  47         GetInformationABoutCreationOfPhantom(pict);
  48
  49         return;
  50 }
  51
  52 void File11_class::CreaInProj()
  53 {
  54
  55         static const INTEGER proj = CHAR2INT('p', 'r', 'o', 'j');
  56
  57         GetName();
  58
  59         // READ ENERGY SPECTRAL DISTRIBUTION
  60         ReadEnergySpectralDistribution();
  61
  62         // READ OBJECTS MAKING UP THE PHANTOM
  63         ReadObjects(proj);
  64
  65         // READ RAYSUM INFORMATION
  66         ReadRaysumInformation();
  67
  68         if (!Creacm.prjflg)
  69         {
  70                 return;
  71         }
  72
  73         // READ IN AVERAGING INFORMATION FOR COMPUTING RAYSUMS
  74         ReadIfoForComputionRaysums();
  75
  76         // INPUT ENTIRE GEOMETRY DATA
  77         ReadGeometryData(proj);
  78
  79         // READ NUMBER OF RAYS AND RAY SPACING
  80         ReadNumberOfRaysAndRaySpacing(proj);
  81
  82         // READ NUMBER OF PROJECTIONS
  83         // READ LIST OF PROJECTION ANGLES INTO WORK AREA
  84         ReadProjectionAngles(TRUE);
  85
  86         // MEASUREMENT STATISTICS  --  NOISE AND SCATTER
  87         GetMesurementStatistics();
  88
  89         // READ BACKGROUND ABSORPTION
  90         ReadBackgroundAbsorption();
  91 }
  92
  93
  94 bool File11_class::ReadProjection(REAL* beta, REAL* ibase, INTEGER usrays)
  95 {
  96  REAL val;
  97
  98  if(fscanf(infile, "%lf", &val) != 1)
  99     return FALSE;
 100
 101  *beta = val;
 102
 103  // skip angle in degrees
 104  if(fscanf(infile, "%lf", &val) != 1)
 105     return FALSE;
 106
 107  for(int j = 0; j < usrays; j++){
 108      if(fscanf(infile, "%lf", &val) != 1){
 109         return FALSE;
 110        }
 111      *(ibase + j) = val;
 112     }
 113
 114  return TRUE;
 115 }
 116
 117 INTEGER File11_class::Open(BOOLEAN wr)
 118 {
 119         if (wr)
 120         {
 121                 infile = fopen("file11", "wb");
 122         }
 123         else
 124         {
 125                 infile = fopen("file11", "rb");
 126         }
 127
 128         if (infile == NULL)
 129         {
 130                 return -1; // error opening file
 131         }
 132         return 0;
 133 }
 134
 135 void File11_class::Close()
 136 {
 137         fclose(infile); // !!! core dump under cygwin
 138 }
 139
 140 void File11_class::Flush()
 141 {
 142         fflush(infile);
 143 }
 144
 145 void File11_class::Rewind()
 146 {
 147         fseek(infile, 0, SEEK_SET);
 148 }
 149
 150 void File11_class::WriteRowOfPhanPixels(REAL* test_aray, INTEGER nelem)
 151 {
 152
 153         for (int k = 0; k < nelem; k++)
 154         {
 155                 if ((k) % 7 == 0)
 156                 {
 157                         fprintf(infile, "\n");
 158                 }
 159
 160                 fprintf(infile, " %25.20f", test_aray[k]);
 161
 162                 if (test_aray[k] != 0)
 163                 {
 164                         fflush(infile);
 165                 }
 166         }
 167 }
 168
 169 void File11_class::WriteEndOfPhan()
 170 {
 171         fprintf(infile, "\n");
 172 }
 173
 174 void File11_class::WriteProj(REAL theta, REAL* nlo, INTEGER usrays)
 175 {
 176         REAL thedeg = theta * 180 / Consts.pi;
 177
 178
 179         fprintf(infile, "\n %25.20f %25.20f", theta, thedeg);
 180
 181         for (int nr = 0; nr < usrays; nr++)
 182         {
 183                 if ((nr) % 3 == 0)
 184                 {
 185                         fprintf(infile, "\n");
 186                 }
 187                 fprintf(infile, " %25.20f", nlo[nr]);
 188         }
 189 }
 190
 191 void File11_class::WriteEndOfProj()
 192 {
 193         fprintf(infile, "\n");
 194 }
 195
 196 void File11_class::ReadPicture(REAL* Pict, INTEGER nelem)
 197 {
 198         REAL* jlo = Pict;
 199         signed CHAR ch;  // added by Lajos, Nov 18, 2004
 200
 201 //#pragma message( "to fix!" )
 202
 203         for(int i = 0; i < nelem; i++){
 204             REAL val;
 205
 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");
 210                    exit(444);
 211                   }
 212                 *jlo++ = val;
 213                }
 214            }
 215
 216         while (((ch = fgetc(infile)) != '\n') && (ch != EOF))
 217         {  // added check for EOF. Lajos, Nov 18, 2004
 218                 ;
 219         }
 220
 221 }
 222
 223 void File11_class::wrhdr(INTEGER word)
 224 {
 225
 226         static const INTEGER keywd[40] =
 227         {
 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(' ', ' ', ' ', ' ')
 249         };
 250
 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
 257                         "    " };
 258
 259         static const CHAR kindStr[5][5] =
 260         { "elip", "rect", "tria", "segm", "sect" };
 261
 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
 283
 284         // DATA SPACE FOR OBJECTS IN PHANTOM AT TOP OF BLANK COMMON
 285
 286         /*
 287          real object(17,1)
 288          equivalence (object, work(2) )
 289          */
 290
 291         BOOLEAN perfct;
 292         BOOLEAN noisfl;
 293         BOOLEAN quanfl;
 294         ///FILE* file;
 295
 296 // VARIABLES INDEXING THE KEYWORDS IN THE ARRAY KEYWD
 297
 298         /*
 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,
 302          3   pixel, equ, pnch
 303
 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/
 311          */
 312
 313         /*
 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/
 317          */
 318
 319         //static const REAL one = 1.0;
 320         int i, j;
 321
 322         INTEGER itype;
 323
 324         quanfl = NoisePar.quanin > 0;
 325         noisfl = quanfl || NoisePar.ultnfl || NoisePar.addnfl;
 326         perfct = !(NoisePar.sctnfl || noisfl);
 327
 328         // SET THE OUTPUT FILE TO FILE11 OR PUNCH
 329         /*
 330          file = file11;
 331
 332          if(word == keywd[KWI_pnch]) {
 333          // ALSO SET VALUES PROPERLY FOR UNAVAILABLE QUANTITIES
 334          ///    file = punch;
 335          Creacm.nobj = 0;
 336          Spctrm.nergy = 1;
 337          Spctrm.energy[0] = -1;
 338          }
 339          */
 340
 341         // WRITE HEADER INFORMATION ON FILE11 FOR PHANTOM & RAYSUMS
 342         // PUNCH OUT THE HEADER INFORMATION FOR THE PHANTOM AND PICTURES
 343
 344         fprintf(infile, "%s", Creacm.name);
 345
 346         // SPECTRUM
 347         if (Spctrm.nergy == 1)
 348         {
 349                 fprintf(infile, "\n%s    %s %4i", keywdStr[KWI_spec],
 350                                 keywdStr[KWI_mono], Spctrm.energy[0]);
 351         }
 352
 353         if (Spctrm.nergy > 1)
 354         {
 355
 356                 fprintf(infile, "\n%s    %4s %4i", keywdStr[KWI_spec],
 357                                 keywdStr[KWI_poly], Spctrm.nergy);
 358                 fprintf(infile, "\n");
 359
 360                 for (i = 0; i < Spctrm.nergy; i++)
 361                 {
 362                         fprintf(infile, "%5i %5i ", Spctrm.energy[i], Creacm.percnt[i]);
 363                 }
 364         }
 365
 366         // OBJECTS
 367
 368         fprintf(infile, "\n%s", keywdStr[KWI_obj]);
 369
 370         if (Creacm.nobj != 0)
 371         {
 372
 373                 for (i = 0; i < Creacm.nobj; i++)
 374                 {
 375                         itype = (int) objects[i].type;
 376
 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]);
 380
 381                         if (Spctrm.nergy > 1)
 382                         {
 383
 384                                 fprintf(infile, "\ndens ");
 385                                 for (j = 1; j < Spctrm.nergy; j++)
 386                                 {
 387                                         fprintf(infile, " %15.10f", objects[i].den1[j]);
 388                                 }
 389                         }
 390                 }
 391         }
 392
 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);
 397
 398         // SKIP PHANTOM AVERAGING AND OTHER INFO SUCH AS NELEM FOR RDPROJ
 399
 400         if (word != keywd[KWI_ray])
 401         {
 402                 if (word == keywd[KWI_phan])
 403                 {
 404
 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);
 409                 }
 410
 411                 return;
 412         }
 413
 414         // HEADER INFO FOR RAYSUMS
 415         // GEOMETRY
 416
 417         fprintf(infile, "\n%4s    %4s %4i\n", keywdStr[KWI_ray], keywdStr[KWI_aver],
 418                         GeoPar.nave2);
 419         for (i = 0; i < GeoPar.nave2; i++)
 420         {
 421                 fprintf(infile, " %4i", GeoPar.naper[i]);
 422         }
 423
 424         fprintf(infile, "\n%4s", keywdStr[KWI_geom]);
 425
 426         if (GeoPar.lino)  /// bug #135. Lajos, 08/02/2005
 427         {
 428                 fprintf(infile, "\n%s", keywdStr[KWI_lino]);
 429         }
 430         else
 431         {
 432                 if (GeoPar.par)
 433                 {
 434                         fprintf(infile, "\n%4s", int2str(Creacm.word1));
 435                         fprintf(infile, " %4s", int2str(Creacm.word2));
 436                         fprintf(infile, " %4s", int2str(Creacm.word3));
 437                 }
 438
 439                 if (GeoPar.div)
 440                 {
 441                         CHAR wrd2[5];
 442
 443                         strcpy(wrd2, int2str(Creacm.word2));
 444                         fprintf(infile,
 445                                         "\n%4s    %4s    source at%10.4f     det dist%10.4f",
 446                                         int2str(Creacm.word1), wrd2, GeoPar.radius, GeoPar.stod);
 447                 }
 448
 449                 fprintf(infile, "\n%s    user    %5i    spacing    %10.4f",
 450                                 keywdStr[KWI_ray], GeoPar.usrays, GeoPar.pinc);
 451
 452
 453                 fprintf(infile, "\n%s    %5i\n", keywdStr[KWI_angl], GeoPar.prjnum);
 454                 for (i = 0; i < GeoPar.prjnum; i++)
 455                 {
 456                         if (((i % 7) == 0) && (i != 0))
 457                         {
 458                                 fprintf(infile, "\n");
 459                         }
 460                         fprintf(infile, " %9.4f", Creacm.pang[i]);
 461                 }
 462
 463         }
 464
 465         // MEASUREMENT STATISTICS
 466
 467         if (perfct)
 468         {
 469                 fprintf(infile, "\n%s    %s", keywdStr[KWI_meas], keywdStr[KWI_perf]);
 470         }
 471
 472
 473         else
 474         {
 475                 fprintf(infile, "\n%s    %s", keywdStr[KWI_meas], keywdStr[KWI_nois]);
 476         }
 477
 478         if (quanfl)
 479         {
 480                 fprintf(infile, "\n%s    %14.4f    %14.4f    cali %2i",
 481                                 keywdStr[KWI_quan], NoisePar.quanmn, NoisePar.quancm,
 482                                 NoisePar.quanin);
 483         }
 484
 485         if (NoisePar.sctnfl)
 486         {
 487                 fprintf(infile, "\n%s    noise peak %20.10f    width %20.10f",
 488                                 keywdStr[KWI_scat], NoisePar.sctnpk, NoisePar.sctnwd);
 489         }
 490
 491         if (NoisePar.ultnfl)
 492         {
 493                 fprintf(infile, "\n%s    noise mean %20.10f    std dev %20.10f",
 494                                 keywdStr[KWI_mult], NoisePar.ultnmn, NoisePar.ultnsd);
 495         }
 496
 497         if (NoisePar.addnfl)
 498         {
 499                 fprintf(infile, "\n%s    noise mean %20.10f    std dev %20.10f",
 500                                 keywdStr[KWI_add], NoisePar.addnmn, NoisePar.addnsd);
 501         }
 502
 503         if (!perfct && noisfl)
 504         {
 505                 fprintf(infile, "\n%s    %4i", keywdStr[KWI_seed], Creacm.iseed);
 506         }
 507
 508         fprintf(infile, "\n%4s    ", keywdStr[KWI_back]);
 509         for (i = 0; i < Spctrm.nergy; i++)
 510         {
 511                 fprintf(infile, " %7.4f", Spctrm.backgr[i]);
 512         }
 513
 514 }
 515