src/snark/rfl.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/rfl.cpp $
   5  $LastChangedRevision: 94 $
   6  $Date: 2014-07-02 18:53:30 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  RHO-FILTERED LAYERGRAM RECONSTRUCTION METHOD
  11  THE CONTROL CARD FOLLOWING THE EXECUTE CARD CONTAINS (IN THE
  12  SNARK FREE FORMAT) THE FOLLOWING INFORMATION:
  13  1) THE DESIRED FILTER ("BAND", "SINC", "COSINE" OR "HAMMING")
  14  2) THE CUTOFF PARAMETER (OR THE VALUE OF ALPHA FOR HAMMING)
  15  3) THE INTERPOLATION METHOD DESIRED
  16  4) BACK PROJECTION REGION SIZE (OPTIONAL)
  17  THE FILTERS ARE CIRCULARLY SYMMETRIC VERSIONS OF THE
  18  CORRESPONDING FILTERS USED IN THE CONVOLUTION ALGORITHM.
  19  THE BAND LIMITING FILTER CORRESPONDS TO THE FILTER USED BY
  20  BRACEWELL + RIDDLE[1] AND RAMACHANDRAN + LAKSHMINARAYANAN[2].
  21  THE SINC FILTER IS A GENERALIZATION OF THE FILTER USED BY SHEPP
  22  + LOGAN[2].
  23  THE COSINE FILTER IS A COMMON LOWPASS FILTER APPLIED TO THE
  24  RECONSTRUCTION PROBLEM.
  25  THE HAMMING FILTER IS THE GENERALIZED HAMMING FILTER, A SPECIAL
  26  CASE OF WHICH WAS USED BY CHESLER + RIEDERER[4].
  27  THE CUTOFF VALUE REPRESENTS THE FRACTION OF THE FOURIER SPACE
  28  POWER RELATIVE TO THE RAY SPACING TO INCLUDE IN THE RECONSTRUCTION
  29  A CUTOFF AT THE FOLDING FREQUENCY OF THE 2D FOURIER TRANSFORM
  30  CORRESPONDS TO CUTOFF=PINC/PIXSIZ.  WHEN "CUTOFF" IS GIVEN AS 0.0,
  31  IT IS SET TO 1.0.  WHEN "CUTOFF" IS NEGATIVE, IT IS SET TO
  32  "2.0*PRJNUM*PINC/(PI*NELEM*PIXSIZ)".
  33  THE BACK PROJECTION REGION SIZE MAY BE USER SPECIFIED AS
  34  ANY EVEN NUMBER GREATER THAN NELEM (WITHIN REASON), OR IF NOT
  35  GIVEN DEFAULTS TO NELEM+1.
  36  EXAMPLE CONTROL CARD SEQUENCE
  37  EXECUTE RFL
  38  BAND LIMITING CUTOFF = -1.0 INTERPOLATION = 2
  39  WARNING.....
  40  THIS VERSION OF THE RFL METHOD WEIGHTS PROJECTION "NP"
  41  BY "(THETA(NP+1)-THETA(NP-1))/2.0".  IN ORDER TO ACCOMPLISH THIS
  42  IT IS ASSUMED THAT THE PROJECTIONS ARE ARRANGED IN INCREASING
  43  ORDER ON THETA AND THAT THETA(PRJNUM)-THETA(1) IS LESS THAN
  44  PI IN THE PARALLEL CASE OR TWOPI IN THE DIVERGENT CASE.
  45  REFERENCES:
  46  [1]  BRACEWELL + RIDDLE, INVERSION OF FAN-BEAM SCANS IN RADIO
  47  ASTROMOMY, ASTROPHYS. J., V.150, PP.427-434 (1967).
  48  [2]  RAMACHANDRAN + LAKSHMINARAYANAN, THREE-DIMENSIONAL
  49  RECONSTRUCTION FROM RADIOGRAPHS AND ELECTRON MICROGRAPHS:
  50  APPLICATION OF CONVOLUTIONS INSTEAD OF FOURIER TRANSFORMS,
  51  PROC. NAT. ACAD. SCI. U.S., V.68, PP.2236-2240 (1971).
  52  [3]  SHEPP + LOGAN, THE FOURIER RECONSTRUCTION OF A HEAD SECTION,
  53  IEEE TRANS. NUCL. SCI., NS-21, PP.21-43 (1974).
  54  [4]  CHESLER + RIEDERER, RIPPLE SUPPRESSION DURING RECONSTRUCTION
  55  IN TRANSVERSE TOMOGRAPHY, PHYS. MED. BIOL., V.20, PP.632-636
  56  (1975).
  57  */
  58
  59 #include <cstdio>
  60 #include <cmath>
  61
  62 #include "blkdta.h"
  63 #include "geom.h"
  64 #include "spctrm.h"
  65 #include "fourie.h"
  66 #include "consts.h"
  67 #include "uiod.h"
  68
  69 #include "int2str.h"
  70 #include "rtfort.h"
  71 #include "bckprj.h"
  72 #include "qfilt.h"
  73 #include "anglst.h"
  74
  75 #include "projfile.h"
  76 #include "infile.h"
  77
  78 #include "rfl.h"
  79
  80 BOOLEAN rfl_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  81 {
  82         REAL* reco = NULL;   //wei, 3/2005
  83         REAL* g;
  84         INTEGER gnum;
  85
  86         INTEGER word;
  87
  88         REAL thmod;
  89         REAL thfac;
  90         INTEGER melen;
  91         INTEGER narea;
  92         INTEGER nn;
  93         INTEGER ndelta;
  94         REAL theta;
  95         REAL sinth;
  96         REAL costh;
  97         REAL theta0;
  98         INTEGER np;
  99         REAL theta1;
 100         INTEGER np2;
 101         REAL theta2;
 102         REAL sinth2;
 103         REAL costh2;
 104         REAL w;
 105         INTEGER nr;
 106         REAL topinc;
 107         REAL scale;
 108         INTEGER nrow;
 109         INTEGER ncol;
 110         INTEGER ind;
 111         INTEGER i;
 112         INTEGER ip;
 113         INTEGER j;
 114         INTEGER jp;
 115         REAL r;
 116         REAL factr;
 117         REAL sum;
 118         REAL diff;
 119         INTEGER index;
 120         INTEGER offset;
 121
 122         BOOLEAN eol;
 123         INTEGER n[3], filter;
 124
 125         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 16, 2004
 126         static const INTEGER fname[4] =
 127         {
 128                         CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'),
 129                         CHAR2INT('c', 'o', 's', 'i'), CHAR2INT('h', 'a', 'm', 'm') };
 130
 131         // INITIALIZATION
 132
 133         ///xalg6=.false.
 134         if (iter > 1)
 135                 return FALSE;
 136         thmod = Consts.pi;
 137         thfac = 2.0;
 138         if (!GeoPar.par)
 139         {
 140                 thmod = Consts.twopi;
 141                 thfac = 4.0;
 142         }
 143         word = InFile.getwrd(TRUE, &eol, fname, 4);
 144
 145         for (filter = 0; filter < 4; filter++)
 146         {
 147                 if (word == fname[filter])
 148                         goto L30;
 149         }
 150
 151         fprintf(output, "\n***** %4s is not a valid filter name", int2str(word));
 152         return TRUE;
 153
 154         L30: Fourie.cutoff = InFile.getnum(FALSE, &eol);
 155         Fourie.interp = InFile.getint(FALSE, &eol);
 156
 157         if (eol)
 158         {
 159                 fprintf(output, "\n *** not enough parameters");
 160                 return TRUE;
 161         }
 162
 163         melen = InFile.getint(FALSE, &eol);
 164         if (Fourie.cutoff < 0.0)
 165                 Fourie.cutoff = (REAL) 2.0 * GeoPar.prjnum * GeoPar.pinc / ((REAL) Consts.pi * GeoPar.nelem * GeoPar.pixsiz);
 166         if (Fourie.cutoff <= Consts.zero)
 167                 Fourie.cutoff = 1.0;
 168         if (Fourie.cutoff > 1.0)
 169                 Fourie.cutoff = 1.0;
 170         fprintf(output, "\n           rho-filtered layergram reconstruction algorithm");
 171         fprintf(output, "\n           using filter %4s", int2str(word));
 172         fprintf(output, "\n           cutoff (or alpha) = %9.4f", Fourie.cutoff);
 173
 174         if ((Fourie.interp < -1) || (Fourie.interp > 6))
 175         {
 176                 fprintf(output, "\n invalid interpolation %4i", Fourie.interp);
 177                 return TRUE;
 178         }
 179
 180         if (Fourie.interp > 0)
 181         {
 182                 fprintf(output, "\n         %2i point lagrange interpolation", Fourie.interp);
 183         }
 184         if (Fourie.interp < 0)
 185         {
 186                 fprintf(output, "\n           cubic spline interpolation");
 187         }
 188         if (Fourie.interp == 0)
 189         {
 190                 fprintf(output, "\n           band limiting (sinc) interpolation");
 191         }
 192         if (!(melen > 0))
 193         {
 194                 melen = GeoPar.nelem + 1;
 195         }
 196
 197         if (melen < GeoPar.nelem)
 198         {
 199                 fprintf(output, "\n***** melen is less than nelem, melen= %4i", melen);
 200                 return TRUE;
 201         }
 202
 203         fprintf(output, "\n           backprojection is performed on a %5i x %5i region", melen, melen);
 204
 205         n[0] = melen;
 206         n[1] = melen / 2;
 207         n[2] = 1;
 208
 209         if ((2 * n[1]) != melen)
 210         {
 211                 fprintf(output, "\n***** melen must be even, melen= %4i", melen);
 212                 return TRUE;
 213         }
 214
 215         narea = (melen + 2) * melen;
 216
 217         reco = new REAL[narea];
 218
 219         for (nn = 0; nn < narea; nn++)
 220         {
 221                 reco[nn] = 0.0;
 222         }
 223
 224         // G IS THE BASE OF PROJECTION DATA
 225
 226         ndelta = 2 * Fourie.interp;
 227         if (Fourie.interp < 0)
 228                 ndelta = 6;
 229         gnum = GeoPar.numray(melen, GeoPar.pixsiz) + ndelta;
 230
 231         g = new REAL[gnum];
 232
 233         // BEGIN THE RECONSTRUCTION
 234
 235         Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
 236
 237         theta0 = theta - thmod;
 238
 239         for (np = 0; np < GeoPar.prjnum; np++)
 240         {
 241
 242                 ProjFile.ReadProj(np, g, gnum);
 243
 244                 Anglst.getang(np, &theta1, &sinth, &costh);
 245
 246                 // COMPUTE THE WEIGHT FOR THE PROJECTION
 247
 248                 np2 = np + 1;
 249
 250                 if (np == (GeoPar.prjnum - 1))
 251                         np2 = 0;
 252
 253                 Anglst.getang(np2, &theta2, &sinth2, &costh2);
 254
 255                 if (np == (GeoPar.prjnum - 1))
 256                         theta2 += thmod;
 257
 258                 w = (theta2 - theta0) / thfac;
 259
 260                 if (w <= 0.0)
 261                 {
 262                         fprintf(output, "\n *** angles are not in increasing order");
 263                         if (reco != NULL)
 264                                 delete[] reco;   //wei 3/2005
 265                         return TRUE;
 266                 }
 267
 268                 Fourie.rinc = GeoPar.pinc;
 269                 if (GeoPar.div)
 270                         Fourie.rinc *= GeoPar.radius / GeoPar.stod;
 271                 if (GeoPar.vri)
 272                         Fourie.rinc = (REAL) (Fourie.rinc * MAX0(fabs(sinth), fabs(costh)));
 273                 if (GeoPar.strip)
 274                         w /= Fourie.rinc;
 275                 theta0 = theta1;
 276
 277                 for (nr = 0; nr < gnum; nr++)
 278                 {
 279                         g[nr] *= w;
 280                 }
 281                 // BACK PROJECT THE DATA
 282                 bckprj(reco, melen, g, gnum, sinth, costh, GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
 283
 284         }
 285
 286         // COMPUTE 2-D FOURIER TRANSFORM
 287         delete[] g;  // bug 92 - Lajos - 03/02/2005
 288
 289         rtfort(reco, n, 1);
 290
 291         // MULTIPLY BY FILTER
 292         topinc = (REAL) (2.0 * GeoPar.pinc);
 293
 294         scale = topinc / (melen * GeoPar.pixsiz);
 295
 296         nrow = melen / 2 + 1;
 297         ncol = melen;
 298         ind = -2;
 299
 300         for (i = 0; i < ncol; i++)
 301         {
 302                 ip = MIN0(i, melen - i);
 303                 for (j = 0; j < nrow; j++)
 304                 {
 305                         jp = j;
 306                         ind += 2;
 307                         r = (REAL) (sqrt(((REAL) (ip * ip + jp * jp))) * scale);
 308                         factr = qfilt(r, Fourie.cutoff, filter) / topinc;
 309
 310                         reco[ind] *= factr;
 311                         reco[ind + 1] *= factr;
 312                 }
 313         }
 314
 315         // INVERSE FOURIER TRANSFORM
 316         rtfort(reco, n, -1);
 317
 318         // SCALE TO ESTIMATED AVERAGE DENSITY
 319         sum = 0.0;
 320         offset = (melen - 1 - GeoPar.nelem) / 2;
 321
 322         ind = 0;
 323         for (i = 0; i < GeoPar.nelem; i++)
 324         {
 325                 index = (i + offset) * melen + offset;
 326                 for (j = 0; j < GeoPar.nelem; j++)
 327                 {
 328                         recon[ind] = reco[index];
 329                         sum += reco[index];
 330                         ind++;
 331                         index++;
 332                 }
 333         }
 334
 335         diff = GeoPar.aveden - sum / GeoPar.area;
 336
 337         for (i = 0; i < GeoPar.area; i++)
 338         {
 339                 recon[i] += diff;
 340         }
 341
 342         if (reco != NULL)
 343                 delete[] reco;   //wei,3/2005
 344                 return FALSE;
 345 }
 346