src/DIGFileSnark/DIGFileSnarkProj.cpp

   1 #include <cstdio>
   2 #include <string>
   3
   4 //#include "DIGFileSnarkProj.h"
   5 #include <DIGFileSnark/DIGFileSnarkProj.h>
   6 \r
   7 #define DIGFILESNARKPROJ_DEBUG_LEVEL 0\r
   8
   9 static char DIGFileSnarkPrjEmptyStr[1] = "";
  10
  11 // Geometry
  12
  13 // divergent
  14 ProjGeometry::ProjGeometry(double pRadius, double pSrcDetDistance, DetectorTypeEnum pDetectorType)
  15 {\r
  16 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
  17   printf("ProjGeometry::ProjGeometry\n");\r
  18 #endif\r
  19
  20         Radius = pRadius;
  21         SrcDetDistance = pSrcDetDistance;
  22         DetectorType = pDetectorType;
  23         GeometryType = GT_DIVERGENT;
  24 }
  25
  26 // parallel
  27 ProjGeometry::ProjGeometry(DetectorTypeEnum pDetectorType, ProjTypeEnum pProjType)
  28 {\r
  29 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
  30   printf("ProjGeometry::ProjGeometry\n");\r
  31 #endif\r
  32
  33         DetectorType = pDetectorType;
  34         ProjType = pProjType;
  35         GeometryType = GT_PARALLEL;
  36 }
  37
  38
  39 // Energy Spectrum
  40
  41 ProjSpectrum::ProjSpectrum(unsigned int eNum)
  42 {\r
  43 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
  44   printf("ProjSpectrum::ProjSpectrum\n");\r
  45 #endif\r
  46
  47         NumberOfEnergies = eNum;
  48         energy = new int[eNum];
  49         ratio = new double[eNum];
  50   background = new double[eNum];
  51 }
  52
  53
  54 int ProjSpectrum::SetEnergy(int i, int e, double r, double bg)
  55 {\r
  56 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
  57   printf("ProjSpectrum::SetEnergy\n");\r
  58 #endif\r
  59
  60   energy[i] = e;
  61   ratio[i] = r;
  62   background[i] = bg;
  63   return 0;
  64 }
  65
  66 /*
  67 int PrjSpectrum_class::SetSpectrum(int* e, double* r, double* bg)
  68 {
  69   for(unsigned int i = 0; i < NumberOfEnergies; i++) {
  70     energy[i] = e[i];
  71     ratio[i] = r[i];
  72     background[i] = bg[i];
  73   }
  74
  75   return 0;
  76 }
  77
  78 // Energy Spectrum
  79
  80 int PrjSpectrum_class::GetNumberOfEnergies(unsigned int* NOE)
  81 {
  82   *NOE = NumberOfEnergies;
  83
  84         return 0;
  85 }
  86
  87
  88 int PrjSpectrum_class::GetEnergy(int index, int* e, double* r, double* bg)
  89 {
  90   *e = energy[index];
  91         *r = ratio[index];
  92         *bg = background[index];
  93
  94         return 0;
  95 }
  96 */
  97
  98 // Noise
  99
 100 ProjNoise::ProjNoise()
 101 {\r
 102 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 103   printf("ProjNoise::ProjNoise\n");\r
 104 #endif\r
 105
 106         IsQuan = IsScat = IsAdd = IsMult = false;
 107 };
 108
 109 // Quantum
 110
 111 void ProjNoise::SetQuantum(double pMean, double pCalibration, int pGantry)
 112 {\r
 113 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 114   printf("ProjNoise::SetQuantum\n");\r
 115 #endif\r
 116
 117         IsQuan = true;
 118
 119         QuantumMean = pMean;
 120         QuantumCalibration = pCalibration;
 121         QuantumGantry = pGantry;
 122 }
 123
 124
 125 // Scatter
 126
 127 void ProjNoise::SetScatter(double pPeak, double pWidth)
 128 {\r
 129 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 130   printf("ProjNoise::SetScatter\n");\r
 131 #endif\r
 132
 133         IsScat = true;
 134
 135         ScatterPeak = pPeak;
 136         ScatterWidth = pWidth;
 137 }
 138
 139
 140 // Additive
 141
 142 void ProjNoise::SetAdditive(double pMean, double pStdDev)
 143 {\r
 144 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 145   printf("ProjNoise::SetAdditive\n");\r
 146 #endif\r
 147
 148         IsAdd = true;
 149
 150         AdditiveMean = pMean;
 151         AdditiveStdDev = pStdDev;
 152 }
 153
 154
 155 // Multiplicative
 156
 157 void ProjNoise::SetMultiplicative(double pMean, double pStdDev)
 158 {\r
 159 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 160   printf("ProjNoise::SetMultiplicative\n");\r
 161 #endif\r
 162
 163         IsMult = true;
 164
 165         MultiplicativeMean = pMean;
 166         MultiplicativeStdDev = pStdDev;
 167 }
 168
 169
 170 // Noise
 171
 172 // Quantum
 173 /*
 174 int PrjNoise_class::GetQuantum(double& m, double& cm, double& g)
 175 {
 176   if(!IsQuan) {
 177     return -1;
 178   }
 179
 180         m = Quantum_Mean;
 181         cm = Quantum_Calibration;
 182         g = Quantum_Gantry;
 183
 184   return 0;
 185 }
 186
 187
 188 // Scatter
 189
 190 int PrjNoise_class::GetScatter(double& p, double& w)
 191 {
 192   if(!IsScat) {
 193     return -1;
 194   }
 195
 196         p = Scatter_Peak;
 197         w = Scatter_Width;
 198
 199   return 0;
 200 }
 201
 202
 203 // Additive
 204
 205 int PrjNoise_class::GetAdditive(double& m, double& sd)
 206 {
 207   if(!IsAdd) {
 208     return -1;
 209   }
 210
 211         m = Additive_Mean;
 212         sd = Additive_Std_Dev;
 213
 214   return 0;
 215 }
 216
 217
 218 // Multiplicative
 219
 220 int PrjNoise_class::GetMultiplicative(double& m, double& sd)
 221 {
 222   if(!IsMult) {
 223     return -1;
 224   }
 225
 226         m = Multiplicative_Mean;
 227         sd = Multiplicative_Std_Dev;
 228
 229   return 0;
 230 }
 231 */
 232
 233 const char* DIGFileSnarkProj::TypeStr = "SNARK05 prjfil";
 234 const char* DIGFileSnarkProj::SchemaStr = "DIGFileSnarkProj.xsd";
 235
 236 bool DIGFileSnarkProj::XMLStringsInitialized = false;
 237
 238 const XMLCh* DIGFileSnarkProj::AppHeaderXMLStr;
 239 const XMLCh* DIGFileSnarkProj::GeomHeaderXMLStr;
 240 const XMLCh* DIGFileSnarkProj::DivergentXMLStr;
 241 const XMLCh* DIGFileSnarkProj::RadiusXMLStr;
 242 const XMLCh* DIGFileSnarkProj::SourceDetXMLStr;
 243 const XMLCh* DIGFileSnarkProj::DetTypeXMLStr;
 244 const XMLCh* DIGFileSnarkProj::ParallelXMLStr;
 245 const XMLCh* DIGFileSnarkProj::ProjTypeXMLStr;
 246 const XMLCh* DIGFileSnarkProj::NoiseHeaderXMLStr;
 247 const XMLCh* DIGFileSnarkProj::QuantumXMLStr;
 248 const XMLCh* DIGFileSnarkProj::MeanXMLStr;
 249 const XMLCh* DIGFileSnarkProj::CalibMeasXMLStr;
 250 const XMLCh* DIGFileSnarkProj::GantryXMLStr;
 251 const XMLCh* DIGFileSnarkProj::ScatterXMLStr;
 252 const XMLCh* DIGFileSnarkProj::PeakXMLStr;
 253 const XMLCh* DIGFileSnarkProj::WidthXMLStr;
 254 const XMLCh* DIGFileSnarkProj::AdditiveXMLStr;
 255 const XMLCh* DIGFileSnarkProj::StdDevXMLStr;
 256 const XMLCh* DIGFileSnarkProj::MultipXMLStr;
 257 const XMLCh* DIGFileSnarkProj::SpectrumHeaderXMLStr;
 258 const XMLCh* DIGFileSnarkProj::NumOfEnerXMLStr;
 259 const XMLCh* DIGFileSnarkProj::EnergyLevelXMLStr;
 260 const XMLCh* DIGFileSnarkProj::EnergyXMLStr;
 261 const XMLCh* DIGFileSnarkProj::RatioXMLStr;
 262 const XMLCh* DIGFileSnarkProj::BackgroundXMLStr;
 263
 264 DIGFileSnarkProj::DIGFileSnarkProj()
 265 {
 266
 267   // XML string initialization
 268   if(!XMLStringsInitialized) {
 269
 270           AppHeaderXMLStr      = XMLString::transcode("application_header");
 271           GeomHeaderXMLStr     = XMLString::transcode("geometric_header");
 272           DivergentXMLStr      = XMLString::transcode("divergent");
 273           RadiusXMLStr         = XMLString::transcode("radius");
 274           SourceDetXMLStr      = XMLString::transcode("source_detector_distance");
 275           DetTypeXMLStr        = XMLString::transcode("detector_type");
 276           ParallelXMLStr       = XMLString::transcode("parallel");
 277           ProjTypeXMLStr       = XMLString::transcode("projection_type");
 278           NoiseHeaderXMLStr    = XMLString::transcode("noise_header");
 279           QuantumXMLStr        = XMLString::transcode("quantum");
 280           MeanXMLStr           = XMLString::transcode("mean");
 281           CalibMeasXMLStr      = XMLString::transcode("calibration_measurement");
 282           GantryXMLStr         = XMLString::transcode("gantry_arrangement");
 283           ScatterXMLStr        = XMLString::transcode("scatter");
 284           PeakXMLStr           = XMLString::transcode("peak");
 285           WidthXMLStr          = XMLString::transcode("width");
 286           AdditiveXMLStr       = XMLString::transcode("additive");
 287           StdDevXMLStr         = XMLString::transcode("standard_deviation");
 288           MultipXMLStr         = XMLString::transcode("multiplicative");
 289           SpectrumHeaderXMLStr = XMLString::transcode("spectrum_header");
 290           NumOfEnerXMLStr      = XMLString::transcode("number_of_energies");
 291           EnergyLevelXMLStr    = XMLString::transcode("energy_level");
 292           EnergyXMLStr         = XMLString::transcode("energy");
 293           RatioXMLStr          = XMLString::transcode("ratio");
 294           BackgroundXMLStr     = XMLString::transcode("background");
 295   }
 296
 297   Spectrum = NULL;
 298   Geometry = NULL;
 299 }
 300
 301
 302 DIGFileSnarkProj::~DIGFileSnarkProj()
 303 {\r
 304 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 305   fprintf(stderr, "DIGFileSnarkProj:~DIGFileSnarkProj\n");\r
 306 #endif\r
 307
 308   if(Geometry != NULL) {
 309           delete Geometry;
 310   }
 311
 312   if(Spectrum != NULL) {
 313           delete Spectrum;
 314   }
 315 }
 316
 317 ////////////////////////////////////////////////////////////////
 318 ////////////////////////////////////////////////////////////////
 319
 320
 321 ///////////////////////////////////////////////
 322 // main header reading
 323 ///////////////////////////////////////////////
 324
 325 int DIGFileSnarkProj::GetMainHeader(ProjFileMH* MainHeader)
 326 {\r
 327 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 328   fprintf(stderr, "DIGFileSnarkProj:GetMainHeader\n");\r
 329 #endif\r
 330
 331   // read title
 332   const char* Title;
 333 \r
 334 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 335   printf("Getting Main Header Title\n");\r
 336 #endif\r
 337
 338   if(GetTitle(&Title) != 0) {
 339     printf("Error reading Main Header Title\n");
 340     return -1;
 341   }
 342
 343   //printf("Main Header Title: %s\n", Title);\r
 344 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 345   printf("Seting Main Header Title\n");\r
 346 #endif\r
 347
 348   if(MainHeader->Title.Set(Title) != 0) {
 349     printf("Error seting Main Header Title\n");
 350     return -1;
 351   }\r
 352 \r
 353 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 354   printf("Seting Main Header Dimensions\n");\r
 355 #endif
 356
 357   if(GetDimensions(&(MainHeader->Dimensions)) != 0) {
 358     printf("Error reading Main Header Dimensions\n");
 359     return -1;
 360   }
 361
 362   //printf("Main Header Dimensions: %d\n", MainHeader->Dimensions);
 363 \r
 364 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 365   printf("Seting Main Header Sampling\n");\r
 366 #endif
 367
 368   if(GetSampling(&(MainHeader->Sampling)) != 0) {
 369     printf("Error reading Main Header Sampling\n");
 370     return -1;
 371   }
 372
 373   //printf("Main Header Sampling: %f\n", MainHeader->Sampling);
 374 \r
 375 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 376   printf("Getting Main Header comments\n");\r
 377 #endif\r
 378
 379   const char* Comment;
 380
 381   if(GetComment(&Comment) != 0) {
 382     printf("Error reading Main Header Comment\n");
 383     return -1;
 384   }
 385 \r
 386 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 387   printf("Seting Main Header Comments\n");\r
 388 #endif\r
 389
 390   //printf("Main Header Comment: %s\n", Comment);
 391
 392   if(MainHeader->Comment.Set(Comment) != 0) {
 393     printf("Error Seting Main Header Comment\n");
 394     return -1;
 395   }
 396
 397   return 0;
 398 }
 399
 400 int DIGFileSnarkProj::GetAppHeader(ProjFileAH* pAppHeader)
 401 {\r
 402 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 403   fprintf(stderr, "DIGFileSnarkProj:GetAppHeader\n");\r
 404 #endif\r
 405
 406   // application header comment
 407   const char* Comment;
 408 \r
 409 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 410   printf("Getting Application Header Comments\n");\r
 411 #endif\r
 412
 413   if(GetAppComment(&Comment) != 0) {
 414     printf("Error reading Application Header Comment\n");
 415     return -1;
 416   }
 417 \r
 418 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 419   printf("Setting Application Header Comments\n");\r
 420 #endif\r
 421
 422   //printf("Application Header Comment: %s\n", Comment);
 423
 424   if(pAppHeader->Comment.Set(Comment) != 0) {
 425     printf("Error Seting Application Header Comment\n");
 426     return -1;
 427   }
 428 \r
 429 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
 430   printf("Getting Application Header Geometry Type\n");\r
 431 #endif\r
 432
 433   // geometry
 434   if(GetGeometryType(&(pAppHeader->GeometryType)) != 0) {
 435     printf("Error reading Application Header Geometry Type\n");
 436     return -1;
 437   }
 438
 439   //printf("Application Header Geometry Type: %s\n", (pAppHeader->GeometryType == GT_DIVERGENT) ? "Divergent": "Parallel");
 440
 441   if(pAppHeader->GeometryType == GT_DIVERGENT) {  // divergent
 442
 443     if(GetGeometry(&(pAppHeader->Radius), &(pAppHeader->SrcDetDistance), &(pAppHeader->DetectorType)) != 0) {
 444       printf("Error reading Application Header Geometry\n");
 445       return -1;
 446     }
 447
 448     //printf("Application Header Geometry Radius: %f\n", pAppHeader->Radius);
 449     //printf("Application Header Geometry SrcDetDistance: %f\n", pAppHeader->SrcDetDistance);
 450
 451   }
 452   else { // parallel
 453
 454     if(GetGeometry(&(pAppHeader->DetectorType), &(pAppHeader->ProjType)) != 0) {
 455       printf("Error reading Application Header Geometry\n");
 456       return -1;
 457     }
 458
 459     //printf("Application Header Geometry Projection Type: %s\n", (pAppHeader->ProjType == PT_STRIP) ? "Strip": "Line");
 460
 461   }
 462
 463   /*
 464   switch(pAppHeader->DetectorType) {
 465
 466   case DT_TANGENT:
 467     printf("Application Header Geometry Detector Type: TANGENT\n");
 468     break;
 469
 470   case DT_ARC:
 471     printf("Application Header Geometry Detector Type: ARC\n");
 472     break;
 473
 474   case DT_UNIFORM:
 475     printf("Application Header Geometry Detector Type: UNIFORM\n");
 476     break;
 477
 478   case DT_VARIABLE:
 479     printf("Application Header Geometry Detector Type: VARIABLE\n");
 480     break;
 481   }
 482   */
 483
 484   // spectrum
 485
 486   if(GetNoOfEnergies(&(pAppHeader->NoOfEnergies)) != 0) {
 487     printf("Error reading Application Header Spectrum Number Of Energies\n");
 488     return -1;
 489   }
 490
 491   //printf("Application Header Spectrum Number Of Energies: %u\n", pAppHeader->NoOfEnergies);
 492
 493   pAppHeader->Energy = new int[pAppHeader->NoOfEnergies];
 494   pAppHeader->Ratio = new double[pAppHeader->NoOfEnergies];
 495   pAppHeader->Background = new double[pAppHeader->NoOfEnergies];
 496
 497   for(unsigned int k = 0; k < pAppHeader->NoOfEnergies; k++) {
 498     if(GetSpectrum(k, &(pAppHeader->Energy[k]), &(pAppHeader->Ratio[k]), &(pAppHeader->Background[k])) != 0) {
 499       printf("Error reading Application Header Spectrum Energy\n");
 500       return -1;
 501     }
 502
 503     //printf("Application Header Spectrum Energy: %u %i %f %f\n", k, pAppHeader->Energy[k], pAppHeader->Ratio[k], pAppHeader->Background[k]);
 504
 505   }
 506
 507   // noise
 508   if(HasNoiseQuantum(&(pAppHeader->QuantumFlag)) != 0) {
 509     printf("Error reading Application Header Noise Quantum Flag\n");
 510     return -1;
 511   }
 512
 513   if(pAppHeader->QuantumFlag) {
 514
 515     if(GetNoiseQuantum(&(pAppHeader->QuantumM), &(pAppHeader->QuantumCM), &(pAppHeader->QuantumGA)) != 0) {
 516       printf("Error reading Application Header Noise Quantum Flag\n");
 517       return -1;
 518     }
 519
 520     //printf("Application Header Noise Quantum: %f %f %i\n", pAppHeader->QuantumM, pAppHeader->QuantumCM, pAppHeader->QuantumGA);
 521   }
 522
 523   if(HasNoiseScatter(&(pAppHeader->ScatterFlag)) != 0) {
 524     printf("Error reading Application Header Noise Scatter Flag\n");
 525     return -1;
 526   }
 527
 528   if(pAppHeader->ScatterFlag) {
 529
 530     if(GetNoiseScatter(&(pAppHeader->ScatterP), &(pAppHeader->ScatterW)) != 0) {
 531       printf("Error reading Application Header Noise Scatter\n");
 532       return -1;
 533     }
 534
 535     //printf("Application Header Noise Scatter: %f %f\n", pAppHeader->ScatterP, pAppHeader->ScatterW);
 536   }
 537
 538   if(HasNoiseAdditive(&(pAppHeader->AdditiveFlag)) != 0) {
 539     printf("Error reading Application Header Noise Additive Flag\n");
 540     return -1;
 541   }
 542
 543   if(pAppHeader->AdditiveFlag) {
 544     if(GetNoiseAdditive(&(pAppHeader->AdditiveM), &(pAppHeader->AdditiveSD)) != 0) {
 545       printf("Error reading Application Header Noise Additive\n");
 546       return -1;
 547     }
 548
 549     //printf("Application Header Noise Additive: %f %f\n", pAppHeader->AdditiveM, pAppHeader->AdditiveSD);
 550   }
 551
 552   if(HasNoiseMultiplicative(&(pAppHeader->MultiplicativeFlag)) != 0) {
 553     printf("Error reading Application Header Noise Multiplicative Flag\n");
 554     return -1;
 555   }
 556
 557   if(pAppHeader->MultiplicativeFlag) {
 558
 559     if(GetNoiseMultiplicative(&(pAppHeader->MultiplicativeM), &(pAppHeader->MultiplicativeSD)) != 0) {
 560       printf("Error reading Application Header Noise Multiplicative\n");
 561       return -1;
 562     }
 563
 564     //printf("Application Header Noise Multiplicative: %f %f\n", pAppHeader->MultiplicativeM, pAppHeader->MultiplicativeSD);
 565   }
 566
 567   return 0;
 568 }
 569
 570 ///////////////////////////////////////////////
 571 // Projection header reading
 572 ///////////////////////////////////////////////
 573
 574 int DIGFileSnarkProj::GetProjHeader(ProjFilePrjH* pProjHeader)
 575 {\r
 576 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 577   fprintf(stderr, "DIGFileSnarkProj:GetProjHeader\n");\r
 578 #endif\r
 579
 580   // read projection number
 581
 582   if(GetProjNo(&(pProjHeader->ProjNo)) != 0) {
 583     printf("Error reading Projection Number\n");
 584     return -1;
 585   }
 586
 587   //printf("Projection Number: %u\n", pProjHeader->ProjNo);
 588
 589   // read angle
 590
 591   if(GetProjAngle(&(pProjHeader->Angle)) != 0) {
 592     printf("Error reading Projection Angle\n");
 593     return -1;
 594   }
 595
 596   //printf("Array Set Parameters: %f\n", pProjHeader->Angle);
 597
 598   // read comment
 599   const char* Comment;
 600
 601   if(GetProjComment(&Comment) != 0) {
 602     printf("Error reading Projection Comment\n");
 603     return -1;
 604   }
 605
 606   //printf("Projection Comment: %s\n", Comment);
 607
 608   if(pProjHeader->Comment.Set(Comment) != 0) {
 609     printf("Error Seting Projection Comment\n");
 610     return -1;
 611   }
 612
 613   return 0;
 614 }
 615
 616 ///////////////////////////////////////////////////////////////////////////////
 617 ///////////////////////////////////////////////////////////////////////////////
 618 // Writing
 619 ///////////////////////////////////////////////////////////////////////////////
 620 ///////////////////////////////////////////////////////////////////////////////
 621
 622 ///////////////////////////////////////////////////////////////////////////////
 623 // Open for writing
 624 ///////////////////////////////////////////////////////////////////////////////
 625
 626 int DIGFileSnarkProj::Open(const char* pFileName, ProjFileMH* pMainHeader)
 627 {\r
 628 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 629   fprintf(stderr, "DIGFileSnarkProj:Open\n");\r
 630 #endif\r
 631
 632   if(
 633     Open(
 634       pFileName,                  // file name
 635       pMainHeader->Title.Get(),   // title
 636       pMainHeader->Dimensions,    // dimensions
 637       pMainHeader->Sampling,      // sampling
 638       pMainHeader->Comment.Get()  // comment
 639     ) != 0
 640   )
 641   {
 642     printf("Error opening file\n");
 643     return -1;
 644   }
 645   return 0;
 646 }
 647
 648 // open for writing
 649 int DIGFileSnarkProj::Open(
 650   const char*   pFileName,     // File Name
 651   const char*   pTitle,        // Title
 652   unsigned int  pNumberOfRays, // USRAYS
 653   double        pSampling,     // PINC
 654   const char*   pComments      // Comments
 655 )
 656 {
 657   //const char*           Schema;
 658   DIGDimensions  dim;
 659   DIGSampling    sampX;
 660   DIGSampling    samp = { 0.0, 0.0, 0.0 };
 661
 662   dim.x = pNumberOfRays;
 663   dim.y = 1;
 664   dim.z = 1;
 665
 666   sampX.x = pSampling;
 667   sampX.y = 0.0;
 668   sampX.z = 0.0;
 669
 670   Noise.IsQuan = false;
 671   Noise.IsAdd = false;
 672   Noise.IsScat = false;
 673   Noise.IsMult = false;
 674
 675   return DIGFile::Open(
 676     pFileName,            // File Name
 677     SchemaStr,            // Schema
 678     pTitle,               // Title
 679     TypeStr,              // Type,
 680     1,                    // Chanels
 681     DIGValueType_REAL,    // ValueType
 682     DIGDataType_DOUBLE,   // DataType
 683     DIGDataFormat_BINARY, // DataFormat
 684     DIGGrid_SC,           // Grid
 685     DIGBasis_VORONOI,     // Basis
 686     DIGUnit_UNSPECIFIED,  // Unit
 687     &dim,                 // Dimensions
 688     &sampX,               // SamplingX
 689     &samp,                // SamplingY
 690     &samp,                // SamplingZ
 691     pComments,            // Comments
 692     ""                    // Other Unit
 693   );
 694 };
 695
 696
 697 ///////////////////////////////////////////////
 698 // application header writing
 699 ///////////////////////////////////////////////
 700
 701 int DIGFileSnarkProj::SetAppHeader(ProjFileAH* pAppHeader)
 702 {\r
 703 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 704   fprintf(stderr, "DIGFileSnarkProj:SetAppHeader\n");\r
 705 #endif\r
 706
 707   const char* Comment;
 708
 709   if(pAppHeader->Comment.Get(&Comment)) {\r
 710 \r
 711 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 712   printf("** Error: Getting Application Header Comments\n");\r
 713 #endif
 714     return -1;
 715   }
 716
 717   SetAppComment(Comment);
 718
 719   switch(pAppHeader->GeometryType) {
 720   case GT_DIVERGENT:
 721     //pOut->SetGeometry(double r, double d, DetectorTypeEnum t);
 722     if(SetGeometry(pAppHeader->Radius, pAppHeader->SrcDetDistance, pAppHeader->DetectorType) != 0) {\r
 723 \r
 724 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 725       printf("** Error: Setting Application Header Geometry\n");\r
 726 #endif
 727       return -1;
 728     }
 729     break;
 730
 731   case GT_PARALLEL:
 732     if(SetGeometry(pAppHeader->DetectorType, pAppHeader->ProjType) != 0) {\r
 733 \r
 734 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 735       printf("** Error: Setting Application Header Geometry\n");\r
 736 #endif
 737       return -1;
 738     }
 739     break;
 740
 741   default:\r
 742 \r
 743 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 744       printf("** Error: Setting Application Header - Unknown Geometry\n");\r
 745 #endif
 746     ;// Error
 747   }
 748
 749   if(SetSpectrum(pAppHeader->NoOfEnergies, pAppHeader->Energy, pAppHeader->Ratio, pAppHeader->Background) != 0) {\r
 750 \r
 751 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 752     printf("** Error: Setting Application Header Spectrum\n");\r
 753 #endif
 754     return -1;
 755   }
 756
 757   if(pAppHeader->QuantumFlag) {
 758           if(SetNoiseQuantum(pAppHeader->QuantumM, pAppHeader->QuantumCM, pAppHeader->QuantumGA) != 0) {\r
 759 \r
 760 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 761     printf("** Error: Setting Application Header Noise - Quantum\n");\r
 762 #endif
 763       return -1;
 764     }
 765   }
 766
 767   if(pAppHeader->ScatterFlag) {
 768           if(SetNoiseScatter(pAppHeader->ScatterP, pAppHeader->ScatterW) != 0) {\r
 769 \r
 770 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 771     printf("** Error: Setting Application Header Noise - Scatter\n");\r
 772 #endif
 773       return -1;
 774     }
 775   }
 776
 777   if(pAppHeader->AdditiveFlag) {
 778           if(SetNoiseAdditive(pAppHeader->AdditiveM, pAppHeader->AdditiveSD) != 0) {\r
 779 \r
 780 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 781     printf("** Error: Setting Application Header Noise - Additive\n");\r
 782 #endif
 783       return -1;
 784     }
 785   }
 786
 787   if(pAppHeader->MultiplicativeFlag) {
 788           if(SetNoiseMultiplicative(pAppHeader->MultiplicativeM, pAppHeader->MultiplicativeSD) != 0) {\r
 789 \r
 790 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 791     printf("** Error: Setting Application Header Noise - Multiplicative\n");\r
 792 #endif
 793       return -1;
 794     }
 795   }
 796
 797   return 0;
 798 }
 799
 800 ///////////////////////////////////////////////
 801 // projection appending
 802 ///////////////////////////////////////////////
 803
 804 int DIGFileSnarkProj::AppendProj(ProjFilePrjH* pProjHeader, const double* pData)
 805 {\r
 806 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 807   fprintf(stderr, "DIGFileSnarkProj:AppendProj\n");\r
 808 #endif\r
 809
 810   return AppendProj(
 811     pProjHeader->ProjNo,
 812     pProjHeader->Angle,
 813     pProjHeader->Comment.Get(),
 814     pData
 815   );
 816 }
 817
 818 int DIGFileSnarkProj::AppendProj(\r
 819   unsigned int  pProjNo,
 820   double        pAngle,
 821   const char*   pComment,
 822   const double* pData
 823 )
 824 {\r
 825 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 826   fprintf(stderr, "DIGFileSnarkProj:AppendProj\n");\r
 827 #endif\r
 828
 829   int ret;
 830   char ProjNoStr[8];   // Projection Number as a String
 831   char AngleStr[256];  // angle as a String
 832
 833         AngleStr[0] = 0;
 834         sprintf(AngleStr, "%.20f", pAngle);
 835
 836         ProjNoStr[0] = 0;
 837         sprintf(ProjNoStr, "%d", pProjNo);
 838
 839   if((ret = AppendArraySet("projection", ProjNoStr, AngleStr, pComment)) != 0) {
 840     return ret;
 841   }
 842
 843   return AppendArray(0, "", pData);
 844 }
 845
 846 ////////////////////////////////////////////////////////////////
 847 ////////////////////////////////////////////////////////////////
 848
 849 // open for reading
 850
 851 int DIGFileSnarkProj::Open(const char* pFileName)\r
 852 {
 853 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 854   fprintf(stderr, "DIGFileSnarkProj:Open\n");\r
 855 #endif\r
 856
 857   Noise.IsQuan = false;
 858   Noise.IsAdd = false;
 859   Noise.IsScat = false;
 860   Noise.IsMult = false;
 861
 862   if(DIGFile::Open(pFileName) != 0) {\r
 863 \r
 864 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 865     printf("** Error: Openning DIG file\n");\r
 866 #endif
 867     return -1; // error opening dig file
 868   }
 869
 870   if(strcmp(MainHeader.Type.Get(), TypeStr) != 0) {\r
 871 \r
 872 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 873     printf("** Error: Not a Snark05 prjfil file\n");\r
 874 #endif
 875     return -2; // not a DIGFileSnarkProj
 876   }
 877
 878   return 0;
 879 }
 880
 881 // comment
 882
 883 int DIGFileSnarkProj::GetAppComment(const char** pComment)
 884 {\r
 885 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 886   fprintf(stderr, "DIGFileSnarkProj:GetAppComment\n");\r
 887 #endif\r
 888
 889   return AppComment.Get(pComment);
 890 }
 891
 892 int DIGFileSnarkProj::GetGeometryType(GeometryTypeEnum* pGeometryType)
 893 {\r
 894 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 895   fprintf(stderr, "DIGFileSnarkProj:GetGeometryType\n");\r
 896 #endif\r
 897
 898   *pGeometryType = Geometry->GeometryType;
 899   return 0;
 900 }
 901
 902 // divergent
 903 int DIGFileSnarkProj::GetGeometry(double* pRadius, double* pSrcDetDistance, DetectorTypeEnum* pDetectorType)
 904 {\r
 905 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 906   fprintf(stderr, "DIGFileSnarkProj:GetGeometry\n");\r
 907 #endif\r
 908
 909   if(Geometry == NULL) {\r
 910 \r
 911 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 912     printf("** Error: Undefines Geometry\n");\r
 913 #endif
 914     return -1; // error
 915   }
 916
 917         *pRadius = Geometry->Radius;
 918         *pSrcDetDistance = Geometry->SrcDetDistance;
 919         *pDetectorType = Geometry->DetectorType;
 920
 921   return 0;
 922 }
 923
 924 // parallel
 925 int DIGFileSnarkProj::GetGeometry(DetectorTypeEnum* pDetectorType, ProjTypeEnum* pProjType)
 926 {\r
 927 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 928   fprintf(stderr, "DIGFileSnarkProj:GetGeometry\n");\r
 929 #endif\r
 930
 931   if(Geometry == NULL) {\r
 932 \r
 933 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 934     printf("** Error: Undefines Geometry\n");\r
 935 #endif
 936     return -1; // error
 937   }
 938
 939         *pDetectorType = Geometry->DetectorType;
 940         *pProjType = Geometry->ProjType;
 941
 942   return 0;
 943 }
 944
 945
 946 int DIGFileSnarkProj::GetNoOfEnergies(unsigned int* NumbOfEnergies)
 947 {\r
 948 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 949   fprintf(stderr, "DIGFileSnarkProj:GetNoOfEnergies\n");\r
 950 #endif\r
 951
 952   *NumbOfEnergies = Spectrum->NumberOfEnergies;
 953
 954   return 0;
 955 }
 956
 957 int DIGFileSnarkProj::GetSpectrum(unsigned int i, int* e, double* r, double* bg)
 958 {\r
 959 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 960   fprintf(stderr, "DIGFileSnarkProj:GetSpectrum\n");\r
 961 #endif\r
 962
 963   if(Spectrum == NULL) {\r
 964 \r
 965 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 966     printf("** Error: GetSpectrum - Undefines Spectrum\n");\r
 967 #endif
 968     return -1; // error
 969   }
 970
 971   if(i >= Spectrum->NumberOfEnergies) {\r
 972 \r
 973 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
 974     printf("** Error: GetSpectrum - invalid parameter\n");\r
 975 #endif
 976     return -2; // error
 977   }
 978
 979   *e = Spectrum->energy[i];
 980   *r = Spectrum->ratio[i];
 981   *bg = Spectrum->background[i];
 982
 983   return 0;
 984 }
 985
 986
 987 int DIGFileSnarkProj::HasNoiseQuantum(bool* a)
 988 {\r
 989 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
 990   fprintf(stderr, "DIGFileSnarkProj:HasNoiseQuantum\n");\r
 991 #endif\r
 992
 993   *a = Noise.IsQuan;
 994   return 0;
 995 }
 996
 997
 998 int DIGFileSnarkProj::GetNoiseQuantum(double* pMean, double* pCalibration, int* pGantry)
 999 {\r
1000 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1001   fprintf(stderr, "DIGFileSnarkProj:GetNoiseQuantum\n");\r
1002 #endif\r
1003
1004   if(Noise.IsQuan != true) {\r
1005 \r
1006 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1007     printf("** Error: GetNoiseQuantum - Quantum Noise Undefined\n");\r
1008 #endif
1009     return -1; // error
1010   }
1011
1012         *pMean = Noise.QuantumMean;
1013   *pCalibration = Noise.QuantumCalibration;
1014   *pGantry = Noise.QuantumGantry;
1015
1016   return 0;
1017 }
1018
1019
1020 // Scatter
1021
1022 int DIGFileSnarkProj::HasNoiseScatter(bool* a)
1023 {\r
1024 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1025   fprintf(stderr, "DIGFileSnarkProj:HasNoiseScatter\n");\r
1026 #endif\r
1027
1028   *a = Noise.IsScat;
1029   return 0;
1030 }
1031
1032 int DIGFileSnarkProj::GetNoiseScatter(double* pPeak, double* pWidth)
1033 {\r
1034 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1035   fprintf(stderr, "DIGFileSnarkProj:GetNoiseScatter\n");\r
1036 #endif\r
1037
1038   if(Noise.IsScat != true) {\r
1039 \r
1040 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1041     printf("** Error: GetNoiseScatter - Scatter Noise Undefined\n");\r
1042 #endif
1043     return -1; // error
1044   }
1045
1046         *pPeak = Noise.ScatterPeak;
1047   *pWidth = Noise.ScatterWidth;
1048
1049   return 0;
1050 }
1051
1052
1053 // Additive
1054
1055 int DIGFileSnarkProj::HasNoiseAdditive(bool* a)
1056 {\r
1057 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1058   fprintf(stderr, "DIGFileSnarkProj:HasNoiseAdditive\n");\r
1059 #endif\r
1060
1061   *a = Noise.IsAdd;
1062   return 0;
1063 }
1064
1065
1066 int DIGFileSnarkProj::GetNoiseAdditive(double* pMean, double* pStdDev)
1067 {\r
1068 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1069   fprintf(stderr, "DIGFileSnarkProj:GetNoiseAdditive\n");\r
1070 #endif\r
1071
1072   if(Noise.IsAdd != true) {\r
1073 \r
1074 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1075     printf("** Error: GetNoiseAdditive - Additive Noise Undefined\n");\r
1076 #endif
1077     return -1; // error
1078   }
1079
1080         *pMean = Noise.AdditiveMean;
1081   *pStdDev = Noise.AdditiveStdDev;
1082
1083   return 0;
1084 }
1085
1086
1087 // Multiplicative
1088
1089 int DIGFileSnarkProj::HasNoiseMultiplicative(bool* a)
1090 {\r
1091 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1092   fprintf(stderr, "DIGFileSnarkProj:HasNoiseMultiplicative\n");\r
1093 #endif\r
1094
1095   *a = Noise.IsMult;
1096   return 0;
1097 }
1098
1099
1100 int DIGFileSnarkProj::GetNoiseMultiplicative(double* pMean, double* pStdDev)
1101 {\r
1102 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1103   fprintf(stderr, "DIGFileSnarkProj:GetNoiseMultiplicative\n");\r
1104 #endif\r
1105
1106   if(Noise.IsMult != true) {\r
1107 \r
1108 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1109     printf("** Error: GetNoiseMultiplicative - Multiplicative Noise Undefined\n");\r
1110 #endif
1111     return -1; // error
1112   }
1113
1114         *pMean = Noise.MultiplicativeMean;
1115   *pStdDev = Noise.MultiplicativeStdDev;
1116
1117   return 0;
1118 }
1119
1120
1121
1122 // comment
1123
1124 int DIGFileSnarkProj::SetAppComment(const char* pComment)
1125 {\r
1126 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1127   fprintf(stderr, "DIGFileSnarkProj:SetAppComment\n");\r
1128 #endif\r
1129
1130   return AppComment.Set(pComment);
1131 }
1132
1133 // divergent
1134 int DIGFileSnarkProj::SetGeometry(double pRadius, double pSrcDetDistance, DetectorTypeEnum pDetectorType)
1135 {\r
1136 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1137   fprintf(stderr, "DIGFileSnarkProj:SetGeometry\n");\r
1138 #endif\r
1139
1140   if(Geometry != NULL) {\r
1141 \r
1142 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1143     printf("** Error: SetGeometry - Geometry Aleready Defined\n");\r
1144 #endif
1145     return -1; // already set
1146   }
1147
1148         Geometry = new ProjGeometry(pRadius, pSrcDetDistance, pDetectorType);
1149
1150   return 0;
1151 }
1152
1153 // parallel
1154 int DIGFileSnarkProj::SetGeometry(DetectorTypeEnum pDetectorType, ProjTypeEnum pProjType)
1155 {\r
1156 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1157   fprintf(stderr, "DIGFileSnarkProj:SetGeometry\n");\r
1158 #endif\r
1159
1160   if(Geometry != NULL) {\r
1161 \r
1162 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1163     printf("** Error: SetGeometry - Geometry Aleready Defined\n");\r
1164 #endif
1165     return -1; // already set
1166   }
1167
1168         Geometry = new ProjGeometry(pDetectorType, pProjType);
1169
1170   return 0;
1171 }
1172
1173 int DIGFileSnarkProj::SetSpectrum(unsigned int eNum, int* e, double* r, double* bg)
1174 {\r
1175 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1176   fprintf(stderr, "DIGFileSnarkProj:SetSpectrum\n");\r
1177 #endif\r
1178
1179   if(Spectrum != NULL) {\r
1180 \r
1181 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1182     printf("** Error: SetSpectrum - Spectrum Aleready Defined\n");\r
1183 #endif
1184     return -1; // already set
1185   }
1186
1187         Spectrum = new ProjSpectrum(eNum);
1188
1189   for(unsigned int i = 0; i < eNum; i++) {
1190     Spectrum->energy[i] = e[i];
1191     Spectrum->ratio[i] = r[i];
1192     Spectrum->background[i] = bg[i];
1193   }
1194
1195   return 0;
1196 }
1197
1198
1199 int DIGFileSnarkProj::SetNoiseQuantum(double pMean, double pCalibration, int pGantry)
1200 {\r
1201 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1202   fprintf(stderr, "DIGFileSnarkProj:SetNoiseQuantum\n");\r
1203 #endif\r
1204
1205   if(Noise.IsQuan) {\r
1206 \r
1207 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1208     printf("** Error: SetNoiseQuantum - Noise Quantum Aleready Defined\n");\r
1209 #endif
1210     return -1; // already set
1211   }
1212
1213         Noise.SetQuantum(pMean, pCalibration, pGantry);
1214
1215   return 0;
1216 }
1217
1218
1219 // Scatter
1220
1221 int DIGFileSnarkProj::SetNoiseScatter(double pPeak, double pWidth)
1222 {\r
1223 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1224   fprintf(stderr, "DIGFileSnarkProj:SetNoiseScatter\n");\r
1225 #endif\r
1226
1227   if(Noise.IsScat) {\r
1228 \r
1229 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1230     printf("** Error: SetNoiseScatter - Noise Scatter Aleready Defined\n");\r
1231 #endif
1232     return -1; // already set
1233   }
1234
1235         Noise.SetScatter(pPeak, pWidth);
1236
1237   return 0;
1238 }
1239
1240
1241 // Additive
1242
1243 int DIGFileSnarkProj::SetNoiseAdditive(double pMean, double pStdDev)
1244 {\r
1245 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1246   fprintf(stderr, "DIGFileSnarkProj:SetNoiseAdditive\n");\r
1247 #endif\r
1248
1249   if(Noise.IsAdd) {\r
1250 \r
1251 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1252     printf("** Error: SetNoiseAdditive - Noise Additive Aleready Defined\n");\r
1253 #endif
1254     return -1; // already set
1255   }
1256
1257         Noise.SetAdditive(pMean, pStdDev);
1258
1259   return 0;
1260 }
1261
1262
1263 // Multiplicative
1264
1265 int DIGFileSnarkProj::SetNoiseMultiplicative(double pMean, double pStdDev)
1266 {\r
1267 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1268   fprintf(stderr, "DIGFileSnarkProj:SetNoiseMultiplicative\n");\r
1269 #endif\r
1270
1271   if(Noise.IsMult) {\r
1272 \r
1273 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 0\r
1274     printf("** Error: SetNoiseMultiplicative - Noise Multiplicative Aleready Defined\n");\r
1275 #endif
1276     return -1; // already set
1277   }
1278
1279         Noise.SetMultiplicative(pMean, pStdDev);
1280
1281   return 0;
1282 }
1283
1284
1285 /////////////////////////////////////////////////////////////////////
1286 //
1287 // Internal Methods
1288 //
1289 /////////////////////////////////////////////////////////////////////
1290
1291 // NOTE: error handling is not included
1292 int DIGFileSnarkProj::ParseApplicationHeader()
1293 {\r
1294 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1295   fprintf(stderr, "DIGFileSnarkProj:ParseApplicationHeader\n");\r
1296 #endif\r
1297 \r
1298 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
1299     printf("Parsing Application Header\n");\r
1300 #endif\r
1301 \r
1302 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
1303     printf("Getting Comments\n");\r
1304 #endif\r
1305
1306   // get Comment Texts
1307   char* tCom = DIGFileSnarkPrjEmptyStr;\r
1308
1309   DOMNodeList* CommentsList = ApplicationHeaderElement->getElementsByTagName(CommentsXMLStr);
1310
1311   int tNumC = CommentsList->getLength();
1312   if(tNumC > 0) {
1313
1314     DOMElement* CommentsElement = (DOMElement*) CommentsList->item(0);
1315
1316     if(GetElementText(CommentsElement, &tCom) != 0) {
1317       ;// allow empty comments
1318     }
1319   }
1320
1321   AppComment.Set(tCom);
1322
1323 \r
1324 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
1325     printf("Getting Geometric Header\n");\r
1326 #endif\r
1327
1328   // get geometric header
1329
1330   DOMNodeList* GeomHeaderList = ApplicationHeaderElement->getElementsByTagName(GeomHeaderXMLStr);
1331   DOMElement* GeomHeaderElement = (DOMElement*) GeomHeaderList->item(0);
1332
1333   // get divergent
1334
1335   DOMNodeList* DivergentList = GeomHeaderElement->getElementsByTagName(DivergentXMLStr);
1336         if(DivergentList->getLength() > 0) {
1337                 DOMElement* DivergentElement = (DOMElement*) DivergentList->item(0);
1338
1339                 // get radius, source detector distance, detector type
1340
1341                 double tRadius = 0;
1342                 GetDoubleAttributeValue(DivergentElement, RadiusXMLStr, &tRadius);
1343
1344                 double tSrcDetDist = 0;
1345                 GetDoubleAttributeValue(DivergentElement, SourceDetXMLStr, &tSrcDetDist);
1346
1347                 DOMNodeList* DetectList = DivergentElement->getElementsByTagName(DetTypeXMLStr);
1348                 DOMElement* DetectElement = (DOMElement*) DetectList->item(0);
1349
1350                 char* tstring;
1351                 DetectorTypeEnum tDet;
1352
1353                 GetElementText(DetectElement, &tstring);
1354
1355     if(strcmp(tstring, "tangent") == 0) {
1356       tDet = DT_TANGENT;
1357     }
1358
1359     if(strcmp(tstring, "arc") == 0) {
1360       tDet = DT_ARC;
1361     }
1362
1363                 Geometry = new ProjGeometry(tRadius, tSrcDetDist, tDet);
1364         }
1365         else {
1366                 // get parallel tag
1367
1368                 DOMNodeList* ParallelList = GeomHeaderElement->getElementsByTagName(ParallelXMLStr);
1369                 DOMElement* ParallelElement = (DOMElement*) ParallelList->item(0);
1370
1371                 // get detector type
1372
1373                 DOMNodeList* DetectList = ParallelElement->getElementsByTagName(DetTypeXMLStr);
1374                 DOMElement* DetectElement = (DOMElement*) DetectList->item(0);
1375
1376                 char *tstring;
1377                 DetectorTypeEnum tDet;
1378
1379                 GetElementText(DetectElement, &tstring);
1380
1381     if(strcmp(tstring, "uniform") == 0) {
1382       tDet = DT_UNIFORM;
1383     }
1384
1385     if(strcmp(tstring, "variable") == 0) {
1386       tDet = DT_VARIABLE;
1387     }
1388
1389                 // get projection type
1390
1391                 DOMNodeList* ProjectList = ParallelElement->getElementsByTagName(ProjTypeXMLStr);
1392                 DOMElement* ProjectElement = (DOMElement*) ProjectList->item(0);
1393
1394                 ProjTypeEnum tProj;
1395
1396                 GetElementText(ProjectElement, &tstring);
1397
1398     if(strcmp(tstring, "line") == 0) {
1399       tProj = PT_LINE;
1400     }
1401     else if(strcmp(tstring, "strip") == 0) {
1402       tProj = PT_STRIP;
1403     }
1404
1405
1406                 Geometry = new ProjGeometry(tDet, tProj);
1407         }
1408 \r
1409 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
1410     printf("Getting Noise Header\n");\r
1411 #endif\r
1412
1413         // get noise header
1414
1415         DOMNodeList* NoiseHeaderList = ApplicationHeaderElement->getElementsByTagName(NoiseHeaderXMLStr);
1416         DOMElement* NoiseHeaderElement = (DOMElement*) NoiseHeaderList->item(0);
1417
1418         // get quantum
1419
1420         DOMNodeList* QuantumList = NoiseHeaderElement->getElementsByTagName(QuantumXMLStr);
1421         if(QuantumList->getLength() > 0) {
1422                 DOMElement* QuantumElement = (DOMElement*) QuantumList->item(0);
1423
1424                 // get mean, calibration measurement, gantry arrangement
1425
1426                 double tMean = 0, tCal = 0;
1427                 int tGan = 0;
1428
1429                 GetDoubleAttributeValue(QuantumElement, MeanXMLStr, &tMean);
1430                 GetDoubleAttributeValue(QuantumElement, CalibMeasXMLStr, &tCal);
1431                 GetIntAttributeValue(QuantumElement, GantryXMLStr, &tGan);
1432
1433                 Noise.SetQuantum(tMean, tCal, tGan);
1434         }
1435
1436         // get scatter
1437
1438         DOMNodeList* ScatterList = NoiseHeaderElement->getElementsByTagName(ScatterXMLStr);
1439         if(ScatterList->getLength() > 0) {
1440                 DOMElement* ScatterElement = (DOMElement*) ScatterList->item(0);
1441
1442                 // get peak, width
1443
1444                 double tPeak = 0, tWidth = 0;
1445
1446                 GetDoubleAttributeValue(ScatterElement, PeakXMLStr, &tPeak);
1447                 GetDoubleAttributeValue(ScatterElement, WidthXMLStr, &tWidth);
1448
1449                 Noise.SetScatter(tPeak, tWidth);
1450         }
1451
1452         // get additive
1453
1454         DOMNodeList* AdditList = NoiseHeaderElement->getElementsByTagName(AdditiveXMLStr);
1455         if(AdditList->getLength() > 0) {
1456                 DOMElement* AdditElement = (DOMElement*) AdditList->item(0);
1457
1458                 // get mean, calibration measurement, gantry arrangement
1459
1460                 double tMean = 0, tStd = 0;
1461
1462                 GetDoubleAttributeValue(AdditElement, MeanXMLStr, &tMean);
1463                 GetDoubleAttributeValue(AdditElement, StdDevXMLStr, &tStd);
1464
1465                 Noise.SetAdditive(tMean, tStd);
1466         }
1467
1468         // get multiplicative
1469
1470         DOMNodeList* MultipList = NoiseHeaderElement->getElementsByTagName(MultipXMLStr);
1471         if(MultipList->getLength() > 0) {
1472                 DOMElement* MultipElement = (DOMElement*) MultipList->item(0);
1473
1474                 // get mean, calibration measurement, gantry arrangement
1475
1476                 double tMean = 0, tStd = 0;
1477
1478                 GetDoubleAttributeValue(MultipElement, MeanXMLStr, &tMean);
1479                 GetDoubleAttributeValue(MultipElement, StdDevXMLStr, &tStd);
1480
1481                 Noise.SetMultiplicative(tMean, tStd);
1482         }
1483         \r
1484 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 1\r
1485     printf("Getting Specrum Header\n");\r
1486 #endif\r
1487
1488         // get spectrum header
1489
1490         DOMNodeList* SpectHeaderList = ApplicationHeaderElement->getElementsByTagName(SpectrumHeaderXMLStr);
1491         DOMElement* SpectHeaderElement = (DOMElement*) SpectHeaderList->item(0);
1492
1493         // get number of energies
1494
1495         unsigned int tNumE = 0;
1496
1497         GetUintAttributeValue(SpectHeaderElement, NumOfEnerXMLStr, &tNumE);
1498         //SetEnergyNumber(tNumE);
1499
1500   Spectrum = new ProjSpectrum(tNumE);
1501
1502         // get energy level triplets: energy, ratio, background
1503
1504         DOMNodeList* EnergyList = SpectHeaderElement->getElementsByTagName(EnergyLevelXMLStr);
1505
1506         // NOTE: Error if EnergyList->getLength() != tNumE
1507
1508         int tEner = 0;
1509         double tRatio = 0, tBack = 0;
1510         DOMElement* EnergyElement;
1511
1512         for(unsigned int i = 0; i < EnergyList->getLength(); i++) {
1513                 EnergyElement = (DOMElement*) EnergyList->item(i);
1514
1515                 GetIntAttributeValue(EnergyElement, EnergyXMLStr, &tEner);
1516                 GetDoubleAttributeValue(EnergyElement, RatioXMLStr, &tRatio);
1517                 GetDoubleAttributeValue(EnergyElement, BackgroundXMLStr, &tBack);
1518
1519                 Spectrum->SetEnergy(i, tEner, tRatio, tBack);
1520         }
1521
1522         return 0;
1523 }
1524
1525
1526 int DIGFileSnarkProj::CreateApplicationHeader()
1527 {\r
1528 #if DIGFILESNARKPROJ_DEBUG_LEVEL > 3\r
1529   fprintf(stderr, "DIGFileSnarkProj:CreateApplicationHeader\n");\r
1530 #endif\r
1531
1532   char valStr[20];
1533
1534   XMLCh* tmpXMLStr;
1535
1536   // create main header node
1537   ApplicationHeaderElement = Document->createElement(AppHeaderXMLStr);
1538
1539   // attributes
1540
1541     // there is no attributes
1542
1543   // children
1544
1545
1546   // define comments
1547   const char* Comment;
1548
1549   AppComment.Get(&Comment);
1550
1551         if(Comment != NULL) {
1552
1553     // create comments node
1554     DOMElement* Comm = Document->createElement(CommentsXMLStr);
1555
1556     // create comments text node
1557     tmpXMLStr = XMLString::transcode(Comment);
1558     DOMText* CommentText = Document->createTextNode(tmpXMLStr);
1559     //delete tmpXMLStr;
1560     Comm->appendChild(CommentText);
1561
1562     // append comments to appication header
1563     ApplicationHeaderElement->appendChild(Comm);
1564   }
1565
1566
1567   // geometry
1568
1569   // create geometry node
1570   DOMElement* Geom = Document->createElement(GeomHeaderXMLStr);
1571
1572         if(Geometry->GeometryType == GT_DIVERGENT) {
1573
1574     DOMElement* Diver = Document->createElement(DivergentXMLStr);
1575
1576     // define radius
1577     sprintf(valStr, "%f", Geometry->Radius);
1578     tmpXMLStr = XMLString::transcode(valStr);
1579     Diver->setAttribute(RadiusXMLStr, tmpXMLStr);
1580     //delete tmpXMLStr;
1581
1582     // define source detector distanse
1583     sprintf(valStr, "%f", Geometry->SrcDetDistance);
1584     tmpXMLStr = XMLString::transcode(valStr);
1585     Diver->setAttribute(SourceDetXMLStr, tmpXMLStr);
1586     //delete tmpXMLStr;
1587
1588     // define detector type
1589     DOMElement* DetT = Document->createElement(DetTypeXMLStr);
1590
1591     // create detector type text node
1592     if(Geometry->DetectorType == DT_TANGENT) {
1593       tmpXMLStr = XMLString::transcode("tangent");
1594     }
1595     else {
1596       tmpXMLStr = XMLString::transcode("arc");
1597     }
1598
1599     DOMText* DetTText = Document->createTextNode(tmpXMLStr);
1600     //delete tmpXMLStr;
1601     DetT->appendChild(DetTText);
1602
1603     // append detector type to divergent
1604     Diver->appendChild(DetT);
1605
1606     // append divergent to geometry
1607     Geom->appendChild(Diver);
1608   }
1609   else {
1610
1611     DOMElement* Para = Document->createElement(ParallelXMLStr);
1612
1613     // define detector type
1614     DOMElement* DetT = Document->createElement(DetTypeXMLStr);
1615
1616     if(Geometry->DetectorType == DT_UNIFORM) {
1617       tmpXMLStr = XMLString::transcode("uniform");
1618     }
1619     else {
1620       tmpXMLStr = XMLString::transcode("variable");
1621     }
1622
1623     DOMText* DetTText = Document->createTextNode(tmpXMLStr);
1624     //delete tmpXMLStr;
1625     DetT->appendChild(DetTText);
1626
1627     // append detector type to parallel
1628     Para->appendChild(DetT);
1629
1630
1631     // define projection type
1632     DOMElement* ProjT = Document->createElement(ProjTypeXMLStr);
1633
1634     if(Geometry->ProjType == PT_STRIP) {
1635       tmpXMLStr = XMLString::transcode("strip");
1636     }
1637     else {
1638       tmpXMLStr = XMLString::transcode("line");
1639     }
1640
1641     DOMText* ProjTText = Document->createTextNode(tmpXMLStr);
1642     //delete tmpXMLStr;
1643     ProjT->appendChild(ProjTText);
1644
1645     // append projection type to parallel
1646     Para->appendChild(ProjT);
1647
1648     // append parallel to geometry
1649     Geom->appendChild(Para);
1650   }
1651
1652   // append geometry to application header
1653   ApplicationHeaderElement->appendChild(Geom);
1654
1655
1656         // noise
1657
1658   // create noise node
1659   DOMElement* Nois = Document->createElement(NoiseHeaderXMLStr);
1660
1661         if(Noise.IsQuan) {
1662     DOMElement* Quan = Document->createElement(QuantumXMLStr);
1663
1664     // define mean
1665     sprintf(valStr, "%f", Noise.QuantumMean);
1666     tmpXMLStr = XMLString::transcode(valStr);
1667     Quan->setAttribute(MeanXMLStr, tmpXMLStr);
1668     //delete tmpXMLStr;
1669
1670     // define calibration
1671     sprintf(valStr, "%f", Noise.QuantumCalibration);
1672     tmpXMLStr = XMLString::transcode(valStr);
1673     Quan->setAttribute(CalibMeasXMLStr, tmpXMLStr);
1674     //delete tmpXMLStr;
1675
1676     // define gantry
1677     sprintf(valStr, "%i", Noise.QuantumGantry);
1678     tmpXMLStr = XMLString::transcode(valStr);
1679     Quan->setAttribute(GantryXMLStr, tmpXMLStr);
1680     //delete tmpXMLStr;
1681
1682     // append quantum to noise
1683     Nois->appendChild(Quan);
1684         }
1685
1686
1687   if(Noise.IsScat) {
1688     DOMElement* Scat = Document->createElement(ScatterXMLStr);
1689
1690     // define peak
1691     sprintf(valStr, "%f", Noise.ScatterPeak);
1692     tmpXMLStr = XMLString::transcode(valStr);
1693     Scat->setAttribute(PeakXMLStr, tmpXMLStr);
1694     //delete tmpXMLStr;
1695
1696     // define width
1697     sprintf(valStr, "%f", Noise.ScatterWidth);
1698     tmpXMLStr = XMLString::transcode(valStr);
1699     Scat->setAttribute(WidthXMLStr, tmpXMLStr);
1700     //delete tmpXMLStr;
1701
1702     // append scatter to noise
1703     Nois->appendChild(Scat);
1704         }
1705
1706
1707         if(Noise.IsAdd) {
1708     DOMElement* Add = Document->createElement(AdditiveXMLStr);
1709
1710     // define mean
1711     sprintf(valStr, "%f", Noise.AdditiveMean);
1712     tmpXMLStr = XMLString::transcode(valStr);
1713     Add->setAttribute(MeanXMLStr, tmpXMLStr);
1714     //delete tmpXMLStr;
1715
1716     // define standard deviation
1717     sprintf(valStr, "%f", Noise.AdditiveStdDev);
1718     tmpXMLStr = XMLString::transcode(valStr);
1719     Add->setAttribute(StdDevXMLStr, tmpXMLStr);
1720     //delete tmpXMLStr;
1721
1722     // append additive to noise
1723     Nois->appendChild(Add);
1724         }
1725
1726
1727         if(Noise.IsMult) {
1728     DOMElement* Mul = Document->createElement(MultipXMLStr);
1729
1730     // define mean
1731     sprintf(valStr, "%f", Noise.MultiplicativeMean);
1732     tmpXMLStr = XMLString::transcode(valStr);
1733     Mul->setAttribute(MeanXMLStr, tmpXMLStr);
1734     //delete tmpXMLStr;
1735
1736     // define standard deviation
1737     sprintf(valStr, "%f", Noise.MultiplicativeStdDev);
1738     tmpXMLStr = XMLString::transcode(valStr);
1739     Mul->setAttribute(StdDevXMLStr, tmpXMLStr);
1740     //delete tmpXMLStr;
1741
1742     // append multiplicative to noise
1743     Nois->appendChild(Mul);
1744         }
1745
1746   // append noise to application header
1747   ApplicationHeaderElement->appendChild(Nois);
1748
1749
1750         // spectrum
1751
1752   if(Spectrum == NULL) {
1753     return -2; // undefined spectrum
1754   }
1755
1756   // create spectrum node
1757   DOMElement* Spec = Document->createElement(SpectrumHeaderXMLStr);
1758
1759   // define number of energies
1760   tmpXMLStr = new XMLCh[5];
1761   XMLString::binToText(Spectrum->NumberOfEnergies, tmpXMLStr, 4, 10);
1762   Spec->setAttribute(NumOfEnerXMLStr, tmpXMLStr);
1763   delete tmpXMLStr;
1764
1765         for(unsigned int i = 0; i < Spectrum->NumberOfEnergies; i++) {
1766
1767     // create energy level node
1768     DOMElement* EnerL = Document->createElement(EnergyLevelXMLStr);
1769
1770     // define energy
1771     tmpXMLStr = new XMLCh[5];
1772     XMLString::binToText(Spectrum->energy[i], tmpXMLStr, 4, 10);
1773     EnerL->setAttribute(EnergyXMLStr, tmpXMLStr);
1774     delete tmpXMLStr;
1775
1776     // define ratio
1777     sprintf(valStr, "%f", Spectrum->ratio[i]);
1778     tmpXMLStr = XMLString::transcode(valStr);
1779     EnerL->setAttribute(RatioXMLStr, tmpXMLStr);
1780     //delete tmpXMLStr;
1781
1782     // define background
1783     sprintf(valStr, "%f", Spectrum->background[i]);
1784     tmpXMLStr = XMLString::transcode(valStr);
1785     EnerL->setAttribute(BackgroundXMLStr, tmpXMLStr);
1786     //delete tmpXMLStr;
1787
1788     // append energy level to spectrum
1789     Spec->appendChild(EnerL);
1790         }
1791
1792
1793   // append spectrum to application header
1794   ApplicationHeaderElement->appendChild(Spec);
1795
1796   // append application header to root
1797   RootElement->appendChild(ApplicationHeaderElement);
1798
1799         return 0;
1800 }