src/snark/art.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/art.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THIS IS A GENERAL INPLEMENTATION OF ALGEBRAIC RECONSTRUUCTION
  11  THIS SUBROUTINE IS CAPABLE OF EXECUTION THE FOLLOWING METHODS:
  12
  13  ART3 FROM
  14  A RELAXATION METHOD FOR RECONSTRUCTING OBJECTS FROM NOISY
  15  X-RAY, G.T. HERMAN, MATH. PROG.,V8(1975),1-19.
  16  ART4 FROM
  17  ITERATIVE QUADRATIC OPTIMIZATION WITH APPLICATION IN DI-
  18  AGNOSTIC RADIOLOGY, G.T. HERMAN AND A. LENT, SUNYAB
  19  DEPARTMENT OF COMPUTER SCIENCE, INTERNAL RE-
  20  PORT NO. 122.
  21  BAYESIAN FROM
  22  A COMPUTER IMPLEMENTATION OF A BAYESIAN ANALYSIS OF IMAGE
  23  RECONSTRUCTION, G.T. HERMAN , A. LENT, H. HURWITZ, AND
  24  H.P. LUNG, SUNYAB DEPARTMENT OF COMPUTER SCIENCE INTERNAL
  25  REPORT NO. 134.
  26  WITH THE FOLLOWING OPTIONAL MODIFICATIONS ON THE RECONSTRUCTION
  27  ART2-WAY CONSTRAINT
  28  REFERENCE IS THE SAME AS ART4.
  29  BINARY PATTERNS (BART)
  30  RECONSTRUCTION OF BINARY PATTERNS FROM A FEW PROJECTIONS,
  31  G.T. HERMAN, INTERNATIONAL COMPUTING SYMPOSIUM 1973,
  32  371-379.
  33
  34  IMPORTANT VARIABLES:
  35  KOUNT - NUMBER OF RAYS USED IN EACH ITERATION
  36  HALFWY- EQUALS TO (UPPER+LOWER)/2.0
  37  THTOL - EQUALS TO (UPPER-LOWER)/1000.0
  38  THRLO - (1) STARTING FROM LOWER, AND IS INCREMENTED BY
  39  (UPPER-LOWER)/(KOUNT*2**(ITER+1)) AFTER EACH RECON-
  40  STRUCTION LOOP IN THIS ROUTINE WHEN BART IS
  41  SPECIFIED (BART WILL BE EXPLAINED LATER)
  42  (2) EQUALS TO LOWER OTHERWISE
  43  THRUP - (1) STARTING FROM UPPER, AND IS DECREMENTED BY THE SAME
  44  AMOUNT AS INDICATED IN THRLO (1) AFTER EACH RECON-
  45  STRUCTION LOOP WHEN BART IS SPECIFIED
  46  (2) EQUALS TO UPPER OTHERWISE
  47
  48  THE FOLLOWING KEYWORDS ARE AVAILABLE FOR THE USER TO IMPLEMENT
  49  THIS SUBROUTINE :
  50
  51  ART3 -
  52  RECONSTRUCTION WILL BE BASED ON ART3 METHOD
  53  ART4 -
  54  RECONSTRUCTION WILL BE BASED ON ART4 METHOD
  55  BAYESIAN -
  56  RECONSTRUCTION WILL BE BASED ON BAYESIAN METHOD
  57  CONSTRAINT -
  58  EITHER "ART2" OR "BOUND" OR "BART" MUST BE SPECIFIED AFTER
  59  THIS KEYWORD.  WHEN "ART2" IS SPECIFIED, THE ART2 WAY
  60  CONSTRAINT IS USED IN THE RECONSTRUCTION.  IF "BOUND" IS
  61  SPECIFIED, THE RECONSTRUCTION IS BOUNDED BY THE BOUNDARIES
  62  "UPPER" AND "LOWER".  WHEN "BART" IS SPECIFIED, THE RECON-
  63  STRUCTED PICTURE IS ASSUMED TO BE A BINARY PATTERN(WITH
  64  GRAY LEVELS "UPPER" AND "LOWER").  AND THE PICTURE ELEMENTS AR
  65  COMPARED WITH THRLO, THRUP, AND HALFWY WITH ASSOCIATED CON-
  66  TRAINT.  WHEN LOFL IS ON, IF THE ORIGINAL PICTURE ELEMENT IS
  67  LESS THAN THRLO AND ITS RECONSTRUCTION IS GREATER THAN HALFWY,
  68  THEN THIS ELEMENT IS CHANGED TO HALFWY-THTOL.  WHEN UPFL IS
  69  ON, IF THE ORIGINAL PICTURE ELEMENT IS GREATER THAN THRUP
  70  AND ITS RECONSTRUCTION IS LESS THAN HALFWY THEN IT IS
  71  CHANGED TO HALFWY+THTOL.
  72  RELAX -
  73  THE USER CAN SPECIFY THE VALUE OF THE RELAXATION PARAMETER TO
  74  BE USED IN THE RECONSTRUCTION BY THIS OPTION.  THE RELAX-
  75  ATION PARAMETER CAN BE EITHER CONSTANT OR VARIABLE.  IF A
  76  CONSTANT IS DESIRED THE USER MUST SPECIFY "RELAX CONSTANT"
  77  FOLLOWED BY A FLOATING POINT NUMBER.  IF THE PARAMETER IS
  78  VARIABLE THEN "RELAX VARIABLE" MUST BE SPECIFIED, AND A
  79  FUNCTION RSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT)
  80  FOR CALCULATING THE VARIABLE PARAMETER MUST
  81  BE SUPPLIED BY THE USER, WHERE THE CALCULATED
  82  RELAXATION PARAMETER IS RETURNED BY ASSIGNING THE VALUE TO
  83  RSET.  THE DEFAULT IS CONSTANT 1.0.
  84  CONRELAX -
  85  THE MODIFIERS FOR CONRELAX IS THE SAME AS OPTION "RELAX".
  86  IN THE VARIABLE CASE THE NAME OF THE USER SUPPLIED
  87  FUNCTION IS CRSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT).
  88  NORM -
  89  AN INTEGER MUST BE GIVEN AFTER "NORM".  THE USER HAS HIS
  90  FREEDOM OF CHOOSING THE TYPE OF NORM (>0) USED IN THE ITER-
  91  ATIVE STEP BY SPECIFYING THE INTEGER NUMBER OF THE NORM.
  92  THE DEFAULT IS 2-NORM.
  93  STEPS -
  94  AN INTEGER MUST BE SPECIFIED AFTER THIS KEYWORD, AND THE
  95  VALUE OF THE INTEGER IS THE NUMBER OF RAYS USED FOR THE RE-
  96  CONSTRUCTION FOR EACH SNARK05 ITERATION.  THE RAYS CHOSEN
  97  ONE AFTER ANOTHER ARE DETERMINED BY THE SNARK05 "SELECT"
  98  COMMAND.  THE DEFAULT IS A NUMBER EQUIVALENT TO
  99  PRJNUM*NRAYS.
 100  NOMLZ -
 101  THE FINAL RECONSTRUCTION OF EACH ITERATION
 102  WILL BE NORMALIZED SO THAT THE AVERAGE DENSITY OF THE PICT-
 103  URE EQUALS TO AVEDEN; OTHERWISE THE RECONSTRUCTION IS
 104  NOT NORMALIZED.
 105  TOLERANCE -
 106  THERE ARE THREE POSSIBILITIES TO DEFINE THE TOLERANCE FOR
 107  ART3 AND ART4:
 108  (1) "TOLERANCE FIXED T", WHERE T IS A USER SUPPLIED FL-
 109  OATING POINT NUMBER, AND IT IS USED AS THE TOLERANCE DUR-
 110  ING THE RECONSTRUCTION.
 111  (2) "TOLERANCE VARIABLE", THE TOLERANCE IS CHANGED FOR EACH
 112  SNARK05 ITERATIVE STEP.  THE USER HAS TO SUPPLY A FUNC-
 113  TION TSET(RECON,ITER,NP,NR,I,RAYSUM,LIST,WEIGHT)
 114  SUCH THAT THE CALCULATED TOLERANCE IS ASSIGNED TO
 115  TSET AT EXIT.
 116  (3) "TOLERANCE NOISE T", THE TOLERANCE IS DEPENDENT ON THE
 117  NOISE IN THE PROJECTION DATA.  IT IS CALCULATED BY
 118  T*SQRT(VARIAN) WHERE VARIAN IS OBTAINED BY
 119  CALLING FUNCTION ERRFAC.
 120
 121  */
 122
 123 #include <cstdio>
 124 #include <cmath>
 125
 126 #include "blkdta.h"
 127 #include "blob.h"
 128 #include "geom.h"
 129 #include "fourie.h"
 130 #include "modefl.h"
 131 #include "consts.h"
 132 #include "uiod.h"
 133 #include "anglst.h"
 134 #include "errfac.h"
 135 #include "errpar.h"
 136 #include "pick.h"
 137 #include "pseudo.h"
 138
 139 #include "art.h"
 140
 141 INTEGER art_class::Init()
 142 {
 143         //Xalg1User3.Open(); //Ran, bug 266
 144         //ArtIn.delta_1 = NULL; //wei, 3/2005 //Ran, bug 266
 145         return 0;
 146 }
 147
 148 INTEGER art_class::Reset()
 149 {
 150         // Xalg1User3.Close(); //Ran, bug 266
 151         // if (ArtIn.delta_1 != NULL) delete[] ArtIn.delta_1; //wei,3/2005 //Ran, bug 266
 152         if (!ArtIn.art3f)
 153                 delete[] ArtIn.totalRays; //wei, bug 266
 154         return 0;
 155 }
 156
 157 BOOLEAN art_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
 158 {
 159
 160         static INTEGER iflagr;
 161         static INTEGER iflagt;
 162         static INTEGER iflagc;
 163         static REAL thrinc;
 164         static REAL dist;
 165         static REAL minprd;
 166         static REAL rinc2;
 167         INTEGER i;
 168         INTEGER np;
 169         INTEGER nr;
 170         REAL raysum;
 171         REAL theta;
 172         REAL sinth;
 173         REAL costh;
 174         REAL v;
 175         REAL psrsum;
 176         INTEGER numb;
 177         REAL snorm;
 178         INTEGER j;
 179         REAL diff;
 180         REAL abdiff;
 181
 182         REAL sum;
 183         REAL temp;
 184         REAL averec;
 185         REAL avedif;
 186
 187         BOOLEAN erf = true;
 188         BOOLEAN erflag;
 189         static BOOLEAN flag; //may not need to be static hstau
 190
 191         REAL pictk;
 192
 193         INTEGER area;
 194
 195         REAL pof2;
 196
 197         if (Blob.pix_basis)
 198         {
 199                 area = GeoPar.area;
 200         }
 201         else
 202         {
 203                 area = Blob.area;
 204         }
 205
 206         if (iter == 1)
 207         {
 208                 flag = FALSE;
 209                 iflagr = 1;
 210                 iflagt = 1;
 211                 iflagc = 1;
 212
 213                 readin(recon, list, weight, &erflag);
 214
 215                 if (erflag)
 216                 {
 217                         return TRUE;
 218                 }
 219
 220                 ArtIn.thrlo = Modefl.lower;
 221                 ArtIn.thrup = Modefl.upper;
 222                 thrinc = 0.0;
 223                 ArtIn.halfwy = (Modefl.upper + Modefl.lower) / (REAL) 2.0;
 224                 ArtIn.thtol = (Modefl.upper - Modefl.lower) / (REAL) 1000.0;
 225                 dist = Modefl.upper - Modefl.lower;
 226
 227                 if (ArtIn.noisef)
 228                 {
 229                         minprd = 0.0;
 230                         Fourie.rinc = GeoPar.pinc;
 231                         rinc2 = GeoPar.pinc * GeoPar.pinc;
 232                         errpar(erf);
 233                 }
 234         } // iter == 1
 235
 236         if (ArtIn.bartf)
 237         {
 238
 239                 pof2 = 1;
 240                 for (i = 0; i < iter + 1; i++)
 241                         pof2 *= 2.0;
 242
 243                 thrinc = dist / (ArtIn.kount * pof2);
 244         }
 245
 246         // ONE ITERATION
 247
 248         for (i = 0; i < ArtIn.kount; i++)
 249         {
 250                 pick(&np, &nr);
 251                 raysum = Anglst.prdta(np, nr);
 252
 253                 if (!ArtIn.rconsf)
 254                 {
 255
 256                         // GET VARIABLE RELAXATION PARAMETER
 257
 258                         ArtIn.relax = rset(recon, iter, np, nr, iflagr, raysum, list,
 259                                         weight, &flag);
 260                         if (flag)
 261                         {
 262                                 return TRUE;
 263                         }
 264                 }
 265                 if (ArtIn.noisef)
 266                 {
 267
 268                         // GET VARIABLE TOLERANCE BASED ON NOISE IN PROJECTION
 269
 270                         if (!(GeoPar.line || GeoPar.uni))
 271                         {
 272                                 Anglst.getang(np, &theta, &sinth, &costh);
 273                                 Fourie.rinc = GeoPar.pinc * MAX0(sinth, costh);
 274                                 rinc2 = Fourie.rinc * Fourie.rinc;
 275                         }
 276                         v = errfac(raysum, Fourie.rinc, rinc2);
 277                         ArtIn.toler = ArtIn.ct * (REAL) sqrt(v);
 278                 }
 279                 else
 280                 {
 281                         if (!ArtIn.tolerf)
 282                         {
 283
 284                                 // GET VARIABLE TOLERANCE
 285
 286                                 ArtIn.toler = tset(recon, iter, np, nr, iflagt, raysum, list,
 287                                                 weight, &flag);
 288                                 if (flag)
 289                                 {
 290                                         return TRUE;
 291                                 }
 292                                 if (ArtIn.toler < 0.0)
 293                                 {
 294                                         fprintf(output,
 295                                                         "\n          ***tolerance calculated at iter= %5i  np= %5i  nr= %5i  tolerance= %9.5f***",
 296                                                         iter, np, nr, ArtIn.toler); //hstau
 297
 298                                         fprintf(output,
 299                                                         "\n          ***error - tolerance can not be less than 0.0***"); //hstau
 300                                         return TRUE;
 301                                 }
 302                         }
 303                 }
 304
 305                 if (!ArtIn.crf)
 306                 {
 307
 308                         // GET VARIABLE CR
 309
 310                         ArtIn.cr = crset(recon, iter, np, nr, iflagc, raysum, list, weight,
 311                                         &flag);
 312                         if (flag)
 313                         {
 314                                 return TRUE;
 315                         }
 316                 }
 317
 318                 if (Blob.pix_basis)
 319                 {
 320                         psrsum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
 321                                         GeoPar.line, ArtIn.art2f);
 322                 }
 323                 else
 324                 {
 325                         psrsum = Blob.bpseudo(recon, np, nr, list, weight, &numb, &snorm,
 326                                         GeoPar.line, ArtIn.art2f);
 327                 }
 328
 329                 if (numb != 0)
 330                 {
 331
 332                         // CALCULATE P-NORM WHEN P IS OTHER THAN 2
 333
 334                         if (ArtIn.p != 2)
 335                         {
 336                                 snorm = 0.0;
 337                                 for (j = 0; j < numb; j++)
 338                                 {
 339                                         snorm += (REAL) pow(weight[j], ArtIn.p);
 340                                 }
 341                         }
 342                         if (snorm > Consts.zero)
 343                         {
 344                                 diff = raysum - psrsum;
 345                                 abdiff = (REAL) fabs(diff);
 346
 347                                 // BACKPROJECTION
 348                                 bkproj(recon, diff, abdiff, ArtIn.totalRays, np, nr, list,
 349                                                 weight, numb, snorm);
 350
 351                         }
 352                 }
 353                 // bug 222 - swr - 04/08/07
 354                 ArtIn.thrlo += thrinc;
 355                 ArtIn.thrup -= thrinc;
 356         }
 357
 358         if (ArtIn.nomlf)
 359         {
 360
 361                 // NORMALIZATION
 362
 363                 sum = 0.0;
 364                 for (i = 0; i < area; i++)
 365                 {
 366                         temp = recon[i];
 367                         temp = clip(temp, ArtIn.thrlo, ArtIn.thrup, 1.0);
 368                         if (!Blob.pix_basis)
 369                         { //if blobs
 370                                 if (i % 2 == Blob.pr)
 371                                 {
 372                                         sum += temp;
 373                                 }
 374                         }
 375                         else
 376                         {
 377                                 sum += temp;
 378                         }
 379                 }
 380
 381                 averec = sum / float(area);
 382                 if (!Blob.pix_basis)
 383                 {
 384                         averec *= 2.0;
 385                 }
 386                 avedif = GeoPar.aveden - averec; //aveden in both basis coincide
 387
 388                 // IN THE FOLLOWING LOOP OF 820 THE COMMENTED LINE WAS
 389                 // REPLACED WITH THE CODE INSIDE LOOP 820....10/13/87
 390                 // Dr. GABOR T HERMAN
 391
 392                 for (i = 0; i < area; i++)
 393                 {
 394
 395                         pictk = recon[i] + avedif;
 396                         if (!ArtIn.art2f)
 397                         {
 398                                 pictk = clip(pictk, ArtIn.thrlo, ArtIn.thrup, 1.0);
 399                                 if (ArtIn.bartf)
 400                                 {
 401                                         if (Modefl.lofl && (recon[i] <= ArtIn.thrlo)
 402                                                         && (pictk >= ArtIn.halfwy))
 403                                                 pictk = ArtIn.halfwy - ArtIn.thtol;
 404
 405                                         if (Modefl.upfl && (recon[i] >= ArtIn.thrup)
 406                                                         && (pictk <= ArtIn.halfwy))
 407                                                 pictk = ArtIn.halfwy + ArtIn.thtol;
 408                                 }
 409                         }
 410                         recon[i] = pictk;
 411                 }
 412         }
 413
 414         return FALSE;
 415 }
 416