src/snark/rdproj.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/rdproj.cpp $
   5  $LastChangedRevision: 91 $
   6  $Date: 2014-07-02 17:34:27 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  RDPROJ GENERATES THE RANDOM ACCESS VERSION OF THE PROJECTION DATA
  11  ON  PRJFIL.
  12  THE MODIFIER ON  PROJECTION CARD IS ONE OF THE FOLLOWING -
  13  REAL     THIS IS THE DEFAULT VALUE.  PROJECTION DATA ARE READ FROM
  14  FILE11  AND PUT ON PRJFIL
  15  PSEUDO   CREATES PSEUDO PROJECTION DATA FROM INFORMATION ON FILE11
  16  AND INPUT  AND PUTS ON PRJFIL.  THIS REQUIRES THAT
  17  PICTURE  TEST  BE ALREADY SPECIFIED.
  18  OLD      THIS ASSUMES THAT PRJFIL ALREADY EXISTS FROM A PREVIOUS
  19  SNARK RUN.  THE GEOMETRY AND NOISE INFORMATION ARE IN THE
  20  FIRST RECORD OF  PRJFIL.
  21
  22  SEE SNARK MANUAL FOR MORE DETAILS.
  23  *
  24  *
  25  *Modification date:  July 13, 2008
  26  *Author: Joanna Klukowska (jk)
  27  *Changes:  Added functionality for beam hardening correction.
  28  *          If the user sets the options for correction in the input
  29  *          file, rdproj() createas nergy temporary prjfil's with to
  30  *          pass to bh_correction::correction() for beam hardening correction
  31  *          which is responsible for creation of the final prjfil
  32  *          that is used in the remainder of the program.
  33  *          The option of beam hardening correction is only appropriate
  34  *          when polychromatic data is provided.
  35  *
  36  */
  37
  38 #include <cstdlib>
  39 #include <cstdio>
  40 #include <cmath>
  41 #include <unistd.h>
  42 #include <sys/stat.h>
  43 #include <sys/types.h>
  44 #include <errno.h>
  45
  46 #include "blkdta.h"
  47 #include "creacm.h"
  48 #include "geom.h"
  49 #include "anglst.h"
  50 #include "spctrm.h"
  51 #include "consts.h"
  52 #include "uiod.h"
  53
  54 #include "raylen.h"
  55 #include "posit.h"
  56 #include "DIGGauss.h"
  57 #include "creaer.h"
  58 #include "bldlst.h"
  59 #include "pseudo.h"
  60 #include "bh_correction.h"
  61
  62 #include "recfile.h"
  63 #include "file11.h"
  64 #include "projfile.h"
  65 #include "infile.h"
  66
  67 #include "rdproj.h"
  68 #include "bh_correction.h"
  69
  70 #include <time.h>
  71
  72 void rdproj()
  73 {
  74
  75         static const INTEGER proj_codes[4] =
  76         { CHAR2INT('o', 'l', 'd', ' '), CHAR2INT('r', 'e', 'a', 'l'), CHAR2INT('p',
  77                         's', 'e', 'u'), CHAR2INT('b', 'h', 'p', 's')    //jk 9/19/2008
  78                         //used only for recursive runs of snark during the
  79                         //beam hardening correction process
  80                         //should not be listed as an option in the manual
  81                         };
  82
  83         // jk 7/13/2008 added for beam hardening correction
  84         static const INTEGER bhc_codes[1] =
  85         { CHAR2INT('b', 'e', 'a', 'm') };
  86
  87         static BOOLEAN called = FALSE;
  88
  89         REAL dtor;
  90         INTEGER nword;
  91         INTEGER nr;
  92         REAL umean;
  93         REAL amean;
  94         REAL totden;
  95         REAL totlen;
  96         REAL hfside;
  97         INTEGER np;
  98         REAL sinth;
  99         REAL theta;
 100         REAL costh;
 101         REAL ang;
 102         REAL the;
 103         REAL angle;
 104         INTEGER numb;
 105         REAL snorm;
 106
 107         INTEGER* lbase = NULL;
 108         REAL* wbase = NULL; //wei 3/2005
 109
 110         int i, j, k;
 111
 112         BOOLEAN eol;
 113         REAL ax, ay, mx, my;
 114
 115         // jk 7/13/2008 added for beam hardening correction
 116
 117         INTEGER bhcword; //BHC keyword read from the input file
 118         SnarkPrjFile * BHC_PrjFile = NULL; //temporary prjfil's
 119         char filename[256];
 120         char command[256];
 121         bh_correction bh_cor;
 122
 123         ///////////////////////////////////////////////////////////
 124
 125         // CHECK FOR MULTIPLE CALLS
 126
 127         if (called)
 128                 creaer(1);
 129
 130         called = TRUE;
 131         dtor = (REAL) (Consts.pi / 180.0);
 132
 133         // CHECK THAT RDPICT HAS BEEN PERFORMED BEFORE THIS
 134
 135         if (GeoPar.nelem <= 0)
 136                 creaer(23);
 137
 138         // DETERMINE THE PROJECTION MODIFIER AND THE INPUT FILE
 139
 140         // changed last parameter from "4" to "3". lajos, Nov 18, 2004
 141         nword = InFile.getwrd(FALSE, &eol, proj_codes, 4);
 142
 143         // IF OLD  IS OPTED JUST READ HEADER & QUIT
 144
 145         if (nword == proj_codes[0]) // OLD
 146         {
 147                 // OPEN PROJECTION FILE FOR PROJECTION OLD CASE
 148
 149                 if (ProjFile.Open((char *) "prjfil") != 0)
 150                 {
 151                         fprintf(output, "\n **** unable to open prjfil");
 152                         fprintf(output, "\n **** program aborted\n");
 153                         exit(-1);
 154                 };
 155
 156                 ProjFile.ReadHeader(Creacm.name, &NoisePar, &Spctrm, &Anglst);
 157
 158                 // compute average density - fixes bug 232 of snark05
 159                 double *projdata = new double[GeoPar.nrays];
 160                 totden = 0.0;
 161                 totlen = 0.0;
 162                 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
 163                 for (np = 0; np < GeoPar.prjnum; ++np)
 164                 {
 165                         ProjFile.ReadProj(np, projdata, GeoPar.nrays);
 166                         for (nr = 0; nr < GeoPar.nrays; nr++)
 167                         {
 168                                 if (GeoPar.strip)
 169                                 {
 170                                         if (GeoPar.vri)
 171                                         {
 172                                                 totden +=
 173                                                                 (REAL) (projdata[nr]
 174                                                                                 / (GeoPar.pinc
 175                                                                                                 * MAX0(fabs(sinth), fabs(costh))));
 176                                         }
 177                                         else
 178                                         {
 179                                                 totden += projdata[nr] / GeoPar.pinc;
 180                                         }
 181                                 }
 182                                 else
 183                                 {
 184                                         totden += projdata[nr];
 185                                 }
 186                                 posit(np, nr, &ax, &ay, &mx, &my);
 187                                 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
 188                                                 hfside, 1.0, 0.0);
 189                         }
 190                 }
 191                 GeoPar.aveden = totden / totlen;
 192                 fprintf(output, "\n         estimate of totlen = %16.6f", totlen);
 193                 fprintf(output, "\n         estimate of totden = %16.6f", totden);
 194                 fprintf(output, "\n         estimate of average density = %10.4f",
 195                                 GeoPar.aveden);
 196                 delete[] projdata;
 197         }
 198
 199         else
 200         {
 201                 // READ PREFACE TO PROJECTION DATA
 202
 203                 if (nword == proj_codes[2] || nword == proj_codes[3]) // PSEUDO
 204                 {
 205                         NoisePar.resetNoise(); // jklukowska, bug 263, in case the noise was set
 206                                                                    // previously, we need to reset it to initial values
 207                         InFile.CreaInProj();
 208                 }
 209                 else
 210                 {
 211                         if (nword == proj_codes[1]) // REAL
 212                         {
 213                                 if (!File11.isOpen)
 214                                 {
 215                                         if (File11.Open(FALSE) != 0)
 216                                         {
 217                                                 fprintf(output,
 218                                                                 "\n **** unable to open file11 for reading");
 219                                                 fprintf(output, "\n **** program aborted\n");
 220                                                 exit(-1);
 221                                         };
 222                                         File11.isOpen = TRUE;
 223                                 };
 224
 225                                 File11.CreaInProj();
 226                         }
 227                         else
 228                         {
 229                                 fprintf(output,
 230                                                 "\n **** OLD, REAL or PSEUDO must be specified");
 231                                 fprintf(output, "\n **** program aborted\n");
 232                                 exit(-1);
 233                         }
 234                 };
 235
 236                 if (Creacm.erflag)
 237                         exit(555);
 238
 239                 // jk 7/13/2008 added for beam hardening correction
 240                 //check if bh_correction flag is set
 241                 bhcword = InFile.getwrd(FALSE, &eol, bhc_codes, 1);
 242                 if (bhcword == bhc_codes[0])
 243                 {        //beam hardening correction flag set
 244                         bh_cor.bhc = TRUE;
 245                         //check if data is polychromatic
 246                         if (Spctrm.nergy == 1) //bhc cannot be performed
 247                         {
 248                                 fprintf(output,
 249                                                 "\n **** Beam hardening correction is not appropriate\n"
 250                                                                 "for monochromatic data (nergy = 1).\n"
 251                                                                 "Program will continue ignoring BHC command.\n\n");
 252                                 bh_cor.bhc = FALSE;
 253                         }
 254                         else
 255                         {
 256
 257                                 //read in the parameters for beam hardening correction
 258                                 bh_cor.readInputFilePoly();
 259
 260                                 //allocate temporary prjfil's
 261                                 BHC_PrjFile = new SnarkPrjFile[1];
 262
 263                                 time_t t;
 264                                 time(&t);   //get the current time
 265                                 //create temporary directory
 266                                 sprintf(bh_cor.tempDir, "bhc%i", (int) t);
 267                                 sprintf(command, "mkdir %s", bh_cor.tempDir);
 268                                 if (system(command) != 0)
 269                                 {
 270                                         fprintf(stdout, "\nERROR creating temporary directory %s\n",
 271                                                         filename);
 272                                         exit(1);
 273                                 }
 274                         }
 275                 }
 276
 277                 if (ProjFile.Open((char*) "prjfil", GeoPar.usrays, GeoPar.prjnum) != 0)
 278                 {
 279                         fprintf(output, "\n **** unable to open prjfil");
 280                         fprintf(output, "\n **** program aborted\n");
 281                         exit(-1);
 282                 };
 283
 284                 // jk 7/13/2008 added for beam hardening correction
 285                 //open temporary prjfil's
 286                 if (bh_cor.bhc)
 287                 {
 288                         for (i = 0; i < 1; i++)
 289                         {
 290                                 sprintf(filename, "%s/q_p0", bh_cor.tempDir);
 291                                 if (BHC_PrjFile[i].Open(filename, GeoPar.usrays, GeoPar.prjnum)
 292                                                 != 0)
 293                                 {
 294                                         fprintf(output, "\n **** unable to open temporary prjfil");
 295                                         fprintf(output, "\n **** program aborted\n");
 296                                         exit(-1);
 297                                 }
 298                         }
 299                 }
 300
 301                 ProjFile.WriteHeader(Creacm.name, &GeoPar, &NoisePar, &Spctrm);
 302
 303                 ProjFile.WriteAngles(Anglst.bth);
 304
 305                 // jk 7/13/2008 added for beam hardening correction
 306                 // write header info to temporary prjfil's
 307                 if (bh_cor.bhc)
 308                 {
 309                         for (i = 0; i < 1; i++)
 310                         {
 311                                 BHC_PrjFile[i].WriteHeader(Creacm.name, &GeoPar, &NoisePar,
 312                                                 &Spctrm);
 313                                 BHC_PrjFile[i].WriteAngles(Anglst.bth);
 314                         }
 315                 }
 316
 317                 // ALLOCATE SPACE FOR ARRAYS LIST,WEIGHT,SPHI,CPHI,PBASE
 318
 319                 bldlst(&lbase, &wbase);
 320
 321                 if (GeoPar.div)
 322                 {
 323                         Anglst.genphi();
 324                 }
 325
 326                 Anglst.pbase = new REAL[GeoPar.nrays];
 327
 328                 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 329                 {
 330                         Anglst.pbase[nr] = 0.0;
 331                 }
 332
 333                 // IF PROJECTION DATA ARE NOISY WITH BIAS DUE TO SPECIFICATION OF
 334                 // SKEWED DISTRIBUTIONS PREPARE HERE TO REMOVE BIAS
 335
 336                 umean = 1.0;
 337
 338                 if (NoisePar.ultnmn != 0.0)
 339                 {
 340                         umean = (REAL) 1.0 / NoisePar.ultnmn;
 341                 }
 342
 343                 amean = 0.0;
 344
 345                 if (NoisePar.addnfl)
 346                 {
 347                         amean = NoisePar.addnmn;
 348                 }
 349
 350                 // INITIALIZE FOR COMPUTATION OF AVERAGE DENSITY
 351
 352                 totden = 0.0;
 353                 totlen = 0.0;
 354                 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
 355
 356                 // CHECK TO SEE IF PHANTOM IS PRESENT WHEN PSEUDO IS OPTED
 357
 358                 if ((nword == proj_codes[2] || nword == proj_codes[3])
 359                                 && (Creacm.test == NULL))
 360                 {
 361                         creaer(24);
 362                 }
 363
 364                 // SELECT PROJECTION BY PROJECTION AND PREPARE  PRJFIL
 365
 366                 for (np = 0; np < GeoPar.prjnum; np++)
 367                 {
 368                         if (nword == proj_codes[1])
 369                         {
 370
 371                                 if (!File11.ReadProjection(&theta, Anglst.pbase + GeoPar.fusray,
 372                                                 GeoPar.lusray - GeoPar.fusray + 1))
 373                                 {
 374                                         fprintf(output,
 375                                                         "\n **** end of file encountered on file11");
 376                                         fprintf(output, "\n **** program aborted\n");
 377                                         if (lbase != NULL)
 378                                                 delete[] lbase;
 379                                         if (wbase != NULL)
 380                                                 delete[] wbase; //wei 3/2005
 381                                         exit(-1);
 382                                 }
 383
 384                                 sinth = (REAL) sin(theta);
 385                                 costh = (REAL) cos(theta);
 386                                 ang = theta / dtor;
 387
 388                                 the = *(Anglst.bth + np);
 389
 390                                 angle = the / dtor;
 391                                 if ((REAL) fabs(theta - the) > (REAL) 0.001)
 392                                 {
 393                                         fprintf(output, "\n **** angles specified for projection "
 394                                                         "%4i do not match", np);
 395                                         fprintf(output, "\n      angle given in header is %10.4f "
 396                                                         "degrees = %10.4f radians", angle, the);
 397                                         fprintf(output, "\n      angle given in data   is %10.4f "
 398                                                         "degrees = %10.4f radians", ang, theta);
 399                                         fprintf(output, "\n **** program aborted\n");
 400                                         if (lbase != NULL)
 401                                                 delete[] lbase;
 402                                         if (wbase != NULL)
 403                                                 delete[] wbase; //wei 3/2005
 404                                         exit(-1);
 405                                 }
 406                         }
 407                         else
 408                         {
 409                                 // PROJECTION  PSEUDO or BHPSEUDO  OPTION
 410
 411                                 theta = Creacm.pang[np];
 412
 413                                 sinth = (REAL) sin(theta);
 414                                 costh = (REAL) cos(theta);
 415
 416                                 // bug221 start - swr - 03/03/07
 417                                 for (nr = 0; nr < GeoPar.nrays; nr++)
 418                                         Anglst.pbase[nr] = 0.0;
 419
 420                                 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 421                                 {
 422                                         REAL temp;
 423                                         if (nword == proj_codes[3])
 424                                         {
 425                                                 //BHPSEUDO option
 426                                                 //jk 07/20/2008 modyfying how pseudo projections are
 427                                                 // computed - they need to be polychromatic for
 428                                                 //beam hardening correction
 429                                                 temp = 0.0;
 430                                                 for (j = 0; j < Spctrm.nergy; j++)
 431                                                 {
 432                                                         temp += Spctrm.engwt[j]
 433                                                                         * NoisePar.raysum(
 434                                                                                         pseudo(phantom_poly_array[j], np,
 435                                                                                                         nr, lbase, wbase, &numb,
 436                                                                                                         &snorm, GeoPar.line,
 437                                                                                                         FALSE));
 438                                                 }
 439                                         }
 440                                         else
 441                                         {
 442                                                 //PSEUDO option
 443                                                 temp = NoisePar.raysum(
 444                                                                 pseudo(Creacm.test, np, nr, lbase, wbase, &numb,
 445                                                                                 &snorm, GeoPar.line, FALSE));
 446                                         }
 447                                         if (NoisePar.quanin > 0 && NoisePar.quanin != 4)
 448                                                 temp *= NoisePar.quanmn;
 449                                         Anglst.pbase[nr] = NoisePar.applyQuantumNoise(np, nr, temp);
 450                                 }
 451
 452                                 // TAKE CARE OF NOISE
 453
 454                                 NoisePar.scatterNoise(Anglst.pbase);
 455
 456                                 for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 457                                 {
 458                                         REAL temp = NoisePar.computeAttenuation(np, nr,
 459                                                         Anglst.pbase[nr]);
 460                                         temp = NoisePar.multiplicativeNoise(temp);
 461                                         Anglst.pbase[nr] = NoisePar.additiveNoise(temp);
 462                                 }
 463                                 // bug221 end - swr - 03/03/07
 464                         }
 465                         // REMOVE BIAS IN PROJECTION DATA POSSIBLY INTRODUCED BY A SKEWED
 466                         // NOISE DISTRIBUTION
 467
 468                         if (NoisePar.ultnfl || NoisePar.addnfl)
 469                         {
 470                                 for (i = GeoPar.fusray; i <= GeoPar.lusray; i++)
 471                                 {
 472                                         Anglst.pbase[i] = (Anglst.pbase[i] - amean) * umean;
 473                                 }
 474                         }
 475
 476                         // ADD UP FOR ESTIMATION OF AVERAGE DENSITY FROM RAYSUMS
 477
 478                         for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 479                         {
 480                                 // Changed here for SNARK89 Version 1.01 7/13/90 to conform to manual
 481                                 if (GeoPar.strip)
 482                                 {
 483                                         if (GeoPar.vri)
 484                                         {
 485                                                 totden +=
 486                                                                 (REAL) (Anglst.pbase[nr]
 487                                                                                 / (GeoPar.pinc
 488                                                                                                 * MAX0(fabs(sinth), fabs(costh))));
 489                                         }
 490                                         else
 491                                         {
 492                                                 totden += Anglst.pbase[nr] / GeoPar.pinc;
 493                                         }
 494                                 }
 495                                 else
 496                                 {
 497                                         totden += Anglst.pbase[nr];
 498                                 }
 499
 500                                 posit(np, nr, &ax, &ay, &mx, &my);
 501
 502                                 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
 503                                                 hfside, 1.0, 0.0);
 504
 505                         }
 506
 507                         // jk 7/13/2008 added for beam hardening correction
 508                         // write data to temporary prjfil's if bhc is set
 509                         if (bh_cor.bhc)
 510                         {
 511                                 for (i = 0; i < 1; i++)
 512                                 {
 513
 514                                         //iterate over values of Anglst.pbase and multiply by
 515                                         //coefficients of the correction polynomial
 516                                         //for (nr = GeoPar.fusray; nr <= GeoPar.lusray; nr++)
 517                                         INTEGER NoOfRays = GeoPar.lusray - GeoPar.fusray + 1;
 518                                         REAL t, tt;
 519                                         REAL * temp = new REAL[NoOfRays];
 520                                         for (k = 0; k < NoOfRays; k++)
 521                                                 temp[k] = 0;
 522                                         for (nr = 0; nr < NoOfRays; nr++)
 523                                         {
 524                                                 t = Anglst.pbase[nr + GeoPar.fusray];
 525                                                 tt = 1;
 526                                                 for (j = 0; j <= bh_cor.noDegree; j++)
 527                                                 {
 528                                                         temp[nr] += bh_cor.a[j] * tt;
 529                                                         tt = tt * t;
 530                                                 }
 531
 532                                         }
 533                                         BHC_PrjFile[i].WriteProj(np, temp); //&(Anglst.pbase[GeoPar.fusray]));
 534                                         delete[] temp;
 535                                 }
 536                         }
 537                         else //write data to the final prjfil
 538                         {
 539
 540                                 // WRITE ON PRJFIL.
 541                                 ProjFile.WriteProj(np, &(Anglst.pbase[GeoPar.fusray]));
 542                         }
 543                 }
 544                 if (!bh_cor.bhc)
 545                 {
 546                         GeoPar.aveden = totden / totlen;
 547                         fprintf(output, "\n\n         estimate of totlen = %16.6f", totlen);
 548
 549                         fprintf(output, "\n         estimate of totden = %16.6f", totden);
 550
 551                         fprintf(output, "\n         estimate of average density = %10.4f",
 552                                         GeoPar.aveden);
 553                 }
 554                 // WRITE UPDATED HEADER ON PRJFIL
 555
 556                 // jk 7/13/2008 added for beam hardening correction
 557                 // write header to temporary prjfil's if bhc is set
 558                 if (bh_cor.bhc)
 559                 {
 560                         for (i = 0; i < 1; i++)
 561                         {
 562                                 BHC_PrjFile[i].WriteHeader(Creacm.name, &GeoPar, &NoisePar,
 563                                                 &Spctrm);
 564                         }
 565                 }
 566                 else //write header to the final prjfil
 567                 {
 568
 569                         ProjFile.WriteHeader(Creacm.name, &GeoPar, &NoisePar, &Spctrm);
 570                 }
 571
 572                 // PREPARE TO EXIT FROM ROUTINE
 573
 574                 //free(Anglst.bth); deleted now in the destructor
 575         }
 576
 577         // CLOSE THE PRJFIL
 578
 579         // jk 7/13/2008 added for beam hardening correction
 580         if (bh_cor.bhc)
 581         {
 582                 for (i = 0; i < 1; i++)
 583                 {
 584                         BHC_PrjFile[i].Close();
 585                 }
 586         }
 587         else
 588         {
 589
 590                 ProjFile.Close();
 591         }
 592
 593         // jk 7/13/2008 added for beam hardening correction
 594         // CALL BEAM HARDENING CORRECTION ROUTINES
 595         if (bh_cor.bhc)
 596         {
 597                 ProjFile.Close();
 598                 bh_cor.correction(BHC_PrjFile);
 599
 600                 // OPEN PROJECTION FILE JUST CORRECTED FOR BEAM HARDENING
 601
 602                 if (ProjFile.Open((char*) "prjfil") != 0)
 603                 {
 604                         fprintf(output,
 605                                         "\n **** unable to open prjfil after beam hardening");
 606                         fprintf(output, "\n **** program aborted\n");
 607                         exit(-1);
 608                 };
 609
 610                 ProjFile.ReadHeader(Creacm.name, &NoisePar, &Spctrm, &Anglst);
 611
 612                 // compute average density
 613                 double *projdata = new double[GeoPar.nrays];
 614                 totden = 0.0;
 615                 totlen = 0.0;
 616                 hfside = (REAL) 0.5 * GeoPar.nelem * GeoPar.pixsiz;
 617                 for (np = 0; np < GeoPar.prjnum; ++np)
 618                 {
 619                         ProjFile.ReadProj(np, projdata, GeoPar.nrays);
 620                         for (nr = 0; nr < GeoPar.nrays; nr++)
 621                         {
 622                                 if (GeoPar.strip)
 623                                 {
 624                                         if (GeoPar.vri)
 625                                         {
 626                                                 totden +=
 627                                                                 (REAL) (projdata[nr]
 628                                                                                 / (GeoPar.pinc
 629                                                                                                 * MAX0(fabs(sinth), fabs(costh))));
 630                                         }
 631                                         else
 632                                         {
 633                                                 totden += projdata[nr] / GeoPar.pinc;
 634                                         }
 635                                 }
 636                                 else
 637                                 {
 638                                         totden += projdata[nr];
 639                                 }
 640                                 posit(np, nr, &ax, &ay, &mx, &my);
 641                                 totlen += raylen(SOT_rect, ax, ay, mx, my, 0.0, 0.0, hfside,
 642                                                 hfside, 1.0, 0.0);
 643                         }
 644                 }
 645                 GeoPar.aveden = totden / totlen;
 646                 fprintf(output, "\n         estimate of totlen = %16.6f", totlen);
 647                 fprintf(output, "\n         estimate of totden = %16.6f", totden);
 648                 fprintf(output, "\n         estimate of average density = %10.4f",
 649                                 GeoPar.aveden);
 650                 delete[] projdata;
 651
 652                 ProjFile.Close();
 653         }
 654         if (!bh_cor.bhc)
 655         {
 656                 fprintf(output, "\n         projection data read");
 657                 fprintf(output, "\n         %s", Creacm.name);
 658         }
 659         // WRITE PROJECTION HEADER INTO RECORD nelem + 3 OF
 660         // RECONSTRUCTION FILE
 661
 662         RecFile.SetProjName(Creacm.name);
 663
 664         //delete the created directory bhc
 665         if (bh_cor.bhc)
 666         {
 667                 sprintf(command, "rm -rf %s", bh_cor.tempDir);
 668                 if (system(command) != 0)
 669                 {
 670                         fprintf(stdout, "\n\n****ERROR deleting temporary directory %s\n\n",
 671                                         bh_cor.tempDir);
 672                 }
 673         }
 674
 675         if (lbase != NULL)
 676                 delete[] lbase;
 677         if (wbase != NULL)
 678                 delete[] wbase; //wei 3/2005
 679         if (BHC_PrjFile != NULL)
 680                 delete[] BHC_PrjFile;
 681
 682         for (i = 6; i >= 0; i--)
 683         {
 684                 if (phantom_poly_array[i] != NULL)
 685                         delete[] phantom_poly_array[i];
 686         }
 687
 688         return;
 689 }