src/snark/projfile.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/projfile.cpp $
   5  $LastChangedRevision: 118 $
   6  $Date: 2014-07-09 14:22:29 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cstdlib>
  12 #include <cstring>
  13 #include <cmath>
  14
  15 #include "uiod.h"
  16
  17 #include "projfile.h"
  18
  19 SnarkPrjFile ProjFile;
  20
  21 SnarkPrjFile::SnarkPrjFile()
  22 {
  23         AnglesBuff = NULL;
  24 }
  25
  26 INTEGER SnarkPrjFile::Open(const char* FileName)
  27 {
  28
  29         if (DIGFileSnarkProj::Open(FileName) != 0)
  30         {
  31                 return -1;
  32         }
  33
  34         Angles = 1;
  35         return 0;
  36 }
  37
  38 INTEGER SnarkPrjFile::ReadHeader(char* name, Noise_class* Noise,
  39                 spctrm_class* Spct, anglst_class* Anglst)
  40 {
  41         const char* tname;
  42         double treal1, treal2;
  43         int tint;
  44         int i;
  45         unsigned int tuint;
  46
  47         //////////////
  48         // get Name //
  49         //////////////
  50
  51         GetTitle(&tname);
  52         strncpy(name, tname, 80);
  53         name[80] = 0;
  54
  55         /////////////////////////////
  56         // get GeoPar's attributes //
  57         /////////////////////////////
  58
  59         GetNoOfProjs(&tuint);
  60         NumberOfProjections = GeoPar.prjnum = tuint;
  61
  62         double tDIGSamp;
  63         GetSampling(&tDIGSamp);
  64         GeoPar.pinc = tDIGSamp;
  65
  66         unsigned int tDIGDim;
  67         GetNoOfRays(&tDIGDim);
  68         NumberOfRays = GeoPar.usrays = tDIGDim;
  69
  70         GeometryTypeEnum tGeomType;
  71         GetGeometryType(&tGeomType);
  72
  73         DetectorTypeEnum tDetType;
  74         ProjTypeEnum tProjType;
  75
  76         if (tGeomType == GT_DIVERGENT)
  77         {
  78                 GeoPar.div = true;
  79                 GeoPar.par = false;
  80                 GeoPar.strip = false; // bug 207 - swr - 3/3/06
  81                 GeoPar.line = true;
  82
  83                 GetGeometry(&treal1, &treal2, &tDetType);
  84
  85                 if (tDetType == DT_TANGENT)
  86                 {
  87                         GeoPar.tang = true;
  88                         GeoPar.arc = false;
  89                 }
  90                 else
  91                 {
  92                         GeoPar.tang = false;
  93                         GeoPar.arc = true;
  94                 };
  95
  96                 GeoPar.radius = treal1;
  97                 GeoPar.stod = treal2;
  98         }
  99         else
 100         {
 101                 GeoPar.div = false;
 102                 GeoPar.par = true;
 103
 104                 GetGeometry(&tDetType, &tProjType);
 105
 106                 if (tDetType == DT_UNIFORM)
 107                 {
 108                         GeoPar.uni = true;
 109                         GeoPar.vri = false;
 110                 }
 111                 else
 112                 {
 113                         GeoPar.uni = false;
 114                         GeoPar.vri = true;
 115                 };
 116
 117                 if (tProjType == PT_STRIP)
 118                 {
 119                         GeoPar.strip = true;
 120                         GeoPar.line = false;
 121                 }
 122                 else
 123                 {
 124                         GeoPar.strip = false;
 125                         GeoPar.line = true;
 126                 };
 127         };
 128
 129         ////////////////////////////
 130         // get Noise's attributes //
 131         ////////////////////////////
 132
 133         //bug221 start - swr - 03/03/07
 134         if (GetNoiseQuantum(&treal1, &treal2, &tint) == 0)
 135         {
 136                 Noise->setQuantumNoise(tint, treal1, treal2);
 137         }
 138         else
 139         {
 140                 Noise->setQuantumNoise(0, 1.0, 1.0);
 141         }
 142
 143         if (GetNoiseScatter(&treal1, &treal2) == 0)
 144         {
 145                 Noise->setScatter(treal1, treal2);
 146         };
 147
 148         if (GetNoiseAdditive(&treal1, &treal2) == 0)
 149         {
 150                 Noise->setAdditiveNoise(treal1, treal2);
 151         };
 152
 153         if (GetNoiseMultiplicative(&treal1, &treal2) == 0)
 154         {
 155                 Noise->setMultiplicativeNoise(treal1, treal2);
 156         };
 157         //bug221 end - swr - 03/03/07
 158
 159         ///////////////////////////
 160         // get Spct's attributes //
 161         ///////////////////////////
 162
 163         GetNoOfEnergies(&tuint);
 164         Spct->nergy = tuint;
 165
 166         for (i = 0; i < tuint; i++)
 167         {
 168                 GetSpectrum(i, &tint, &treal1, &treal2);
 169
 170                 Spct->energy[i] = tint;
 171                 Spct->engwt[i] = treal1;
 172                 Spct->backgr[i] = treal2;
 173         };
 174
 175         ///////////////////////////////////
 176         // calculate GeoPar's attributes //
 177         ///////////////////////////////////
 178
 179         GeoPar.snrays = GeoPar.numray(GeoPar.nelem, GeoPar.pixsiz);
 180
 181         if (GeoPar.usrays > GeoPar.snrays)
 182         {
 183                 GeoPar.nrays = GeoPar.usrays;
 184         }
 185         else
 186         { // bug194 - swr - 1/6/06
 187                 GeoPar.nrays = GeoPar.snrays;
 188         }
 189
 190
 191         GeoPar.midray = GeoPar.nrays / 2;
 192         GeoPar.fsnray = (GeoPar.nrays - GeoPar.snrays) / 2;
 193         GeoPar.lsnray = GeoPar.fsnray + GeoPar.snrays - 1;
 194         GeoPar.fusray = (GeoPar.nrays - GeoPar.usrays) / 2;
 195         GeoPar.lusray = GeoPar.fusray + GeoPar.usrays - 1;
 196
 197         //////////////////////////////
 198         // get Anglst's  attributes //
 199         //////////////////////////////
 200
 201         Anglst->Init(GeoPar.prjnum);
 202
 203         for (i = 0; i < GeoPar.prjnum; i++)
 204         {
 205
 206                 if (SelectProj(i) != 0)
 207                 {
 208                         return -1;
 209                 }
 210
 211                 if (GetProjAngle(&(Anglst->bth[i])) != 0)
 212                 {
 213                         return -1;
 214                 }
 215
 216                 Anglst->bsin[i] = (REAL) sin(Anglst->bth[i]);
 217                 Anglst->bcos[i] = (REAL) cos(Anglst->bth[i]);
 218         };
 219
 220         Anglst->pbase = new REAL[GeoPar.nrays];
 221         Anglst->incore = -1;
 222
 223         if (GeoPar.div)
 224         {
 225                 Anglst->genphi();
 226         }
 227
 228         return 0;
 229 }
 230
 231 INTEGER SnarkPrjFile::ReadProj(INTEGER np, REAL* data, unsigned INTEGER nsize)
 232 {
 233         unsigned int i;
 234         unsigned int tNum;
 235         INTEGER tSkip;
 236         REAL *tArray;
 237
 238         if (SelectProj(np) != 0)
 239         {
 240                 return -1;
 241         }
 242
 243         GetNoOfRays(&tNum);
 244         tArray = new REAL[tNum];
 245         GetProjData(tArray);
 246
 247         if (nsize <= NumberOfRays)
 248         {
 249
 250                 // NSIZE IS FEWER THAN OR EQUAL TO USRAYS
 251                 tSkip = (NumberOfRays - nsize) / 2;
 252                 for (i = 0; i < nsize; i++)
 253                 {
 254                         data[i] = tArray[i + tSkip];
 255                 }
 256         }
 257         else
 258         {
 259
 260                 // NSIZE IS GREATER THAN USRAYS
 261                 for (i = 0; i < nsize; i++)
 262                 {
 263                         data[i] = 0.0;
 264                 }
 265
 266                 tSkip = (nsize - NumberOfRays) / 2;
 267
 268                 for (i = 0; i < NumberOfRays; i++)
 269                 {
 270                         data[i + tSkip] = tArray[i];
 271                 }
 272         };
 273
 274         delete[] tArray;
 275
 276         return 0;
 277 }
 278
 279 INTEGER SnarkPrjFile::Open(const char* FileName, INTEGER usrays, INTEGER prjnum)
 280 {
 281         NumberOfRays = usrays;
 282         NumberOfProjections = prjnum;
 283
 284         AnglesBuff = new REAL[prjnum];
 285
 286         if (DIGFileSnarkProj::Open(FileName, "unknown", usrays, GeoPar.pinc,
 287                         "This projection file has been created by snark14") != 0)
 288                 return -1;
 289
 290         Angles = 0;
 291
 292         return 0;
 293 }
 294
 295 INTEGER SnarkPrjFile::WriteHeader(char* name, Geometry_class* Geo,
 296                 Noise_class* Noise, spctrm_class* Spct)
 297 {
 298         //////////////
 299         // set Name //
 300         //////////////
 301
 302         SetTitle(name);
 303
 304         /////////////////////////////
 305         // set GeoPar's attributes //
 306         /////////////////////////////
 307
 308         if (Geo->div)
 309         {
 310                 if (Geo->tang)
 311                 {
 312                         SetGeometry(Geo->radius, Geo->stod, DT_TANGENT);
 313                 }
 314                 else
 315                 {
 316                         SetGeometry(Geo->radius, Geo->stod, DT_ARC);
 317                 }
 318         };
 319
 320         if (Geo->par)
 321         {
 322                 if (Geo->uni)
 323                 {
 324                         if (Geo->strip)
 325                         {
 326                                 SetGeometry(DT_UNIFORM, PT_STRIP);
 327                         }
 328                         else
 329                         {
 330                                 SetGeometry(DT_UNIFORM, PT_LINE);
 331                         }
 332                 }
 333                 else
 334                 {
 335                         if (Geo->strip)
 336                         {
 337                                 SetGeometry(DT_VARIABLE, PT_STRIP);
 338                         }
 339                         else
 340                         {
 341                                 SetGeometry(DT_VARIABLE, PT_LINE);
 342                         }
 343                 };
 344         };
 345
 346         ////////////////////////////
 347         // set Noise's attributes //
 348         ////////////////////////////
 349
 350         if (Noise->quanin > 0)
 351         {
 352                 SetNoiseQuantum(Noise->quanmn, Noise->quancm, Noise->quanin);
 353         }
 354
 355         if (Noise->sctnfl == true)
 356         {
 357                 SetNoiseScatter(Noise->sctnpk, Noise->sctnwd);
 358         }
 359
 360         if (Noise->addnfl == true)
 361         {
 362                 SetNoiseAdditive(Noise->addnmn, Noise->addnsd);
 363         }
 364
 365         if (Noise->ultnfl == true)
 366         {
 367                 SetNoiseMultiplicative(Noise->ultnmn, Noise->ultnsd);
 368         }
 369
 370         ///////////////////////////
 371         // set Spct's attributes //
 372         ///////////////////////////
 373
 374         int i;
 375         double* treal1 = new double[Spct->nergy];
 376         double* treal2 = new double[Spct->nergy];
 377
 378         for (i = 0; i < Spct->nergy; i++)
 379         {
 380                 treal1[i] = Spct->engwt[i];
 381                 treal2[i] = Spct->backgr[i];
 382         };
 383
 384         SetSpectrum(Spct->nergy, Spct->energy, treal1, treal2);
 385
 386         delete[] treal1;
 387         delete[] treal2;
 388
 389         return 0;
 390 }
 391
 392 INTEGER SnarkPrjFile::WriteAngles(REAL* bth)
 393 {
 394
 395         unsigned int i;
 396
 397
 398         // if angles already stored report error
 399         if (Angles == 1)
 400         {
 401                 return -1;
 402         }
 403
 404         // store angles into buffer
 405         for (i = 0; i < NumberOfProjections; i++)
 406         {
 407                 AnglesBuff[i] = bth[i];
 408         };
 409
 410         Angles = 1;
 411         return 0;
 412 }
 413
 414 INTEGER SnarkPrjFile::WriteProj(INTEGER np, REAL* iprjec)
 415 {
 416
 417         if (Angles == 0)
 418         {
 419                 fprintf(output, "\n **** PRJFile error: projection can not "
 420                                 "be stored while angles are not given !!!\n");
 421                 return -1;
 422         };
 423
 424         if (AppendProj(np, AnglesBuff[np], "Projection created by snark14", iprjec)
 425                         != 0)
 426         {
 427                 return -1;
 428         }
 429
 430         return 0;
 431 }
 432
 433 INTEGER SnarkPrjFile::Close()
 434 {
 435
 436         if (AnglesBuff != NULL)
 437         {
 438                 delete[] AnglesBuff; // bug 92 - Lajos - 03/02/2005
 439                 AnglesBuff = NULL;
 440         }
 441
 442         // closes the file
 443
 444         if (DIGFileSnarkProj::Close() != 0)
 445         {
 446                 return 1;
 447         }
 448
 449         return 0;
 450
 451 }