src/snark/bh_correction.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/bh_correction.cpp $
   5  $LastChangedRevision: 118 $
   6  $Date: 2014-07-09 14:22:29 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  * Description:
  11  *
  12  * The class bh_correction is responsible for performing beam hardening
  13  * correction for polychromatic input data. Given the prjfil that was calculated
  14  * based on weighted average of attenuations at different energy levels and
  15  * corrections coefficients [1,2] it approximates monochromatic projection data
  16  * which is then used by the reconstruction algorithms.
  17  *
  18  *
  19  *
  20  * References:
  21  * [1] G. T. Herman, Image Reconstruction from Projections, 1980
  22  * [2] snark14 manual
  23  *
  24  */
  25
  26 #include <cstdlib>
  27 #include <cstdio>
  28 #include <cmath>
  29
  30 #include <unistd.h>
  31 #include <sys/types.h>
  32
  33 #include <fstream>
  34 #include <string.h>
  35
  36 #include "bh_correction.h"
  37
  38 using namespace std;
  39
  40 int bh_correction::readInputFilePoly()
  41 {
  42
  43         BOOLEAN eol, eol1, eol2;
  44         int i, j;   //loop counter variables
  45
  46         //read in the number of iterations
  47         noIter = InFile.getint(TRUE, &eol);
  48         if (eol == TRUE) //too few values provided
  49         {
  50                 fprintf(output, "\n **** ERROR: number of iterations missing\n"
  51                                 " **** program aborted\n\n");
  52                 exit(-1);
  53         }
  54         fprintf(output, "\n Number of iterations for beam hardening "
  55                         "correction set to r = %i\n\n", noIter);
  56
  57         //read in the degree of the polynomial that approximates monochromatic
  58         //projections given polychromatic projections
  59         noDegree = InFile.getint(FALSE, &eol);
  60         if (eol == TRUE) //too few values provided
  61         {
  62                 fprintf(output, "\n **** ERROR: the degree is missing\n"
  63                                 " **** program aborted\n\n");
  64                 exit(-1);
  65         }
  66         fprintf(output, "\n The degree of the polynomial for beam hardening "
  67                         "correction set to  %i\n\n", noDegree);
  68
  69         //read in the  coefficients of polynomial that approximates monochromatic
  70         //projections given polychromatic projections
  71         a = new REAL[noDegree + 1];
  72         a[0] = InFile.getnum(TRUE, &eol);
  73         if (eol == TRUE) //too few values provided
  74         {
  75                 fprintf(output, "\n **** ERROR: too few polynomial coefficients \n"
  76                                 " **** program aborted\n\n");
  77                 exit(-1);
  78         }
  79         for (i = 1; i <= noDegree; i++)
  80         {
  81                 a[i] = InFile.getnum(FALSE, &eol);
  82                 if (eol == TRUE) //too few values provided
  83                 {
  84                         fprintf(output, "\n **** ERROR: too few polynomial coefficients \n"
  85                                         " **** program aborted\n\n");
  86                         exit(-1);
  87                 }
  88         }
  89
  90         if (noIter != 0)
  91         {
  92                 //read in the number of energy levels
  93                 //this should be the same as Spctrm.nergy, otherwise exit
  94                 nerg = InFile.getint(TRUE, &eol);
  95                 if (eol == TRUE) //too few values provided
  96                 {
  97                         fprintf(output, "\n **** ERROR: the nerg is missing\n"
  98                                         " **** program aborted\n\n");
  99                         exit(-1);
 100                 }
 101
 102                 if (nerg != Spctrm.nergy) //these two should be the same
 103                 {
 104                         fprintf(output, "\n **** ERROR: the value of nerg does not agree\n"
 105                                         "with the nomber of energy levels written on file11\n"
 106                                         " **** program aborted\n\n");
 107                         exit(-1);
 108                 }
 109                 fprintf(output, "\n The number of energy levels is %i\n\n", nerg);
 110
 111                 //read in the number of points for piecewise linear conversion of
 112                 //phantom data from energy[0] to other energy levels
 113                 noPointsConv = InFile.getint(TRUE, &eol);
 114                 if (eol == TRUE) //too few values provided
 115                 {
 116                         fprintf(output, "\n **** ERROR: the noPointsConv is missing\n"
 117                                         " **** program aborted\n\n");
 118                         exit(-1);
 119                 }
 120                 fprintf(output,
 121                                 "\n The number of points for conversion of the polynomial "
 122                                                 "for beam hardening correction set to %i\n\n",
 123                                 noPointsConv);
 124
 125                 //increase noPointsConv by 1 since we need point (0,0)
 126                 noPointsConv++;
 127                 //read in the  points used for constant linear conversion of
 128                 //phantom data from energy[0] to other energy levels
 129                 for (j = 1; j < Spctrm.nergy; j++)
 130                 {
 131                         //skip j=0 since there is no conversion from the energy[0] to energy[0]
 132                         c[j] = new point[noPointsConv];
 133
 134                         c[j][0].x = 0.0; //these are always 0.0
 135                         c[j][0].y = 0.0; //because vacum attenuates nothing at all energy levels
 136
 137                         for (i = 1; i < noPointsConv; i++)
 138                         {
 139                                 c[j][i].x = InFile.getnum(TRUE, &eol1);
 140                                 c[j][i].y = InFile.getnum(FALSE, &eol2);
 141                                 if (eol1 == TRUE || eol2 == TRUE) //too few values provided
 142                                 {
 143                                         fprintf(output, "\n **** ERROR: too few point for "
 144                                                         "energy level %i \n **** program aborted\n\n", j);
 145                                         exit(-1);
 146                                 }
 147                         }
 148                 }
 149         }
 150
 151         return 0;
 152 }
 153
 154 ///////////////////////////////////////////////////////////////////////////////
 155
 156 int bh_correction::correction(SnarkPrjFile * BHC_PrjFile)
 157 {
 158         INTEGER i, j, k; //loop counters
 159         INTEGER r = 0;
 160         REAL x, y, u, den; //store elementary object descriptions
 161         char command[256]; //store command line to be used in system() call
 162         char filename[256]; //store file name created on the run
 163         char line[256]; //store temporarly a line of a file
 164
 165         if (noIter == 0)
 166         {
 167                 sprintf(command, "cp %s/q_p0 prjfil", tempDir);
 168                 if (system(command) != 0)
 169                 {
 170                         fprintf(output, "\n **** ERROR: could not execute  command"
 171                                         "%s\n  **** program aborted\n\n", command);
 172                         exit(-1);
 173                 }
 174                 return 0;
 175         }
 176
 177         //an input file for reconstruction part of iteration,
 178         //uses existing prjfil to generate recfil
 179         ofstream bhcInputFromPrjfil;
 180         //an array of input files to be created based on the recfil created by
 181         //reconstruction step during an iterative run
 182         //[0] will be used for monochromatic input file
 183         //[1] will be used for polychromatic input file
 184         ofstream bhcInputFiles[2];
 185
 186         //a copy of the original input file for the main snark14 run
 187         //it is needed here for information about projection generation, etc.
 188         ifstream inputFile;
 189
 190         //a recfil from the itermediate reconstruction during the iterative process
 191         DIGFileSnarkRec bhcRecfil;
 192
 193         unsigned int NoOfPixels, NoOfRecSets;
 194
 195         //a reconstruction image for the intermediate reconstruction during
 196         //the iterative process
 197         REAL * bhcRecBuffer;
 198
 199         /* make sure that the nergy prjfil's are present,
 200          * this projfil's contain the origila projection data corrected by multi-
 201          * plying by the correction coefficient provided by the user in the input
 202          * file
 203          */
 204
 205         /* create nergy input files that use the same structures as the original
 206          * input file, use PROJECTION OLD and run convolution algorithm for
 207          * reconstruction
 208          */
 209
 210         //create and open input files
 211         //copy the original input file to the temp directory and open it for reading
 212         sprintf(command, " cp file11 %s/inputFile.in                    ", tempDir);
 213         if (system(command) != 0)
 214         {
 215                 fprintf(output,
 216                                 "\n **** ERROR: could not copy input file to inputFile.in\n "
 217                                                 " **** program aborted\n\n");
 218                 exit(-1);
 219         }
 220
 221         //open a new input file for writing
 222
 223         sprintf(filename, "%s/bhc_in.in", tempDir);
 224         bhcInputFromPrjfil.open(filename);
 225         if (!bhcInputFromPrjfil.is_open())
 226         {
 227                 fprintf(output, "\n **** ERROR: could not open file %s from "
 228                                 "temp directory\n  **** program aborted\n\n", filename);
 229                 exit(-1);
 230         }
 231
 232         //create input files that use the existing prjfil's to perform the reconstruction
 233         bhcInputFromPrjfil << "PICTURE RECONSTUCTION " << GeoPar.nelem << " "
 234                         << GeoPar.pixsiz << "\n";
 235         bhcInputFromPrjfil << "PROJECTION OLD \n";
 236         bhcInputFromPrjfil << "EXECUTE CONV \n";
 237         bhcInputFromPrjfil << "monochromatic intermediate reconstruction\n";
 238         bhcInputFromPrjfil << "hamming 0.54 2 \n";
 239         bhcInputFromPrjfil << "END\n";
 240
 241         //close the input file
 242         bhcInputFromPrjfil.close();
 243
 244         //allocalte needed memory
 245         //bhcRecfils = new DIGFileSnarkRec[Spctrm.nergy];
 246         //bhcRec_buffer = new REAL * [Spctrm.nergy];
 247         NoOfPixels = GeoPar.nelem * GeoPar.nelem;
 248
 249         /* repeat noIter times
 250          */
 251         for (r = 0; r < noIter; r++)
 252         {
 253                 sprintf(filename, "%s/inputFile.in", tempDir);
 254                 inputFile.open(filename);
 255                 if (!inputFile.is_open())
 256                 {
 257                         fprintf(output,
 258                                         "\n **** ERROR: could not read the inputFile.in from "
 259                                                         "temp directory\n  **** program aborted\n\n");
 260                         exit(-1);
 261                 }
 262
 263                 /* run snark14 and create new input files based on reconstruction images,
 264                  * then create monochromatic projections (nergy of them) and polychromatic
 265                  * projections in order to perform the next step of iterative data refinement
 266                  *
 267                  * correct the data in projfil's and prepare the input files for the next
 268                  * round of iterative data refinement step
 269                  */
 270
 271                 //run snark14 to obtain reconstruction in recfil
 272                 //for (i = 0; i < Spctrm.nergy; i++) {
 273                 //rename the existing q_p%i to prjfil
 274                 sprintf(command, "cd %s\n cp q_p%i  prjfil", tempDir, r);
 275                 if (system(command) != 0)
 276                 {
 277                         fprintf(output, "\n **** ERROR: could not execute  command "
 278                                         "%s\n  **** program aborted\n\n", command);
 279                         exit(-1);
 280                 }
 281                 //run snark14 with input file in"i".in
 282                 sprintf(command, "cd %s\n snark14 bhc_in.in >out.out", tempDir);
 283                 if (system(command) != 0)
 284                 {
 285                         fprintf(output, "\n **** ERROR: could not execute  command "
 286                                         "%s\n  **** program aborted\n\n", command);
 287                         exit(-1);
 288                 }
 289                 //rename the recfil, so it is not lost at the next run of snark14
 290                 sprintf(command, "cd %s\n mv recfil recfil%i", tempDir, r);
 291                 if (system(command) != 0)
 292                 {
 293                         fprintf(output, "\n **** ERROR: could not execute  command "
 294                                         "%s\n  **** program aborted\n\n", command);
 295                         exit(-1);
 296                 }
 297                 //}
 298
 299                 //open just created recfil's and prepare data for reading
 300
 301                 sprintf(filename, "%s/recfil%i", tempDir, r);
 302                 if (bhcRecfil.Open(filename) != 0)
 303                 {
 304                         fprintf(output, "\n **** ERROR: could not open file "
 305                                         "%s (r = %i)\n  **** program aborted\n\n", filename, r);
 306                         exit(-1);
 307                 }
 308
 309                 // verify that there is only one reconstruction set (convolution
 310                 // produces only one)
 311
 312                 if (bhcRecfil.GetNoOfRecSets(&NoOfRecSets) != 0)
 313                 {
 314                         fprintf(output,
 315                                         "\n **** ERROR: could not read the number of reconstruction"
 316                                                         " sets from file %s\n  **** program aborted\n\n",
 317                                         filename);
 318                         exit(-1);
 319                 }
 320                 if (NoOfRecSets != 1)
 321                 {
 322                         fprintf(output, "\n **** ERROR: the number (= %i) of reconstruction"
 323                                         " sets from file %s is incorrect, should be =1\n  "
 324                                         "**** program aborted\n\n", NoOfRecSets, filename);
 325                         exit(-1);
 326                 }
 327
 328                 //allocate space for reconstuction data
 329                 bhcRecBuffer = new REAL[NoOfPixels];
 330
 331                 //select the reconstruction
 332                 if (bhcRecfil.SelectRec(0) != 0)
 333                 {
 334                         fprintf(output, "\n **** ERROR: could not select reconstruction #0 "
 335                                         " from file %s "
 336                                         "**** program aborted\n\n", filename);
 337                         exit(-1);
 338                 }
 339
 340                 //read in the data from the reconstruction pixel by pixel
 341                 bhcRecfil.GetRecData(bhcRecBuffer);
 342
 343                 //close the recfil
 344                 bhcRecfil.Close();
 345
 346                 //}
 347
 348                 //  create new set of input files for monochromatic projections (overwrite the old ones)
 349                 //with structures being pixel by pixel
 350                 //description based on the reconstruction images from recfil's above;
 351                 //  create an input file for polychromatic projection with density
 352                 //values based on densities obtained in the reconstruction recfil's
 353                 //and distibution based on the original input file
 354
 355                 //open the monochromatic input file for writing
 356                 //and the polychromatic input file for writing
 357                 for (i = 0; i < 2; i++)
 358                 {
 359                         sprintf(filename, "%s/%iin_rec%i.in", tempDir, r, i);
 360                         bhcInputFiles[i].open(filename);
 361                         if (!bhcInputFiles[i].is_open())
 362                         {
 363                                 fprintf(output, "\n **** ERROR: could not open file %s from "
 364                                                 "temp directory\n  **** program aborted\n\n", filename);
 365                                 exit(-1);
 366                         }
 367                 }
 368
 369                 //start copying initial lines from the original input file to
 370                 //the new input files
 371                 //  until lines
 372                 //       SPECTRUM POLYCHROMATIC imax
 373                 //       d1 p1 d2 p2  ... d1imax  pimax
 374                 // are encountered
 375                 for (i = 0; i < 2; i++)
 376                 {
 377                         bhcInputFiles[i] << "CREATE\n";
 378                 }
 379
 380                 while (!inputFile.eof())
 381                 {
 382                         inputFile.getline(line, 255);
 383                         if (!strncasecmp(line, "SPEC", 4))
 384                         {        //if line starts with the word "SPECTRUM"
 385                                 break; //jump out of the loop
 386                         }
 387                         //otherwise write the line to all the newly created input files
 388                         for (i = 0; i < 2; i++)
 389                         {
 390                                 bhcInputFiles[i] << line << "\n";
 391                         }
 392                 }
 393                 //write the line
 394                 //       SPECTRUM MONOCHROMATIC "Spctrm.energy[i]
 395                 // to the monochromatic input files
 396
 397                 bhcInputFiles[0] << "SPECTRUM MONOCHROMATIC   " << Spctrm.energy[0]
 398                                 << "\n";
 399
 400                 //write the lines
 401                 //       SPECTRUM POLYCHROMATIC "Spctrm.nergy"
 402                 //       Spctrm.energy[0]  Spctrm.engwt[0] ... Spctrm.energy[nergy-1]  Spctrm.engwt[nergy-1]
 403                 // to the polychromatic input file
 404                 bhcInputFiles[1] << "SPECTRUM POLYCHROMATIC  " << Spctrm.nergy << "\n";
 405                 for (i = 0; i < Spctrm.nergy; i++)
 406                 {
 407                         bhcInputFiles[1] << Spctrm.energy[i] << "  "
 408                                         << 100 * Spctrm.engwt[i] << "   ";
 409                 }
 410                 bhcInputFiles[1] << "\n";
 411                 //write the line
 412                 //       OBJECTS
 413                 // to the input files
 414                 for (i = 0; i < 2; i++)
 415                 {
 416                         bhcInputFiles[i] << "OBJECTS" << "\n";
 417                 }
 418
 419                 //write the pixel by pixel description of the reconstruction image
 420                 // to the input files
 421                 INTEGER sizeHalf = 0.5 * GeoPar.nelem;
 422
 423                 for (k = 0; k < NoOfPixels; k++)
 424                 {
 425                         x = (k % GeoPar.nelem - sizeHalf) * GeoPar.pixsiz;
 426                         y = (k / GeoPar.nelem - sizeHalf) * GeoPar.pixsiz;
 427                         u = GeoPar.pixsiz / 2;
 428                         den = bhcRecBuffer[k];
 429
 430                         bhcInputFiles[0] << "RECT\t" << x << "\t" << y << "\t" << u << "\t"
 431                                         << u << "\t" << 0 << "\t" << den << "\n";
 432
 433                         bhcInputFiles[1] << "RECT\t" << x << "\t" << y << "\t" << u << "\t"
 434                                         << u << "\t" << 0 << "\t";
 435
 436                         bhcInputFiles[1] << den << "\nDENS";
 437
 438                         //compute the densities using piecewise linear conversion function
 439                         for (i = 1; i < Spctrm.nergy; i++)
 440                         {
 441                                 //REAL denden = 1.0;
 442                                 REAL nextden = 0.0;
 443                                 REAL m, mm;     //slope and y intercept of the line segment
 444
 445                                 for (j = 1; j < noPointsConv; j++)
 446                                 {
 447                                         if (den < c[i][j].x
 448                                                         || (j == noPointsConv - 1 && den > c[i][j].x))
 449                                         {
 450                                                 m = (c[i][j].y - c[i][j - 1].y)
 451                                                                 / (c[i][j].x - c[i][j - 1].x);
 452                                                 mm = c[i][j].y - m * c[i][j].x;
 453                                                 nextden = m * den + mm;
 454                                                 break;
 455                                         }
 456                                         if (den == c[i][j].x)
 457                                         {
 458                                                 nextden = c[i][j].y;
 459                                                 break;
 460                                         }
 461
 462                                 }
 463                                 bhcInputFiles[1] << "\t" << nextden;
 464                         }
 465
 466                         bhcInputFiles[1] << "\n";
 467                 }
 468
 469                 //write the lines in input files up to line
 470                 //      PROJECTION PSEUDO
 471                 for (i = 0; i < 2; i++)
 472                 {
 473                         bhcInputFiles[i] << "LAST 1 \n" << "PHANTOM  AVERAGE "
 474                                         << Creacm.nave1 << "\n" << GeoPar.nelem << "    "
 475                                         << GeoPar.pixsiz << "\n" << "RAYSUM \n" << "PICTURE TEST \n"
 476                                         << "PROJECTION BHPSEUDO \n"
 477                                         << "psudo run to create prjfil for correction\n";
 478                 }
 479
 480                 //skip through lines of the original input file until line
 481                 //   GEOMETRY ....
 482                 //is encountered
 483                 while (!inputFile.eof())
 484                 {
 485                         inputFile.getline(line, 256);
 486                         if (!strncasecmp(line, "GEOM", 4))
 487                         {        //if line starts with the word "GEOM"
 488                                 break; //jump out of the loop
 489                         }
 490                 }
 491
 492                 //copy the lines from the original input file to all the new
 493                 // input files until the line
 494                 //    BACKGROUND ...
 495                 //is encountered
 496                 for (i = 0; i < 2; i++)
 497                 {
 498                         bhcInputFiles[i] << line << "\n";
 499                 }
 500                 while (!inputFile.eof())
 501                 {
 502                         inputFile.getline(line, 255);
 503                         //write the line to all the newly created input files
 504                         for (i = 0; i < 2; i++)
 505                         {
 506                                 bhcInputFiles[i] << line << "\n";
 507                         }
 508                         if (!strncasecmp(line, "BACK", 4))
 509                         {        //if line starts with the word "BACKGROUND"
 510                                 break; //jump out of the loop
 511                         }
 512                 }
 513
 514                 //finish the files with END line and close the files
 515                 for (i = 0; i < 2; i++)
 516                 {
 517                         bhcInputFiles[i] << "END \n";
 518                         bhcInputFiles[i].close();
 519                 }
 520                 //close the original input file
 521                 inputFile.close();
 522
 523                 //obtain projection data from the input files created above
 524
 525                 for (i = 0; i < 2; i++)
 526                 {
 527
 528                         //run snark14 with input file %rin"i".in
 529                         sprintf(command, "cd %s\n snark14 %iin_rec%i.in >out%i.out",
 530                                         tempDir, r, i, i);
 531                         if (system(command) != 0)
 532                         {
 533                                 fprintf(output, "\n **** ERROR: could not execute  command"
 534                                                 "%s\n  **** program aborted\n\n", command);
 535                                 exit(-1);
 536                         }
 537                         //rename the prjfil, so it is not lost at the next run of snark14
 538                         if (i == 1)
 539                                 sprintf(command, "cd %s\n mv prjfil p%i", tempDir, r);
 540                         else
 541                                 sprintf(command, "cd %s\n mv prjfil m%i", tempDir, r);
 542                         if (system(command) != 0)
 543                         {
 544                                 fprintf(output, "\n **** ERROR: could not execute  command"
 545                                                 "%s\n  **** program aborted\n\n", command);
 546                                 exit(-1);
 547                         }
 548                 }
 549
 550                 //correct the projection data following the formula
 551                 //m[i] = miu q[p][i] + (m'[i] - miu q(p')[i])
 552                 idr_correction(r);
 553
 554         } //end of repeat noIter times
 555
 556         fprintf(output, "\n **** Finished the iterations\n\n");
 557
 558         /* write the final value of m[0] to prjfil that will be passed to the remainder
 559          * of the run of snark14
 560          */
 561
 562         sprintf(command, "cp %s/q_p%i prjfil", tempDir, noIter);
 563         if (system(command) != 0)
 564         {
 565                 fprintf(output, "\n **** ERROR: could not execute  command"
 566                                 "%s\n  **** program aborted\n\n", command);
 567                 exit(-1);
 568         }
 569
 570         //if (bhcRecfil != NULL) delete [] bhcRecfil;
 571         //if (bhc_InputFiles != NULL) delete [] bhc_InputFiles;
 572         if (bhcRecBuffer != NULL)
 573                 delete[] bhcRecBuffer;
 574         return -1;
 575 }
 576
 577 ///////////////////////////////////////////////////////////////////////////////
 578 ///////////////////////////////////////////////////////////////////////////////
 579 ///////////////////////////////////////////////////////////////////////////////
 580
 581 int bh_correction::idr_correction( INTEGER r)
 582 {
 583         //correct the projection data following the formula
 584         //m[i] = miu q[p][i] + (m'[i] - miu q(p')[i])
 585         REAL miu = 1.0; //the relaxation parameter, for now permanently set to 1.0
 586         unsigned INTEGER i, j, k;  //loop counters
 587         DIGFileSnarkProj Org, NewM, NewP, Corrected;
 588         char filename[256];
 589         //char command[256];
 590         REAL** Buffer = new REAL*[4]; //buffer to store a projection data
 591
 592         unsigned int NoOfRays;
 593         unsigned int NoOfProjections;
 594         ProjFileMH MainHeader;
 595         ProjFileAH ApplicationHeader;
 596         ProjFilePrjH ProjectionHeader[3];
 597
 598         //open original prjfil corrected by coefficient at energy level 0
 599         //   bhc/q_p0
 600         //read Main Header and Application Header of the file
 601         sprintf(filename, "%s/q_p0", tempDir);
 602         if (Org.Open(filename) != 0)
 603         {
 604                 fprintf(output, "\n **** ERROR: could not open file "
 605                                 "%s (Org)\n  **** program aborted\n\n", filename);
 606                 exit(-1);
 607         }
 608
 609         if (Org.GetNoOfRays(&NoOfRays) != 0)
 610         {
 611                 fprintf(output,
 612                                 "\n **** ERROR: could not read the number of rays from  file "
 613                                                 "%s\n  **** program aborted\n\n", filename);
 614                 exit(-1);
 615         }
 616
 617         if (Org.GetNoOfProjs(&NoOfProjections) != 0)
 618         {
 619                 fprintf(output,
 620                                 "\n **** ERROR: could not read the number of projections from  file "
 621                                                 "%s\n  **** program aborted\n\n", filename);
 622                 exit(-1);
 623         }
 624
 625         if (Org.GetMainHeader(&MainHeader) != 0)
 626         {
 627                 fprintf(output, "\n **** ERROR: could not read Main Header of  file "
 628                                 "%s\n  **** program aborted\n\n", filename);
 629                 exit(-1);
 630         }
 631
 632         if (Org.GetAppHeader(&ApplicationHeader) != 0)
 633         {
 634                 fprintf(output,
 635                                 "\n **** ERROR: could not read Application Header of  file "
 636                                                 "%s\n  **** program aborted\n\n", filename);
 637                 exit(-1);
 638         }
 639
 640         //open just created monochromatic prjfil at energy level i
 641         //   bhc/m%r
 642         sprintf(filename, "%s/m%i", tempDir, r);
 643         if (NewM.Open(filename) != 0)
 644         {
 645                 fprintf(output, "\n **** ERROR: could not open file "
 646                                 "%s\n  **** program aborted\n\n", filename);
 647                 exit(-1);
 648         }
 649
 650         //open just created polychromatic prjfil
 651         //   bhc/p%r
 652         sprintf(filename, "%s/p%i", tempDir, r);
 653         if (NewP.Open(filename) != 0)
 654         {
 655                 fprintf(output, "\n **** ERROR: could not open file "
 656                                 "%s\n  **** program aborted\n\n", filename);
 657                 exit(-1);
 658         }
 659
 660         //open prjfil that will contain corrected data,
 661         //   bhc/q_p%r
 662         //set its Main Header and Application Header
 663         sprintf(filename, "%s/q_p%i", tempDir, r + 1);
 664         if (Corrected.Open(filename, &MainHeader) != 0)
 665         {
 666                 fprintf(output, "\n **** ERROR: could not open file "
 667                                 "%s\n  **** program aborted\n\n", filename);
 668                 exit(-1);
 669         }
 670
 671         if (Corrected.SetAppHeader(&ApplicationHeader) != 0)
 672         {
 673                 fprintf(output,
 674                                 "\n **** ERROR: could not set Application Header of file "
 675                                                 "%s\n  **** program aborted\n\n", filename);
 676                 exit(-1);
 677         }
 678
 679         //allocate memory to store a single projection data
 680         for (j = 0; j < 4; j++)
 681         {
 682                 (Buffer[j]) = new REAL[NoOfRays];
 683         }
 684
 685         for (j = 0; j < NoOfProjections; j++)
 686         {
 687
 688                 // array set selection
 689                 if (Org.SelectProj(j) != 0)
 690                 {
 691                         fprintf(output, "\n **** ERROR: could not select projection of file"
 692                                         "Org\n  **** program aborted\n\n");
 693                         exit(-1);
 694                 }
 695                 //get projection header
 696                 if (Org.GetProjHeader(&ProjectionHeader[0]) != 0)
 697                 {
 698                         fprintf(output,
 699                                         "\n **** ERROR: could not get projection header from file"
 700                                                         "Org\n  **** program aborted\n\n");
 701                         exit(-1);
 702                 }
 703                 //get projection data
 704                 if (Org.GetProjData(Buffer[0]) != 0)
 705                 {
 706                         fprintf(output,
 707                                         "\n **** ERROR: could not get projection data from file"
 708                                                         "Org\n  **** program aborted\n\n");
 709                 }
 710
 711                 // array set selection
 712                 if (NewM.SelectProj(j) != 0)
 713                 {
 714                         fprintf(output, "\n **** ERROR: could not select projection of file"
 715                                         "NewM\n  **** program aborted\n\n");
 716                         exit(-1);
 717                 }
 718
 719                 //get projection data
 720                 if (NewM.GetProjData(Buffer[1]) != 0)
 721                 {
 722                         fprintf(output,
 723                                         "\n **** ERROR: could not get projection data from file"
 724                                                         "NewM\n  **** program aborted\n\n");
 725                 }
 726
 727                 // array set selection
 728                 if (NewP.SelectProj(j) != 0)
 729                 {
 730                         fprintf(output, "\n **** ERROR: could not select projection of file"
 731                                         "NewP\n  **** program aborted\n\n");
 732                         exit(-1);
 733                 }
 734
 735                 //get projection data
 736                 if (NewP.GetProjData(Buffer[2]) != 0)
 737                 {
 738                         fprintf(output,
 739                                         "\n **** ERROR: could not get projection data from file"
 740                                                         "NewP\n  **** program aborted\n\n");
 741                 }
 742
 743                 //calculate corrected value for each ray in the projection
 744
 745                 for (k = 0; k < NoOfRays; k++)
 746                 {
 747
 748                         REAL den = Buffer[2][k];
 749                         REAL denden = 1.0;
 750                         REAL newden = 0.0;
 751                         for (i = 0; i <= noDegree; i++)
 752                         {
 753                                 newden += denden * a[i];
 754                                 denden = denden * den;
 755                         }
 756
 757                         Buffer[3][k] = miu * Buffer[0][k] + Buffer[1][k] - miu * newden;
 758                 }
 759
 760                 //add corrected projection to the Corrected file
 761                 if (Corrected.AppendProj(&ProjectionHeader[0], Buffer[3]) != 0)
 762                 {
 763                         fprintf(output, "\n **** ERROR: could not append projection to file"
 764                                         "Corrected\n  **** program aborted\n\n");
 765                 }
 766         }
 767
 768         //close the files
 769         Org.Close();
 770         NewM.Close();
 771         NewP.Close();
 772         Corrected.Close();
 773
 774         //delete allocated memory
 775         if (Buffer[0] != NULL)
 776                 delete[] Buffer[0];
 777         if (Buffer[1] != NULL)
 778                 delete[] Buffer[1];
 779         if (Buffer[2] != NULL)
 780                 delete[] Buffer[2];
 781         if (Buffer[3] != NULL)
 782                 delete[] Buffer[3];
 783         delete Buffer;
 784
 785         return 0;
 786 }
 787