src/snark/lino.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/lino.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ******************************************************************
  11
  12  - LINOGRAM RECONSTRUCTION ALGORITHM -
  13
  14  THIS PROGRAM IMPLEMENTS AN IMAGE RECONSTRUCTION TECHNIQUE
  15  DERIVED BY PAUL EDHOLM AND GABOR T. HERMAN.  DATA IS COLLECTED
  16  IN A SPECIAL GEOMETRY,AND THE IMAGE IS THEN RECONSTRUCTED
  17  (ESSENTIALLY) BY USING 3 ONE-DIMENSIONAL FOURIER TRANSFORMS.
  18
  19  VERSION 1.0 :     LYLE BERKOWITZ  AUGUST, 1986
  20  1.1 :     DAVID A. ROBERTS   MAY, 1987
  21  :     RESTRUCTURED PROGRAM
  22  ADDED GRIDDING ALGORITHM
  23  1.2 :     DAVID A. ROBERTS  JUNE, 1987
  24  :     USE OF RTFORT TO IMPROVE SPEED
  25  1.3 :     DAVID A. ROBERTS  JULY, 1987
  26  :     USE OF CHIRP-Z ALGORITHM
  27  TO COMPUTE SECOND DFT
  28  1.4 :     D. ODHNER  JAN, 1989
  29  ******************************************************************
  30  */
  31
  32 #include <cstdio>
  33 #include <cmath>
  34
  35 #include "blkdta.h"
  36 #include "geom.h"
  37 #include "consts.h"
  38 #include "uiod.h"
  39 #include "int2str.h"
  40 #include "anglst.h"
  41 #include "rtfort.h"
  42 #include "projfile.h"
  43 #include "infile.h"
  44 #include "qfilt.h"
  45
  46 #include "lino.h"
  47
  48 BOOLEAN lino_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  49 {
  50
  51         static const REAL padwgt = 3.0;
  52
  53         //     NEW VARIABLES
  54
  55         INTEGER i, j, k, r1, ii, iii, iflag, i3, kk;
  56         INTEGER b, word, filter, n3[3], np, np1, np2;
  57         INTEGER nnew;
  58         REAL* img;
  59         INTEGER glnum, gl1num, imgnum;
  60         REAL w, r, theta, sinth, costh, dold, dnew, dist, rinc;
  61         REAL sum, diff, c, cutoff, factr1, factr2, factr3;
  62         BOOLEAN eol;
  63
  64         REAL* g;
  65         REAL* gl;
  66         REAL *gl1;
  67         REAL* otemp;
  68         REAL* temp;
  69
  70         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
  71         static const INTEGER fname[3] =
  72         { CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'), CHAR2INT('c',
  73                         'o', 's', 'i') };
  74
  75         INTEGER intflag = 2;
  76
  77         //commented out since it is never used. Lajos, Dec 16, 2004
  78         //static const INTEGER pad = CHAR2INT('p','a','d',' ');
  79
  80         int n;
  81
  82 //-------------------------------------------------------------------
  83 //...              START MAIN PROGRAM
  84 //------------------------------------------------------------------
  85
  86         // INITIALIZATION
  87
  88         factr2 = 1.0;
  89
  90         // READ IN FILTER NAME FROM USER INPUT FILE
  91
  92         word = InFile.getwrd(TRUE, &eol, fname, 3);
  93         for (filter = 0; filter < 3; filter++)
  94         {
  95                 if (word == fname[filter])
  96                         goto L0003;
  97         }
  98
  99         fprintf(output, "\n **** %4s is not a valid filter name", int2str(word));
 100         fprintf(output, "\n **** LINOGRAM execution aborted\n");
 101         return TRUE;
 102
 103         L0003: cutoff = InFile.getnum(FALSE, &eol);
 104
 105         if (eol)
 106         {
 107                 fprintf(output, "\n **** not enough parameters");
 108                 fprintf(output, "\n **** LINOGRAM execution aborted\n");
 109                 return TRUE;
 110         }
 111
 112         if (cutoff <= 0)
 113                 cutoff = (REAL) 0.0001;
 114
 115         fprintf(output, "\n         linogram reconstruction algorithm");
 116         fprintf(output, "\n          using filter %4s", int2str(word));
 117         fprintf(output, "\n          cutoff = %9.4f", cutoff);
 118
 119         // INITIALIZE LOOP INVARIANTS
 120
 121         nnew = (GeoPar.nelem - 1) / 2;
 122         n3[0] = 1;
 123         n3[2] = 1;
 124
 125         // Decide whether to pad
 126         k = 0;
 127
 128         while ((1 << k) < (GeoPar.nelem + 1))
 129         {
 130                 k++;
 131         }
 132
 133         if (((1 << k) + padwgt * 2 * k)
 134                         <= (GeoPar.nelem + 1 + padwgt * sumfac(GeoPar.nelem + 1)))
 135         {
 136                 n3[1] = 2;
 137
 138                 while (n3[1] < (GeoPar.nelem + 1))
 139                 {
 140                         n3[1] *= 2;
 141                 }
 142
 143         }
 144         else
 145         {
 146                 n3[1] = GeoPar.nelem + 1;
 147         }
 148
 149         dold = (REAL) 2.0 / ((REAL) (GeoPar.usrays));
 150         dist = 1 / (2 * n3[1] * GeoPar.pixsiz);
 151         gl1num = 2 * n3[1] + 2;
 152         glnum = gl1num * GeoPar.usrays;
 153         imgnum = GeoPar.nelem * gl1num;
 154
 155         // G IS THE BASE OF PROJECTION DATA
 156
 157         g = new REAL[GeoPar.usrays + 2];
 158
 159         //******************************************************C
 160         //...  BEGIN THE RECONSTRUCTION
 161         //******************************************************C
 162
 163         // FIRST LOOP ON SETS OF PROJECTION ANGLES
 164         //   IFLAG = 0 FOR THE  FIRST SET OF ANGLES (T)
 165         //   IFLAG = 1 FOR THE SECOND SET OF ANGLES (C)
 166
 167         for (iflag = 0; iflag <= 1; iflag++)
 168         {
 169                 // GET THE FIRST AND LAST PROJECTION NUMBERS
 170
 171                 if (iflag == 0)
 172                 {
 173                         np1 = 0;
 174                         np2 = GeoPar.usrays;
 175                 }
 176                 else
 177                 {
 178                         np1 = GeoPar.usrays;
 179                         np2 = 2 * GeoPar.usrays;
 180                 }
 181
 182                 img = new REAL[imgnum];
 183                 for (ii = 0; ii < imgnum; ii++)
 184                 {
 185                         img[ii] = 0.0;
 186                 }
 187
 188                 gl = new REAL[glnum];
 189                 for (ii = 0; ii < glnum; ii++)
 190                 {
 191                         gl[ii] = 0.0;
 192                 }
 193
 194                 gl1 = new REAL[gl1num];
 195                 for (ii = 0; ii < gl1num; ii++)
 196                 {
 197                         gl1[ii] = 0.0;
 198                 }
 199
 200                 b = 2;
 201                 for (np = np1; np < np2; np++)
 202                 {
 203                         for (ii = 0; ii < GeoPar.usrays + 2; ii++)
 204                         {
 205                                 g[ii] = 0.0;
 206                         }
 207
 208                         ProjFile.ReadProj(np, g, GeoPar.usrays);
 209
 210                         Anglst.getang(np, &theta, &sinth, &costh);
 211
 212                         //-------------------------------------------------------------------C
 213                         //...  1) READ IN THE PROJECTION DATA BACKWARDS
 214                         //...  2) MAKE SURE GL1 ARRAY IS IN PROPER FORMAT BEFORE FFT
 215                         //...     (I.E. WITH THE ZERO VALUE AT THE BEGINNING)
 216                         //-------------------------------------------------------------------C
 217
 218                         i3 = 0;
 219                         for (i = GeoPar.nelem; i >= 0; i--)
 220                         {
 221                                 gl1[i3] = g[i];
 222                                 i3++;
 223                         }
 224
 225                         for (i = GeoPar.usrays + 1; i <= 2 * n3[1]; i++)
 226                         {
 227                                 gl1[i3] = 0.0;
 228                                 i3++;
 229                         }
 230
 231                         for (i = GeoPar.usrays; i >= GeoPar.nelem + 2; i--)
 232                         {
 233                                 gl1[i3] = g[i - 1];
 234                                 i3++;
 235                         };
 236
 237                         //------------------------------------------------------------C
 238                         //...  TAKE FIRST 1-D FOURIER TRANSFORM; USRAY-DIRECTION
 239                         //...  ** N3(2) IS ASSIGNED HALF THE NUMBER OF ELEMENTS **
 240                         //------------------------------------------------------------C
 241
 242                         rtfort(gl1, n3, 1);
 243
 244                         //------------------------------------------------------------C
 245                         //...  MULTIPLY BY FILTER AND WEIGHT (AND SAMPLE SPACING)
 246                         //------------------------------------------------------------C
 247
 248                         w = (REAL) fabs(sinth);
 249                         if (iflag == 1)
 250                                 w = (REAL) fabs(costh);
 251                         c = cutoff / (2 * GeoPar.pinc * w);
 252                         j = 0;
 253                         rinc = dist / w;
 254                         r = 0.0;
 255                         for (r1 = 0; r1 <= n3[1]; r1++)
 256                         {
 257                                 if (r <= Consts.zero)
 258                                 {
 259                                         gl1[j] = 0.0;
 260                                         gl1[j + 1] = 0.0;
 261                                 }
 262                                 else
 263                                 {
 264                                         factr1 = GeoPar.pinc * w * w;
 265                                         if (intflag == 1)
 266                                                 factr2 = qfilt(r, c, 1) / r;
 267                                         if (intflag == 2)
 268                                         {
 269                                                 factr2 = qfilt(r, c, 1) / r;
 270                                                 factr2 = factr2 * factr2;
 271                                         }
 272                                         factr3 = qfilt(r, c, filter);
 273                                         gl1[j] = factr1 * factr2 * factr3 * gl1[j];
 274                                         gl1[j + 1] = factr1 * factr2 * factr3 * gl1[j + 1];
 275                                 }
 276
 277                                 //-------------------------------------------------------------------C
 278                                 //...  STORE DATA IN NEW ARRAY (2-DIMENSIONAL): [NELEM+2,4NELEM+2]
 279                                 //-------------------------------------------------------------------C
 280
 281                                 gl[(r1 + 1) * 2 * GeoPar.usrays - b] = gl1[j];
 282                                 gl[(r1 + 1) * 2 * GeoPar.usrays - b + 1] = gl1[j + 1];
 283                                 j += 2;
 284                                 r += rinc;
 285                         }
 286                         b += 2;
 287                 }
 288
 289                 //************************************C
 290                 //  END OF LOOP ON PROJECTIONS
 291                 //************************************C
 292
 293                 //------------------------------------------------------------------
 294                 //...  PERFORM GRIDDING (INCORPORATES A 1-D FFT IN ANGULAR VARIABLE)
 295                 //...  STORE OUTPUT OF GRIDDING IN TWO DIMENSIONAL ARRAY:
 296                 //...  [NELEM+1,2*NELEM]
 297                 //-------------------------------------------------------------------
 298
 299                 otemp = new REAL[2 * GeoPar.nelem];
 300
 301                 iii = 0;
 302                 for (n = 0; n < n3[1] + 1; n++)
 303                 {
 304                         if (n == 0)
 305                         {
 306
 307                                 // SPECIAL CASE WHEN n = 1
 308                                 // DUE TO NATURE OF FILTERING, DC VALUE WILL ALWAYS BE ZERO.
 309
 310                                 for (kk = 0; kk < 2 * GeoPar.nelem; kk++)
 311                                 {
 312                                         otemp[kk] = 0.0;
 313                                 }
 314                         }
 315                         else
 316                         {
 317
 318                                 // ****  CHANGED DNEW FOR RTFORT!     ****
 319                                 //   (DUE TO PADDING WITH AN EXTRA ZERO)
 320
 321                                 dnew = ((REAL) (n)) / (2 * n3[1]);
 322
 323                                 // ****  USE CHIRP-Z ****
 324
 325                                 temp = &(gl[(n) * 2 * GeoPar.usrays]);
 326                                 czt(temp, otemp, dold, dnew, GeoPar.nelem, nnew);
 327                         }
 328
 329                         //-------------------------------------------------------------------
 330                         //...  WRITE 'OTEMP' INTO THE 2D ARRAY 'IMG' ON A COLUMN-BY-COLUMN BASIS
 331                         //-------------------------------------------------------------------
 332
 333                         ii = 0;
 334                         for (k = 0; k < GeoPar.nelem; k++)
 335                         {
 336                                 img[k * gl1num + iii] = otemp[ii];
 337                                 img[k * gl1num + iii + 1] = otemp[ii + 1];
 338                                 ii += 2;
 339                         }
 340                         iii += 2;
 341                 }
 342                 //-----------------------------------------------------------------*
 343                 //...  PERFORM AN INVERSE FFT ALONG ROWS IN THE ARRAY 'IMG'.
 344                 //...  THE FOURIER SAMPLES ARE ARRANGED IN EACH ROW STARTING FROM
 345                 //...  0 (DC).
 346                 //-----------------------------------------------------------------*
 347
 348                 // LOOP ON ROWS OF 'IMG' (SHOULD BE 'NELEM' ROWS)
 349
 350                 for (n = 0; n < GeoPar.nelem; n++)
 351                 {
 352                         rtfort(&(img[n * gl1num]), n3, -1);
 353
 354                         //----------------------------------------------C
 355                         //     MULTIPLY BY SAMPLE SPACING AFTER FFT:
 356                         //----------------------------------------------C
 357
 358                         for (k = 0; k < 2 * n3[1]; k++)
 359                         {
 360                                 img[(n) * gl1num + k] /= (GeoPar.pixsiz);
 361                         }
 362                 }
 363
 364                 //-----------------------------------------------------------------*
 365                 //...  WE NOW HAVE THE VALUES ALONG ROWS/COLUMNS IN THE IMAGE
 366                 //...  IN THE ROWS OF 'IMG'.  ALL WE HAVE TO DO NOW IS OUTPUT
 367                 //...  THESE VALUES IN THE PROPER ORDER INTO THE ARRAY 'RECON'.
 368                 //-----------------------------------------------------------------*
 369
 370                 // LOOP ON ROWS OF 'IMG' (SHOULD BE 'NELEM' ROWS)
 371
 372                 iii = 0;
 373                 for (n = GeoPar.nelem - 1; n >= 0; n--)
 374                 {
 375
 376                         // LOOP ON LAST 'N' COLUMNS IN EACH ROW [N = (NELEM-1)/2]
 377
 378                         ii = 0;
 379                         for (k = 2 * n3[1] - nnew; k < 2 * n3[1]; k++)
 380                         {
 381                                 if (iflag == 0)
 382                                 {
 383                                         recon[(n) * GeoPar.nelem + ii] += img[(iii) * gl1num + k];
 384                                 }
 385                                 else
 386                                 {
 387                                         recon[(GeoPar.nelem - ii - 1) * GeoPar.nelem + n] +=
 388                                                         img[(iii) * gl1num + k];
 389                                 }
 390                                 ii++;
 391                         }
 392
 393                         // LOOP ON FIRST 'N+1' COLUMNS IN EACH ROW
 394                         // (USING THE SAME INDEX, 'ii', TO INDEX RECON ARRAY)
 395
 396                         for (k = 0; k < nnew + 1; k++)
 397                         {
 398                                 if (iflag == 0)
 399                                 {
 400                                         recon[(n) * GeoPar.nelem + ii] += img[(iii) * gl1num + k];
 401                                 }
 402                                 else
 403                                 {
 404                                         recon[(GeoPar.nelem - ii - 1) * GeoPar.nelem + n] +=
 405                                                         img[(iii) * gl1num + k];
 406                                 }
 407                                 ii++;
 408                         }
 409                         iii++;
 410                 }
 411                 delete[] otemp;
 412                 delete[] img;
 413                 delete[] gl;
 414                 delete[] gl1;
 415         }
 416         delete[] g;
 417
 418         //*************************************************C
 419         //  END OF LOOP ON SETS OF PROJECTION ANGLES
 420         //*************************************************C
 421
 422         sum = 0.0;
 423
 424
 425         for (n = 0; n < GeoPar.area; n++)
 426         {
 427                 sum += recon[n];
 428         }
 429         sum /= (REAL) (GeoPar.nelem * GeoPar.nelem);
 430
 431         diff = (REAL) fabs(GeoPar.aveden - sum);
 432         for (n = 0; n < GeoPar.area; n++)
 433         {
 434                 recon[n] += diff;
 435         }
 436         return FALSE;
 437 }