src/snark/sirt.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/sirt.cpp $
   5  $LastChangedRevision: 97 $
   6  $Date: 2014-07-02 20:07:41 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  IMPLEMENTATION OF THE SIRT ALGORITHM OF PETER GILBERT
  11  TOGETHER WITH EXTENSIONS AND MODIFICATIONS BY
  12  LAKSHMINARAYANAN AND LENT OF SUNY BUFFALO
  13  MAY, 1977;  REVISED JANUARY, 1978.
  14  SEE TECHNICAL REPORT 95 OF DEPT OF COMPUTER SCIENCE
  15  PROGRAMMED IN ANSI FORTRAN FOR PORTABILITY
  16  USES AUXILIARY FILES USER1 AND USER2 TO REDUCE CORE
  17
  18  METHOD...CONTROLS ALGORITHM SELECTION
  19  IF 'GSIRT' USE LAK-LENT GENERALIZATION OF GILBERT SIRT
  20  SEE EQUATION (20) OF REPORT 95
  21  RHO(Q+1) = RHO(Q) + 1.0/SIGMA *(B-(DN)**-1 * PTRANS * P *
  22  RHO(Q)
  23  WHERE
  24  B = (DL)**-1 * P TRANS * R
  25  DN = P TRANS * P * U, U A VECTOR OF ONES
  26  DL IS A WEIGHTED SUM OF RAY LENGTHS
  27  R IS VECTOR OF PROJECTION DATA
  28
  29  FOR 'LSIRT', MODIFIED SIRT ALGORITHMS OF LAK AND LENT
  30  ARE SELECTED
  31  GENERAL FORM OF MODIFIED SIRT IS (SEE EQUATION (23) )
  32  RHO(Q+1) = RHO(Q) + 1.0/SIGMA * (B - DV**-1 * PTRANS *
  33  Z * V**2 * P * RHO(Q))
  34  DEFINITIONS FOLLOW
  35  N(J) = NUMBER OF PIXELS IN RAY J
  36  L(J) = LENGTH OF RAY J
  37  W(J) = WIDTH OF RAY J
  38  FOR THE LSIRT ALGORITHMS, (Z * V**2)(J) TAKES THE FOLLOWING FORM
  39  N(J)**POWER2
  40  DV(I) TAKES THE FOLLOWING FORM
  41  SUM OVER ALL RAYS J PASSING THROUGH PIXEL I OF N(J)**POWER1
  42  *** DV IS DESIGNATED THROUGHOUT THIS PROGRAM BY THE SYMBOL DN ***
  43  B (IN THE TECHNICAL REPORT) HAS THE FOLLOWING FORM
  44  DV**-1 * PTRANS * Z * V * R
  45  FOR EACH RAY J
  46  (Z * V * R)(J) = R(J)/(L(J)*W(J)) * N(J)**POWER1
  47  VALUES OF LSIRT MODIFIER
  48  LSIRT = 1       POWER1=1            POWER2=0
  49  LSIRT = 2       POWER1 = 0          POWER2 = -1
  50  LSIRT = 3       POWER1 = -1         POWER2 = -2
  51  TO GET CONSTRAINED SIRT, USE THE SNARK MODE CARD
  52
  53  ***INPUT SEQUENCE
  54
  55  THE FOLLOWING SEQUENCE OF KEYWORDS IS AVAILABLE TO THE USER
  56
  57  METHOD
  58  THIS KEYWORD MUST BE PRESENT
  59  IT MAY BE FOLLOWED BY EITHER 'GSIRT' OR 'LSIRT TYPE'.
  60  TYPE IS AN INTEGER MODIFIER WHICH MUST TAKE ON VALUES
  61  1, 2, OR 3.
  62
  63  RELAX
  64  IF THIS IS PRESENT, IT MUST BE FOLLOWED BY THE FLOATING-
  65  POINT VALUE OF THE RELAXATION PARAMETER;
  66  SIGMA
  67  IF THIS IS PRESENT, IT MUST BE FOLLOWED BY THE FLOATING-
  68  POINT VALUE OF SIGMA;
  69  RELAXATION FACTOR = 1.0/SIGMA
  70
  71  *****  RELAX AND SIGMA MUST NOT BOTH BE USED  *****
  72
  73  IF NEITHER RELAX OR SIGMA IS USED, RELAXATION FACTOR =
  74  SIGMA = 1.0
  75
  76  START
  77  MAY BE FOLLOWED BY 'BACKPROJECTED'
  78  THIS CAUSES A BACKPROJECTED IMAGE TO BE USED AS THE
  79  INITIAL GUESS
  80  IF START IS NOT PRESENT, THE EXEC COMMAND OPTION IS USED
  81
  82  NORMAL
  83  THIS CAUSES ALL RECONSTRUCTIONS (OTHER THAN THE INITIAL GUESS)
  84  TO BE NORMALIZED TO AVEDEN.  THIS STEP IS OMITTED IF
  85  AVEDEN IS TOO SMALL OR IF THE TOTAL DENSITY OF THE
  86  ESTIMATED IMAGE IS LOW
  87
  88  EXAMPLE INPUT SEQUENCES FOLLOW
  89
  90  METHOD GSIRT SIGMA = .5 START NORMAL
  91  METHOD LSIRT 1 RELAX = 3.
  92
  93  **** THIS SIRT IS FOR STRIP PROJECTION DATA ONLY  ****
  94
  95  ADDED FORMAT 1090 FOR WRITING INTO USER1 BAND USER2..7/21/88.SADA
  96  */
  97
  98 #include <cstdio>
  99 #include <cmath>
 100
 101 #include "blkdta.h"
 102 #include "geom.h"
 103 #include "spctrm.h"
 104 #include "fourie.h"
 105 #include "raysel.h"
 106 #include "modefl.h"
 107 #include "consts.h"
 108 #include "uiod.h"
 109 #include "raylen.h"
 110 #include "posit.h"
 111 #include "anglst.h"
 112 #include "ray.h"
 113 #include "projfile.h"
 114 #include "infile.h"
 115 #include "blob.h"
 116
 117 #include "sirt.h"
 118
 119 #define  nrdis(rinc) MIN0(((INTEGER)(size * (fabs(costh) + fabs(sinth))/(rinc) - 0.001)), limit)
 120
 121 INTEGER sirt_class::Init()
 122 {
 123
 124         hmeth = CHAR2INT('m', 'e', 't', 'h');
 125         hgsir = CHAR2INT('g', 's', 'i', 'r');
 126         hlsir = CHAR2INT('l', 's', 'i', 'r');
 127         hstar = CHAR2INT('s', 't', 'a', 'r');
 128         hnorm = CHAR2INT('n', 'o', 'r', 'm');
 129         hrela = CHAR2INT('r', 'e', 'l', 'a');
 130         hsigm = CHAR2INT('s', 'i', 'g', 'm');
 131         hback = CHAR2INT('b', 'a', 'c', 'k');
 132
 133         w1 = NULL;
 134         w2 = NULL;
 135         w3 = NULL;
 136         gen_inv_dn = NULL;
 137         vector_b = NULL;
 138         pbase = NULL;
 139         sbase = NULL;             //wei, 3/2005
 140
 141         return 0;
 142 }
 143
 144 INTEGER sirt_class::Reset()
 145 {
 146         if (w1 != NULL)
 147                 delete[] w1;
 148         if (w2 != NULL)
 149                 delete[] w2;
 150         if (w3 != NULL)
 151                 delete[] w3;
 152         if (gen_inv_dn != NULL)
 153                 delete[] gen_inv_dn;
 154         if (vector_b != NULL)
 155                 delete[] vector_b;
 156         if (pbase != NULL)
 157                 delete[] pbase;
 158         if (sbase != NULL)
 159                 delete[] sbase;               //wei, 3/2005
 160
 161         return 0;
 162 }
 163
 164 BOOLEAN sirt_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
 165 {
 166
 167         static REAL totden;
 168         REAL size;
 169         INTEGER i;
 170         INTEGER np;
 171         REAL theta;
 172         REAL sinth;
 173         REAL costh;
 174         INTEGER nr;
 175         INTEGER numb;
 176         REAL snorm;
 177         REAL ax;
 178         REAL ay;
 179         REAL cth;
 180         REAL sth;
 181         REAL wl;
 182         REAL raynar;
 183         REAL raylar;
 184         REAL dl;
 185         INTEGER nb;
 186         INTEGER k;
 187         REAL raysum;
 188         REAL w;
 189         REAL psum;
 190         REAL ronew;
 191         REAL ratio;
 192
 193         static BOOLEAN normal;
 194         BOOLEAN bkstrt, gsirt, sirt, eol;
 195
 196         INTEGER limit, first, last;
 197         INTEGER word;
 198         static INTEGER lsirt; // set static hstau
 199         REAL relax;
 200         static REAL sigma;
 201
 202         INTEGER area;
 203
 204         if (Blob.pix_basis)
 205         {
 206                 area = GeoPar.area;
 207         }
 208         else
 209         {
 210                 area = Blob.area;
 211         }
 212
 213         // STATEMENT FUNCTION ENABLING TO DISCARD RAYS THAT DONT
 214         // INTERSECT THE RECONSTRUCTION REGION IN THE SUMMATIONS
 215         // DESCRIBED IN THE PREFACE
 216
 217         if (iter <= 1)
 218         {
 219
 220                 //major changes in the input reading sequence, using new getwrd funcion. hstau 1/2004
 221                 lsirt = 0;
 222                 //for stripe mode only!
 223                 if (GeoPar.line)
 224                 {
 225                         fprintf(output, "\n *** this sirt algorithm is for strip data only, try quad algorithm for line data**");
 226                         return TRUE;
 227                 }
 228
 229                 // SELECT THE SNARK OR USER SET OF RAYS
 230                 limit = GeoPar.lusray - GeoPar.midray;
 231
 232                 if (GeoPar.snrays != GeoPar.usrays)
 233                 {
 234
 235                         //  CHECK FOR SNARK OR USER
 236                         if (!RaySel.useray)
 237                         {
 238
 239                                 // SNARK OPTION IN EFFECT
 240
 241                                 fprintf(output, "\n          snark rays selected");
 242
 243                                 limit = GeoPar.lsnray - GeoPar.midray;
 244                         }
 245                         else
 246                         {
 247                                 // USER OPTION IN EFFECT
 248                                 fprintf(output, "\n          user rays selected");
 249                         }
 250                 }
 251
 252                 //input phase
 253                 INTEGER exp0[1] =
 254                 { hmeth };
 255                 word = InFile.getwrd(TRUE, &eol, exp0, 1);
 256
 257                 if (eol)
 258                 {
 259                         fprintf(output, "\n          **** error - keyword method is missing***");
 260                         return TRUE;
 261                 }
 262
 263                 INTEGER exp1[2] =
 264                 { hgsir, hlsir };
 265                 word = InFile.getwrd(FALSE, &eol, exp1, 2);
 266
 267                 if (eol)
 268                 {
 269                         fprintf(output, "\n          **** error - one and only one of the following methods must be specified: gsirt or lsirt***");
 270                         return TRUE;
 271                 }
 272
 273                 gsirt = word == hgsir;
 274                 sirt = word == hlsir;
 275
 276                 if (word == hgsir)
 277                 {
 278                         fprintf(output, "\n          gsirt method");
 279                 }
 280                 else
 281                 {  //hlsir
 282                         lsirt = InFile.getint(FALSE, &eol);
 283                         if (eol)
 284                         {
 285                                 fprintf(output, "\n          **** error - must specify lsirt type ***");
 286                                 return TRUE;
 287                         }
 288                         if (!(lsirt > 0 && lsirt < 4))
 289                         {
 290                                 return TRUE;
 291                         }
 292                         fprintf(output, "\n          lsirt method");
 293                 }
 294
 295                 //default values
 296                 sigma = 1.0;
 297                 relax = 1.0;
 298                 bkstrt = FALSE;
 299                 normal = FALSE;
 300
 301                 // optional inputs
 302                 INTEGER exp2[4] =
 303                 { hrela, hsigm, hstar, hnorm };
 304                 word = InFile.getwrd(FALSE, &eol, exp2, 4);
 305
 306                 if (eol)
 307                         goto L30;
 308                 // finished reading input
 309
 310                 if (word == hrela)
 311                 {
 312                         relax = InFile.getnum(FALSE, &eol);
 313                         if (eol)
 314                         {
 315                                 fprintf(output, "\n          **** error - must specify relax ***");
 316                                 return TRUE;
 317                         }
 318                         if (relax < Consts.zero)
 319                         {
 320                                 fprintf(output, "\n          **** error - relax too small or negative ***");
 321                                 return TRUE;
 322                         }
 323                         sigma = (REAL) 1.0 / relax;
 324
 325                         INTEGER exp2a[2] =
 326                         { hstar, hnorm };
 327                         word = InFile.getwrd(FALSE, &eol, exp2a, 2);
 328                         if (word == hstar)
 329                                 goto L10;
 330                         else if (word == hnorm)
 331                                 goto L20;
 332                         else
 333                                 goto L30;
 334                         // eol: finished
 335                 }
 336                 else if (word == hsigm)
 337                 {
 338                         sigma = InFile.getnum(FALSE, &eol);
 339                         if (eol)
 340                         {
 341                                 fprintf(output, "\n          **** error - must specify sigma ***");
 342                                 return TRUE;
 343                         }
 344                         if (sigma < Consts.zero)
 345                         {
 346                                 fprintf(output, "\n          **** error - sigma too small or negative ***");
 347                                 return TRUE;
 348                         }
 349                         INTEGER exp2b[2] =
 350                         { hstar, hnorm };
 351                         word = InFile.getwrd(FALSE, &eol, exp2b, 2);
 352                         if (word == hstar)
 353                                 goto L10;
 354                         else if (word == hnorm)
 355                                 goto L20;
 356                         else
 357                                 goto L30;
 358                         // eol: finished
 359                 }
 360
 361                 else if (word == hstar)
 362                 {
 363                         L10: bkstrt = TRUE;
 364                         INTEGER exp3[2] =
 365                         { hback, hnorm };
 366                         word = InFile.getwrd(FALSE, &eol, exp3, 2);
 367                         if (eol || word == hback)
 368                                 goto L30;
 369                         else
 370                                 goto L20;
 371                         //normal
 372                 }
 373                 else
 374                 {
 375                         L20: normal = TRUE;
 376                 }
 377
 378                 L30:
 379                 // ECHOING INPUT
 380
 381                 fprintf(output, "\n          sigma = inverse of relaxation = %10.4f\n",
 382                                 sigma);
 383                 // removed from the echoing hstau.2003
 384
 385                 // NONCONVERGENCE WARNING ON SIGMA
 386                 if (sigma < .4999)
 387                 {
 388                         fprintf(output, "\n ***small sigma value.  probable divergence.***"); //hstau
 389                 }
 390
 391                 if (Modefl.lofl || Modefl.upfl)
 392                 {
 393                         fprintf(output, "\n          this run is constrained sirt");
 394                 }
 395                 else
 396                 {
 397                         fprintf(output, "\n          this run is unconstrained sirt");
 398                 }
 399
 400                 if (normal)
 401                 {
 402                         fprintf(output, "\n          this run is normalized sirt");
 403                 }
 404                 else
 405                 {
 406                         fprintf(output, "\n          this run is unnormalized sirt");
 407                 }
 408
 409                 // CHECK ON SMALL TOTAL DENSITY AND ISSUE WARNING
 410                 totden = ((REAL) (GeoPar.area)) * GeoPar.aveden;
 411
 412                 if (fabs(totden) <= Consts.zero)
 413                 {
 414                         normal = FALSE;
 415                         fprintf(output, "\n ***average density too small; no normalization***"); //hstau
 416                 }
 417
 418                 if (bkstrt)
 419                 {
 420                         fprintf(output, "\n          starting picture is back-projected density \n");
 421                 }
 422                 else
 423                 {
 424                         fprintf(output, "\n          starting picture is as in execute command");
 425                 }
 426
 427                 // END INPUT CHECK AND ECHO
 428                 //
 429                 // *       START WORK
 430                 //
 431
 432                 w1 = new REAL[area];
 433                 w2 = new REAL[area];
 434                 w3 = new REAL[area];
 435                 gen_inv_dn = new REAL[area];
 436                 vector_b = new REAL[area];
 437
 438                 pbase = new REAL[GeoPar.nrays];
 439
 440                 sbase = new INTEGER[GeoPar.prjnum];
 441
 442                 size = GeoPar.nelem * GeoPar.pixsiz / (REAL) 2.0;
 443                 Fourie.rinc = GeoPar.pinc;
 444
 445                 // GET DN AND B ON USER1
 446                 // GENERATE DN AND DL
 447                 // IF 'LSIRT' DONT GENERATE DL,
 448                 // GET P TRANSPOSE N*R/L DIRECTLY
 449                 // WRITTEN FOR VARIABLE RAY WIDTH
 450
 451                 for (i = 0; i < area; i++)
 452                 {
 453                         w2[i] = 0.0;
 454                         w1[i] = 0.0;
 455                 }
 456
 457                 for (np = 0; np < GeoPar.prjnum; np++)
 458                 {
 459                         Anglst.getang(np, &theta, &sinth, &costh);
 460                         if (sirt)
 461                         {
 462                                 ProjFile.ReadProj(np, pbase, GeoPar.nrays);
 463                         }
 464
 465                         if (GeoPar.vri)
 466                                 Fourie.rinc = GeoPar.pinc
 467                                                 * (REAL) MAX0(fabs(costh), fabs(sinth));
 468
 469                         first = GeoPar.midray - nrdis(Fourie.rinc);
 470
 471                         sbase[np] = first;
 472                         //last = GeoPar.nrays + 1 - first;
 473                         last = GeoPar.nrays - 1 - first; //should be 2 less. hstau
 474
 475                         for (nr = first; nr <= last; nr++)
 476                         {
 477                                 if (Blob.pix_basis)
 478                                 {
 479                                         ray(np, nr, list, weight, &numb, &snorm);
 480                                 }
 481                                 else
 482                                 {
 483                                         Blob.bray(np, nr, list, weight, &numb, &snorm);
 484                                 }
 485                                 if (numb != 0)
 486                                 {
 487                                         posit(np, nr, &ax, &ay, &cth, &sth);
 488                                         wl = Fourie.rinc
 489                                                         * raylen(SOT_rect, ax, ay, cth, sth, 0., 0., size,
 490                                                                         size, 1., 0.);
 491                                         if (!sirt)
 492                                         {
 493                                                 raynar = (REAL) (numb);
 494                                                 raylar = wl;
 495                                         }
 496                                         else
 497                                         {
 498                                                 dl = pbase[nr] / wl;
 499                                                 switch (lsirt)
 500                                                 {
 501                                                 case 1:
 502                                                         raynar = float(numb);
 503                                                         raylar = dl * raynar;
 504                                                         break;
 505
 506                                                 case 2:
 507                                                         raynar = 1.0;
 508                                                         raylar = dl;
 509                                                         break;
 510
 511                                                 case 3:
 512                                                         raynar = ((REAL) 1.0) / numb;
 513                                                         raylar = dl * raynar;
 514                                                         break;
 515                                                 }
 516                                         }
 517
 518                                         for (nb = 0; nb < numb; nb++)
 519                                         {
 520                                                 k = list[nb];
 521                                                 w2[k] += raynar;
 522                                                 w1[k] += raylar;
 523                                         }
 524                                 }
 525                         }
 526
 527                 }
 528
 529                 if (!sirt)
 530                 {
 531
 532                         // GSIRT
 533                         // TRYING FOR COMPATIBILITY IN CASE SOME PIXELS DO NOT
 534                         // CONTRIBUTE TO THE PROJECTION DATA
 535
 536                         for (i = 0; i < area; i++)
 537                         {
 538                                 if (w1[i] == 0.0)
 539                                         w2[i] = 0.0;
 540                                 if (w2[i] > Consts.zero)
 541                                         w2[i] = ((REAL) 1.0) / w2[i];
 542                                 if (w2[i] == 0.0)
 543                                         w1[i] = 0.0;
 544                         }
 545
 546                         // W2 NOW HOLDS GENERALIZED INVERSE OF DN
 547                         // WORK HOLDS DL
 548
 549                         for (i = 0; i < area; i++)
 550                         {
 551                                 gen_inv_dn[i] = w2[i];
 552                         }
 553
 554                         // GET P TRANSPOSE R
 555
 556                         for (i = 0; i < area; i++)
 557                         {
 558                                 w2[i] = 0.0;
 559                         }
 560
 561                         for (np = 0; np < GeoPar.prjnum; np++)
 562                         {
 563                                 ProjFile.ReadProj(np, pbase, GeoPar.nrays);
 564
 565                                 first = sbase[np];
 566                                 last = GeoPar.nrays - 1 - first;  //should be 2 less. hstau
 567
 568
 569                                 for (nr = first; nr <= last; nr++)
 570                                 {
 571                                         raysum = pbase[nr];
 572                                         if (Blob.pix_basis)
 573                                         {
 574                                                 ray(np, nr, list, weight, &numb, &snorm);
 575                                         }
 576                                         else
 577                                         {
 578                                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 579                                         }
 580                                         if (numb != 0)
 581                                         {
 582                                                 for (nb = 0; nb < numb; nb++)
 583                                                 {
 584                                                         k = list[nb];
 585                                                         w2[k] += raysum;
 586                                                 }
 587                                         }
 588                                 }
 589
 590                         }
 591
 592                         for (i = 0; i < area; i++)
 593                         {
 594                                 w = w1[i];
 595                                 if (w == 0.0)
 596                                         w2[i] = 0.0;
 597                                 if (w > Consts.zero)
 598                                         w2[i] /= w;
 599                         }
 600                 }
 601                 else
 602                 {
 603
 604                         // LSIRT OPTIONS
 605
 606                         // WORK HOLDS PART OF B; GET GENERALIZED INVERSE IN w2
 607
 608                         for (i = 0; i < area; i++)
 609                         {
 610                                 if (w2[i] > Consts.zero)
 611                                         w2[i] = ((REAL) 1.0) / w2[i];
 612                         }
 613
 614                         for (i = 0; i < area; i++)
 615                         {
 616                                 gen_inv_dn[i] = w2[i];
 617                         }
 618
 619                         // SET UP B VECTOR FOR 'LSIRT'
 620
 621                         for (i = 0; i < area; i++)
 622                         {
 623                                 w2[i] = w1[i] * w2[i];
 624                         }
 625                 }
 626
 627                 for (i = 0; i < area; i++)
 628                 {
 629                         vector_b[i] = w2[i];
 630                 }
 631
 632                 if (bkstrt)
 633                 {
 634                         for (i = 0; i < area; i++)
 635                                 recon[i] = w2[i];
 636                 }
 637
 638                 // ITERATIONS BEGIN HERE
 639         }
 640
 641         for (i = 0; i < area; i++)
 642         {
 643                 w1[i] = vector_b[i];
 644         }
 645
 646         // WORK HOLDS B
 647
 648         for (i = 0; i < area; i++)
 649         {
 650                 w1[i] = recon[i] + w1[i] / sigma;
 651         }
 652
 653         for (i = 0; i < area; i++)
 654         {
 655                 w3[i] = w1[i];
 656         }
 657
 658         // GOING TO GENERATE P TRANSPOSE P RHO IN WORK
 659
 660         for (i = 0; i < area; i++)
 661         {
 662                 w1[i] = 0.0;
 663         }
 664
 665         for (np = 0; np < GeoPar.prjnum; np++)
 666         {
 667                 first = sbase[np];
 668                 last = GeoPar.nrays - 1 - first;    // should be 2 less. hstau
 669
 670
 671                 for (nr = first; nr <= last; nr++)
 672                 {
 673
 674                         if (Blob.pix_basis)
 675                         {
 676                                 ray(np, nr, list, weight, &numb, &snorm);
 677                         }
 678                         else
 679                         {
 680                                 Blob.bray(np, nr, list, weight, &numb, &snorm);
 681                         }
 682
 683                         if (numb != 0)
 684                         {
 685                                 psum = 0.;
 686                                 for (nb = 0; nb < numb; nb++)
 687                                 {
 688                                         k = list[nb];
 689                                         psum += recon[k];
 690                                 }
 691
 692                                 if (lsirt == 2)
 693                                         psum /= float(numb);
 694                                 if (lsirt == 3)
 695                                         psum /= float(numb * numb);
 696
 697                                 for (nb = 0; nb < numb; nb++)
 698                                 {
 699                                         w1[list[nb]] += psum;
 700                                 }
 701                         }
 702                 }
 703         }
 704
 705         for (i = 0; i < area; i++)
 706         {
 707                 recon[i] = gen_inv_dn[i];
 708         }
 709
 710
 711         // RECON HOLDS GENERALIZED INVERSE OF DN
 712         // WORK HOLDS PTRANSPOSE P RHO
 713
 714         for (i = 0; i < area; i++)
 715         {
 716                 w1[i] *= recon[i];
 717         }
 718
 719         for (i = 0; i < area; i++)
 720         {
 721                 recon[i] = w3[i];
 722         }
 723
 724         psum = 0.;
 725         for (i = 0; i < area; i++)
 726         {
 727                 ronew = recon[i] - w1[i] / sigma;
 728                 if (Modefl.lofl)
 729                         ronew = MAX0(ronew, Modefl.lower);
 730                 if (Modefl.upfl)
 731                         ronew = MIN0(ronew, Modefl.upper);
 732                 recon[i] = ronew;
 733
 734                 if (!Blob.pix_basis)
 735                 {  //if blobs
 736                         if (i % 2 == Blob.pr)
 737                         {
 738                                 psum += ronew;
 739                         }
 740                 }
 741                 else
 742                 {
 743                         psum += ronew;
 744                 }
 745         }
 746
 747         if (!normal)
 748         {
 749
 750                 return FALSE;
 751         }
 752
 753         if (fabs(psum) > Consts.zero)
 754         {
 755                 ratio = totden / psum;
 756
 757                 for (i = 0; i < area; i++)
 758                 {
 759                         recon[i] *= ratio;
 760                 }
 761
 762                 return FALSE;
 763         }
 764
 765         fprintf(output,
 766                         "\n ***current total density too small; no normalization for this iteration***"); //hstau
 767
 768         return FALSE;
 769 }