src/snark/mart.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/mart.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  IMPLEMENTATIONS OF VARIOUS VERSIONS OF MULTIPLICATIVE ART ALGORITHM
  11  PROGRAMMED BY A. V. LAKSHMINARAYANAN AND A. LENT
  12  SUNY AT BUFFALO, MAY, 1977;  REVISED JAN 78
  13  IT IS ASSUMED THAT PROJECTION DATA ARE INTRINSICALLY NONNEGATIVE
  14  NONNEGATIVE RECONSTRUCTIONS ARE SOUGHT
  15
  16  DESCRIPTION OF PARAMETERS
  17
  18  INPUT CONSISTS OF THE KEYWORD 'METHOD' FOLLOWED BY ONE OF FOUR
  19  OPTIONS 'GBH', 'LENT', 'LAK1', 'LAK2', FOLLOWED BY AN INTEGER
  20  MODIFIER KOUNT, TWO FLOATING POINT MODIFIERS  RELAX AND DATZRO,
  21  OPTIONALLY FOLLOWED BY WORD MODIFIERS 'NORMAL', 'ENTROPY'
  22  THE KEYWORD 'METHOD' MUST BE PRESENT; ITS MODIFIERS HAVE THE
  23  FOLLOWING MEANING
  24
  25  'GBH', USE MULTIPLICATIVE ART ALGORITHM OF GORDON AND HERMAN
  26  INT. REV. CYTOL. 38, 1974
  27  START WITH PICTURE UNIFORMLY AVEDEN
  28  RELAXATION PARAMETER ADDED TO LITERATURE VERSION OF
  29  ALGORITHM
  30  'LENT', USE MULTIPLICATIVE ART ALGORITHM OF LENT
  31  PROC. CONF ON IMAGE ANALYSIS AND EVALUATION, JULY,1976, TORON
  32  SOCIETY OF PHOTOGRAPHIC SCIENTISTS AND ENGINEERS
  33  START WITH PICTURE AVEDEN EXCEPT FOR PIXELS IN RAYS WITH ZERO
  34  PROJECTION DATA
  35  THIS STARTING POINT IS IMMEDIATELY CHANGED BY APPLICATION
  36  OF NORMALIZATION CONSTRAINT
  37  NORMALIZE COEFFICIENT MATRIX BY ROWS
  38
  39  'LAK1', USE (APPROX) MULTIPLICATIVE ART ALGORITHM OF LAKSHMINARAYANAN
  40  LINEAR APPROXIMATION TO THE NONINTEGRAL POWER OF MULT ART.
  41  FIRST TWO TERMS OF THE BINOMIAL EXPANSION USED.
  42  USE SAME STARTING POINT AS UNDER 'LENT'
  43  NORMALIZE COEFFICIENTS BY MAX ELEMENT OF PROJECTION MATRIX
  44
  45  'LAK2', USE (APPROX) MULT ART OF LAKSHMINARAYANAN
  46  QUADRATIC APPROXIMATION TO THE NONINTEGRAL POWER OF MULT ART.
  47  FIRST THREE TERMS OF THE BINOMIAL EXPANSION USED.
  48  SAME STARTING POINT AND COEFFICIENT NORMALIZATION AS 'LAK1'
  49
  50  IF KOUNT .NE. 0, USE ONLY THE FIRST KOUNT RAYS OF THE PICK SEQUENCE
  51  THE EFFECT OF KOUNT IS DETERMINED BY THE SNARK SELECT CARD
  52  DEFAULT...IF KOUNT .LE. 0, KOUNT RESET  TO COMPLETE SET OF PROJ DATA
  53
  54  RELAX IS THE UNDERRELXATION PARAMETER.
  55  VALUES LESS THAN 1.0 ARE SUGGESTED FOR NOISY PROJECTION DATA.
  56  VALUES BIGGER THAN 1.0 ARE NOT PROHIBITED BY THEORY AND MAY SPEED
  57  CONVERGENCE.
  58  THE    PRECEDING REMARKS ARE NOT MEANINGFUL FOR 'GBH'
  59  FOR 'GBH', RELAX MUST BE CONSIDERED TO HAVE DIMENSIONS
  60
  61  DATZRO..ANY PROJECTION DATA WITH VALUE .LE. THAN DATZRO WILL BE
  62  CONSIDERED ZERO BY THE PROGRAM.
  63
  64  'NORMAL' , IF PRESENT, NORMALIZE RECON TO DENSITY AVEDEN
  65  AT END OF EACH ITERATION
  66
  67  'ENTROPY' ,IF PRESENT, CALCULATE AND PRINT ENTROPY FUNCTIONAL
  68  DEFINITION
  69  LET T = SIGMA RECON(I) OVER ALL PIXELS
  70  LET Y(I) = RECON(I)/T
  71  ENTROPY = - SIGMA Y(I)*ALOG(Y(I))
  72
  73  EXAMPLE INPUT SEQUENCES FOLLOW
  74  METHOD GBH 0 RELAX = 1.0 IGNORE RAYSUMS LT 0.0
  75  METHOD LAK1 KOUNT = 0  0.8  0.0 ENTROPY
  76  */
  77
  78 #include <cstdio>
  79 #include <cmath>
  80
  81 #include "blkdta.h"
  82 #include "geom.h"
  83 #include "consts.h"
  84 #include "uiod.h"
  85 #include "int2str.h"
  86 #include "pick.h"
  87 #include "ray.h"
  88 #include "wray.h"
  89 #include "anglst.h"
  90 #include "pseudo.h"
  91 #include "infile.h"
  92 #include "blob.h"
  93 #include "raysel.h"
  94
  95 #include "projfile.h"  //for debugg purpose
  96
  97 #include "mart.h"
  98
  99 BOOLEAN mart_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
 100 {
 101         static const INTEGER hmeth = CHAR2INT('m', 'e', 't', 'h');
 102         static const INTEGER hnorm = CHAR2INT('n', 'o', 'r', 'm');
 103         static const INTEGER hentr = CHAR2INT('e', 'n', 't', 'r');
 104
 105         static const INTEGER gbh = CHAR2INT('g', 'b', 'h', ' ');
 106         static const INTEGER lent = CHAR2INT('l', 'e', 'n', 't');
 107         static const INTEGER lak1 = CHAR2INT('l', 'a', 'k', '1');
 108         static const INTEGER lak2 = CHAR2INT('l', 'a', 'k', '2');
 109
 110         BOOLEAN xalg0_tmp;
 111         BOOLEAN eol;
 112         static INTEGER nalg;
 113         INTEGER word;
 114         static REAL datzro; //set static
 115         static INTEGER lnorm;
 116         static INTEGER lfnal;
 117         INTEGER i;
 118         INTEGER np;
 119         INTEGER nr;
 120         REAL raysum;
 121         INTEGER numb;
 122         REAL snorm;
 123         INTEGER j;
 124         INTEGER l;
 125         INTEGER isum;
 126         REAL fact;
 127         static REAL aijmax; // set static
 128         REAL wj;
 129         INTEGER k;
 130         REAL trysum;
 131         REAL corfac;
 132         REAL delu;
 133         REAL ajmax;
 134         REAL corr;
 135         REAL roaij;
 136         REAL reldif;
 137         REAL aij;
 138         REAL hafrel;
 139         REAL hafabs;
 140         REAL fac;
 141         REAL tot;
 142         REAL entrop;
 143         REAL rot;
 144         static INTEGER kount;
 145         static REAL relax; //set static
 146
 147         INTEGER area;
 148
 149         if (Blob.pix_basis)
 150         {
 151                 area = GeoPar.area;
 152                 //    num_rays = GeoPar.nrays;
 153         }
 154
 155         else
 156         {
 157                 area = Blob.area;
 158                 //    num_rays = GeoPar.usrays;
 159         }
 160
 161         xalg0_tmp = FALSE;
 162
 163         if (!(iter > 1))
 164         {
 165
 166                 // GET THE METHOD
 167
 168                 INTEGER exp0[1] =
 169                 { hmeth };
 170                 word = InFile.getwrd(TRUE, &eol, exp0, 1);
 171
 172                 if (eol)
 173                 {
 174                         fprintf(output,
 175                                         "\n          **** error - keyword method is missing***");
 176                         return TRUE;
 177                 }
 178
 179                 INTEGER exp1[4] =
 180                 { gbh, lent, lak1, lak2 };
 181
 182                 word = InFile.getwrd(FALSE, &eol, exp1, 4);
 183
 184                 if (eol)
 185                 {
 186                         fprintf(output, "\n          **** error - one and only one of the following methods must be specified : gbh, lent, lak1, or lak2***");
 187                         return TRUE;
 188                 }
 189
 190                 if (word == gbh)
 191                 {
 192                         nalg = 1;
 193                         fprintf(output, "\n          multiplicative art version %s",
 194                                         int2str(gbh));
 195                 }
 196                 else if (word == lent)
 197                 {
 198                         nalg = 2;
 199                         fprintf(output, "\n          multiplicative art version %s",
 200                                         int2str(lent));
 201                 }
 202                 else if (word == lak1)
 203                 {
 204                         nalg = 3;
 205                         fprintf(output, "\n          multiplicative art version %s",
 206                                         int2str(lak1));
 207                 }
 208                 else
 209                 {  //lak2
 210                         nalg = 4;
 211                         fprintf(output, "\n          multiplicative art version %s",
 212                                         int2str(lak2));
 213                 }
 214
 215                 // GET KOUNT, RELAXATION PARAMETER, RAYSUM LIMIT, NORMAL, ENTROPY
 216                 kount = InFile.getint(FALSE, &eol);
 217
 218                 relax = InFile.getnum(FALSE, &eol);
 219
 220                 datzro = InFile.getnum(FALSE, &eol);
 221
 222                 if (eol)
 223                 {
 224                         fprintf(output,
 225                                         "\n          **** some data is missing: either kount, relax, or data_zero***");
 226                         return TRUE;
 227                 }
 228
 229                 // optional inputs
 230
 231                 INTEGER exp3[2] =
 232                 { hnorm, hentr };
 233                 word = InFile.getwrd(FALSE, &eol, exp3, 2);
 234
 235                 if (eol)
 236                         goto L10;
 237                 // finished reading input
 238
 239                 lnorm = FALSE;
 240                 lfnal = FALSE;
 241
 242                 if (word == hnorm)
 243                 {
 244                         lnorm = TRUE;
 245                         INTEGER exp3a[1] =
 246                         { hentr };
 247                         word = InFile.getwrd(FALSE, &eol, exp3a, 1);
 248                         if (word == hentr)
 249                                 goto L5;
 250                         else
 251                                 goto L10;
 252                 }
 253
 254                 else
 255                 {
 256                         L5: lfnal = TRUE;
 257                 }
 258
 259                 //***  CHECKING INPUTS AND SETTING DEFAULTS
 260                 L10: if (kount <= 0)
 261                 {
 262                         if (RaySel.useray)  // select user
 263                                 kount = GeoPar.prjnum * GeoPar.usrays; // change nrays to usrays. hstau 7/2003
 264                         else
 265                                 //select snark
 266                                 kount = GeoPar.prjnum * GeoPar.snrays;
 267                 } // set value to kount for select user or select snark. Wei   11/2004
 268
 269                 if (relax <= 0)
 270                         relax = 1.0;
 271                 if ((nalg == 3) && (relax > 1.0))
 272                         relax = 1.0;
 273                 if (datzro < Consts.zero)
 274                         datzro = Consts.zero;
 275
 276                 fprintf(output, "\n          no of rays used in each iteration %10i",
 277                                 kount);
 278                 fprintf(output, "\n          underrelaxation factor %10.5f", relax);
 279                 fprintf(output,
 280                                 "\n          all projection data with values less than %20.10e are ignored",
 281                                 datzro);
 282
 283                 if (lnorm)
 284                 {
 285                         fprintf(output, "\n          reconstruction is normalized");
 286                 }
 287                 else
 288                 {
 289                         fprintf(output, "\n          reconstruction not normalized");
 290                 }
 291
 292                 if (lfnal)
 293                 {
 294                         fprintf(output, "\n          entropy functional is calculated");
 295                 }
 296                 else
 297                 {
 298                         fprintf(output, "\n          entropy functional not calculated");
 299                 }
 300
 301                 if (GeoPar.aveden <= Consts.zero)
 302                 {
 303                         fprintf(output,
 304                                         "\n ***aveden too small.  probably incorrect use of mart, no execution***");
 305                         return TRUE;
 306                 }
 307
 308                 for (i = 0; i < area; i++)
 309                 {
 310                         recon[i] = GeoPar.aveden;
 311                 }
 312
 313                 //***  NO CHOICE IN INITIAL GUESS; EXALG OPTIONS DONT OPERATE
 314
 315                 if (nalg != 1)
 316                 {
 317                         //  WILL ZERO OUT ANY PIXEL IN A RAY WITH ZERO PROJ DATA
 318                         for (np = 0; np < GeoPar.prjnum; np++)
 319                         {
 320                                 for (nr = 0; nr < GeoPar.nrays; nr++)
 321                                 {
 322                                         raysum = Anglst.prdta(np, nr);
 323                                         if (raysum <= datzro)
 324                                         {
 325                                                 if (GeoPar.strip)
 326                                                 {
 327                                                         if (Blob.pix_basis)
 328                                                         {
 329                                                                 ray(np, nr, list, weight, &numb, &snorm);
 330                                                         }
 331                                                         else
 332                                                         {
 333                                                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 334                                                         }
 335                                                 }
 336                                                 if (GeoPar.line)
 337                                                 {
 338                                                         if (Blob.pix_basis)
 339                                                         {
 340                                                                 wray(np, nr, list, weight, &numb, &snorm);
 341                                                         }
 342                                                         else
 343                                                         {
 344                                                                 Blob.bwray(np, nr, list, weight, &numb, &snorm);
 345                                                         }
 346                                                 }
 347                                                 if (numb != 0)
 348                                                 {
 349                                                         for (j = 0; j < numb; j++)
 350                                                         {
 351                                                                 l = list[j];
 352                                                                 recon[l] = 0.0;
 353                                                         }
 354                                                 }
 355                                         }
 356                                 }
 357
 358                         }
 359
 360                         isum = 0;
 361
 362                         for (i = 0; i < area; i++)
 363                         {
 364                                 if (recon[i] != 0.0)
 365                                 {
 366
 367                                         if (!Blob.pix_basis)
 368                                         { //if blobs
 369                                                 if (i % 2 == Blob.pr)
 370                                                 {
 371                                                         isum++;
 372                                                 }
 373                                         }
 374                                         else
 375                                         {
 376                                                 isum++;
 377                                         }
 378                                 }
 379                         }
 380
 381                         if (isum == 0)
 382                         {
 383
 384                                 fprintf(output,
 385                                                 "\n  ***warning-  all projection data near zero. zero reconstruction.***");
 386                                 return TRUE;
 387                         }
 388
 389                         //  NORMALIZE
 390
 391                         fact = area * GeoPar.aveden / isum;
 392                         if (!Blob.pix_basis)
 393                         { //if blobs
 394                                 fact = fact / 2.0;
 395                         }
 396                         for (i = 0; i < area; i++)
 397                         {
 398                                 if (recon[i] != 0.0)
 399                                         recon[i] = fact;
 400                         }
 401
 402                         //  LOCATE BIGGEST ELEMENT IN PROJECTION MATRIX FOR LAK1 AND LAK2 ONLY
 403                         //  THIS STRATEGY IS CONSERVATIVE IN THAT COEFFICIENTS OF PIXELS
 404                         //  WHOSE VALUES ARE KNOWN TO BE ZERO ARE INCLUDED IN THE SEARCH
 405                         if (nalg != 2)
 406                         {
 407                                 aijmax = 0.0;
 408                                 for (np = 0; np < GeoPar.prjnum; np++)
 409                                 {
 410                                         for (nr = 0; nr < GeoPar.nrays; nr++)
 411                                         {
 412                                                 if (GeoPar.strip)
 413                                                 {
 414                                                         if (Blob.pix_basis)
 415                                                         {
 416                                                                 ray(np, nr, list, weight, &numb, &snorm);
 417                                                         }
 418                                                         else
 419                                                         {
 420                                                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 421                                                         }
 422                                                 }
 423                                                 if (GeoPar.line)
 424                                                 {
 425                                                         if (Blob.pix_basis)
 426                                                         {
 427                                                                 wray(np, nr, list, weight, &numb, &snorm);
 428                                                         }
 429                                                         else
 430                                                         {
 431                                                                 Blob.bwray(np, nr, list, weight, &numb, &snorm);
 432                                                         }
 433                                                 }
 434                                                 if (numb != 0)
 435                                                 {
 436                                                         for (j = 0; j < numb; j++)
 437                                                         {
 438                                                                 wj = weight[j];
 439                                                                 if (wj > aijmax)
 440                                                                         aijmax = wj;
 441                                                         }
 442                                                 }
 443                                         }
 444
 445                                 }
 446
 447                                 //  BASIC UPDATE PART OF PROGRAM
 448                         }
 449                 }
 450         }
 451
 452         switch (nalg)
 453         {
 454
 455         //  GORDON,BENDER,HERMAN MULT ART
 456         case 1:
 457                 for (k = 0; k < kount; k++)
 458                 {
 459                         pick(&np, &nr);
 460
 461                         raysum = Anglst.prdta(np, nr);
 462                         if (Blob.pix_basis)
 463                         {
 464                                 trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
 465                                                 GeoPar.line, FALSE);
 466                         }
 467                         else
 468                         {
 469                                 trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
 470                                                 &snorm, GeoPar.line, FALSE);
 471                         }
 472                         if (numb != 0)
 473                         {
 474                                 if (raysum <= datzro)
 475                                 {
 476                                         for (j = 0; j < numb; j++)
 477                                         {
 478                                                 l = list[j];
 479                                                 recon[l] = 0.0;
 480                                         }
 481                                 }
 482                                 else
 483                                 {
 484
 485                                         if ((trysum <= Consts.zero) && (trace >= 1))
 486                                         {
 487                                                 fprintf(output,
 488                                                                 "\n zero pseudo proj dat %5i %5i %10.5f %10.5f",
 489                                                                 np, nr, raysum, trysum);
 490                                         }
 491
 492                                         if (trysum > Consts.zero)
 493                                         {
 494                                                 corfac = relax * (REAL) log(raysum / trysum);
 495                                                 for (j = 0; j < numb; j++)
 496                                                 {
 497                                                         l = list[j];
 498                                                         recon[l] *= (REAL) exp(corfac * weight[j]);
 499                                                 }
 500                                         }
 501                                 }
 502                         }
 503                 }
 504
 505                 //  END GBH UBDATE
 506                 break;
 507
 508                 //  LENT MULT ART
 509
 510         case 2:
 511
 512                 for (k = 0; k < kount; k++)
 513                 {
 514                         pick(&np, &nr);
 515                         raysum = Anglst.prdta(np, nr);
 516                         if (raysum > datzro)
 517                         {
 518                                 if (Blob.pix_basis)
 519                                 {
 520                                         trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
 521                                                         GeoPar.line, FALSE);
 522                                 }
 523                                 else
 524                                 {
 525                                         trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
 526                                                         &snorm, GeoPar.line, FALSE);
 527                                 }
 528                                 if (numb != 0)
 529                                 {
 530                                         if (trysum != 0.0)
 531                                         {
 532                                                 delu = (REAL) log(raysum / trysum);
 533
 534                                                 // FIND ROW SCALE FACTOR
 535                                                 ajmax = 0.0;
 536
 537                                                 for (j = 0; j < numb; j++)
 538                                                 {
 539                                                         if (weight[j] > ajmax)
 540                                                                 ajmax = weight[j];
 541                                                 }
 542
 543                                                 corr = relax * delu / ajmax;
 544
 545                                                 for (j = 0; j < numb; j++)
 546                                                 {
 547                                                         l = list[j];
 548                                                         recon[l] *= (REAL) exp(corr * weight[j]);
 549                                                 }
 550
 551                                         }
 552                                 }
 553                         }
 554                 }
 555
 556                 // END LENT UPDATE
 557                 break;
 558
 559                 // LAK APPROX MULT ART - LINEAR APPROXIMATION
 560         case 3:
 561
 562                 roaij = relax / aijmax;
 563
 564                 for (k = 0; k < kount; k++)
 565                 {
 566                         pick(&np, &nr);
 567                         raysum = Anglst.prdta(np, nr);
 568                         if (raysum > datzro)
 569                         {
 570                                 if (Blob.pix_basis)
 571                                 {
 572                                         trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
 573                                                         GeoPar.line, FALSE);
 574                                 }
 575                                 else
 576                                 {
 577                                         trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
 578                                                         &snorm, GeoPar.line, FALSE);
 579                                 }
 580                                 if (!((numb == 0) || (trysum <= Consts.zero)))
 581                                 {
 582                                         reldif = (raysum - trysum) / MAX0(raysum, trysum);
 583
 584                                         // LINEAR APPROX
 585                                         for (j = 0; j < numb; j++)
 586                                         {
 587                                                 l = list[j];
 588                                                 aij = weight[j] * roaij;
 589                                                 recon[l] *= ((REAL) 1.0 + reldif * aij);
 590                                         }
 591                                 }
 592                         }
 593                 }
 594                 break;
 595
 596                 //END LAK1 UPDATE
 597                 // LAK APPROX MULT ART - QUADRATIC APPROXIMATION
 598
 599         case 4:
 600
 601                 roaij = relax / aijmax;
 602
 603                 for (k = 0; k < kount; k++)
 604                 {
 605                         pick(&np, &nr);
 606                         raysum = Anglst.prdta(np, nr);
 607                         if (raysum > datzro)
 608                         {
 609                                 if (Blob.pix_basis)
 610                                 {
 611                                         trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
 612                                                         GeoPar.line, FALSE);
 613                                 }
 614                                 else
 615                                 {
 616                                         trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
 617                                                         &snorm, GeoPar.line, FALSE);
 618                                 }
 619                                 if (!((numb == 0) || (trysum <= Consts.zero)))
 620                                 {
 621                                         reldif = (raysum - trysum) / MAX0(raysum, trysum);
 622
 623                                         // QUADRATIC APPROX
 624                                         hafrel = (REAL) 0.5 * reldif;
 625                                         hafabs = (REAL) fabs(hafrel);
 626                                         for (j = 0; j < numb; j++)
 627                                         {
 628                                                 l = list[j];
 629                                                 aij = weight[j] * roaij;
 630                                                 fac = (REAL) 1.0
 631                                                                 + reldif * aij
 632                                                                                 * ((REAL) 1.0 + hafrel * aij + hafabs);
 633                                                 recon[l] *= fac;
 634                                         }
 635                                 }
 636                         }
 637                 }
 638                 break;
 639
 640         }
 641         //  END LAK2 UPDATE
 642
 643         //  CHECK OPTIONS
 644
 645         if (!(lnorm || lfnal))
 646                 return xalg0_tmp;
 647
 648         tot = 0.0;
 649
 650         for (i = 0; i < area; i++)
 651         {
 652                 if (!Blob.pix_basis)
 653                 { //if blobs
 654                         if (i % 2 == Blob.pr)
 655                         {
 656                                 tot += recon[i];
 657                         }
 658                 }
 659                 else
 660                 {
 661                         tot += recon[i];
 662                 }
 663         }
 664
 665         if (tot > Consts.zero)
 666         {
 667                 if (lnorm)
 668                 {
 669                         fact = GeoPar.aveden * area / tot;  // aveden in both basis concide
 670                         if (!Blob.pix_basis)
 671                         {
 672                                 fact /= 2.0;
 673                         }
 674                         for (i = 0; i < area; i++)
 675                         {
 676                                 recon[i] *= fact;
 677                         }
 678
 679                         tot = GeoPar.aveden * area;
 680                         if (!Blob.pix_basis)
 681                         {
 682                                 tot /= 2.0;
 683                         }
 684                 }
 685
 686                 if (!lfnal)
 687                         return xalg0_tmp;
 688
 689                 entrop = 0.0;
 690                 for (i = 0; i < area; i++)
 691                 {
 692                         if (recon[i] > 0.0)
 693                         {
 694                                 rot = recon[i] / tot;
 695                                 entrop -= rot * (REAL) log(rot); //remove the minus sign here, wei
 696                         }
 697                 }
 698
 699                 fprintf(output,
 700                                 "\n                                                        iter %5i entropy %20.10e",
 701                                 iter, entrop);
 702
 703                 return xalg0_tmp;
 704         }
 705         fprintf(output,
 706                         "\n ***total density %20.10e  is too small, no normalization or entropy calculation***",
 707                         tot);
 708         return xalg0_tmp;
 709 }