src/snark/conv.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/conv.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 IMPLEMENTATION OF THE THREE BEST KNOWN
  11  CONVOLUTION METHODS, NAMELY THOSE OF BRACEWELL + RIDDLE[1],
  12  RAMACHANDRAN + LAKSHMINARAYANAN[2], AND SHEPP + LOGAN[3].
  13  THE CONTROL CARD FOLLOWING THE EXECUTE CARD CONTAINS (IN THE
  14  SNARK FREE FORMAT) THE FOLLOWING INFORMATION:
  15  1) THE DESIRED FILTER ("BAND", "SINC", "COSINE" OR "HAMMING")
  16  2) THE CUTOFF PARAMETER (OR THE VALUE OF ALPHA FOR HAMMING)
  17  3) THE INTERPOLATION METHOD DESIRED
  18  THE BAND LIMITING FILTER CORRESPONDS TO THE FILTER USED BY
  19  BRACEWELL + RIDDLE(1) AND RAMACHANDRAN + LAKSHMINARAYANAN(2).
  20  THE SINC FILTER IS GENERALIZATION OF THE FILTER USED BY SHEPP
  21  + LOGAN(2).
  22  THE COSINE FILTER IS A COMMON LOWPASS FILTER APPLIED TO THE
  23  RECONSTRUCTION PROBLEM.
  24  THE HAMMING FILTER IS THE GENERALIZED HAMMMING FILTER, A SPECI
  25  CASE OF WHICH WAS USED BY CHESLER + REIDERER(4).
  26  THE CUTOFF VALUE REPRESENTS THE FRACTION OF FOURIER SPACE
  27  POWER TO INCLUDE IN THE RECONSTRUCTION.  WHEN 'CUTOFF' IS ZERO,
  28  IT IS TAKEN TO 1.0.  WHEN 'CUTOFF' IS NEGATIVE, IT IS SET
  29  TO '2.0*PRJNUM*PINC/(PI*NELEM*PIXSIZ)'.  IF 'CUTOFF' IS GREATER
  30  THAN 1.0, IT IS SET TO 1.0.
  31  IF 'UNI' = .TRUE. (UNIFORM SPACING SPECIFIED),
  32  THEN THE CUTOFF FORMS A CIRCULAR REGION IN FOURIER SPACE
  33  WRERE THE DIAMETER OF THE CIRCLE IS 'CUTOFF/PIXSIZ'.
  34  IF 'VRI' = .TRUE. THEN THE CUTOFF FORMS A SQUARE REGION WHERE
  35  THE LENGTH OF A SIDE IS 'CUTOFF/PIXSIZ'.
  36  EXAMPLE CONTROL CARD SEQUENCE
  37  EXECUTE CONVOLUTION
  38  BAND LIMITING CUTOFF = -1.0 INTERPOLATION = 2
  39  WARNING.....
  40  THIS VERSION OF THE CONVOLUTION 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
  46  REFERENCES:
  47  [1]  BRACEWELL + RIDDLE, INVERSION OF FAN-BEAM SCANS IN RADIO
  48  ASTROMOMY, ASTROPHYS. J., V.150, PP.427-434 (1967).
  49  [2]  RAMACHANDRAN + LAKSHMINARAYANAN, THREE-DIMENSIONAL
  50  RECONSTRUCTION FROM RADIOGRAPHS AND ELECTRON MICROGRAPHS:
  51  APPLICATION OF CONVOLUTIONS INSTEAD OF FOURIER TRANSFORMS,
  52  PROC. NAT. ACAD. SCI. U.S., V.68, PP.2236-2240 (1971).
  53  [3]  SHEPP + LOGAN, THE FOURIER RECONSTRUCTION OF A HEAD SECTION,
  54  IEEE TRANS. NUCL. SCI., NS-21, PP.21-43 (1974).
  55  [4]  CHESLER + RIEDERER, RIPPLE SUPPRESSION DURING RECONSTRUCTION
  56  IN TRANSVERSE TOMOGRAPHY, PHYS. MED. BIOL., V.20, PP.632-636
  57  (1975).
  58  */
  59
  60 #include <cstdlib>
  61 #include <cstdio>
  62 #include <cmath>
  63
  64 #include "blkdta.h"
  65 #include "geom.h"
  66 #include "spctrm.h"
  67 #include "fourie.h"
  68 #include "consts.h"
  69 #include "uiod.h"
  70 #include "int2str.h"
  71 #include "anglst.h"
  72 #include "bckprj.h"
  73 #include "qinit.h"
  74 #include "projfile.h"
  75 #include "infile.h"
  76
  77 #include "conv.h"
  78
  79 BOOLEAN conv_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
  80 {
  81
  82         INTEGER filter;
  83
  84         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
  85         static const INTEGER fname[4] =
  86         {
  87                         CHAR2INT('b', 'a', 'n', 'd'),
  88                         CHAR2INT('s', 'i', 'n', 'c'),
  89                         CHAR2INT('c', 'o', 's', 'i'),
  90                         CHAR2INT('h', 'a', 'm', 'm')
  91         };
  92
  93         REAL thmod;
  94         REAL thfac;
  95         INTEGER word;
  96         BOOLEAN eol;
  97         INTEGER ndelta;
  98         INTEGER gpnum;
  99         REAL* g = NULL;        //wei 3.2005
 100         INTEGER qnum;
 101         REAL* q = NULL;        //wei 3.2005
 102         REAL* gp = NULL;        //wei 3.2005
 103         REAL theta;
 104         REAL sinth;
 105         REAL costh;
 106         REAL theta0;
 107         INTEGER np;
 108         REAL theta1;
 109         INTEGER np2;
 110         REAL theta2;
 111         REAL sinth2;
 112         REAL costh2;
 113         REAL w;
 114         INTEGER izro;
 115         REAL* gpi;
 116         INTEGER nr;
 117         REAL* gi;
 118         REAL* qi;
 119         REAL sum;
 120         INTEGER iz;
 121         INTEGER i;
 122
 123         // INITIALIZATION
 124
 125         if (iter > 1)
 126                 return FALSE;
 127         thmod = Consts.pi;
 128         thfac = 2.0;
 129         if (!GeoPar.par)
 130         {
 131                 thmod = Consts.twopi;
 132                 thfac = 4.0;
 133         }
 134
 135         word = InFile.getwrd(TRUE, &eol, fname, 4);
 136         for (filter = 0; filter < 4; filter++)
 137         {
 138                 if (word == fname[filter])
 139                         goto L30;
 140         }
 141
 142         fprintf(output, "\n***** %s is not a valid filter name", int2str(word)); //(JD 12/31/03)
 143
 144         return TRUE;
 145
 146         L30: Fourie.cutoff = InFile.getnum(FALSE, &eol);
 147         Fourie.interp = InFile.getint(FALSE, &eol);
 148
 149         if (eol)
 150         {
 151                 fprintf(output, "\n *** not enough arguments");
 152
 153                 return TRUE;
 154         }
 155
 156         if (Fourie.cutoff < 0.0)
 157                 Fourie.cutoff = (REAL) 2.0 * GeoPar.prjnum * GeoPar.pinc
 158                                 / ((REAL) Consts.pi * GeoPar.nelem * GeoPar.pixsiz);
 159         if (Fourie.cutoff <= Consts.zero)
 160                 Fourie.cutoff = 1.0; // bug 151 - swr - 9/24/05
 161         if (Fourie.cutoff > 1.0)
 162                 Fourie.cutoff = 1.0;
 163
 164         fprintf(output, "\n          convolution reconstruction algorithm");
 165         fprintf(output, "\n          using filter %s", int2str(word));
 166         fprintf(output, "\n          cutoff (or alpha) = %10.4f", Fourie.cutoff);
 167
 168         if ((Fourie.interp < -1) || (Fourie.interp > 6))
 169         {
 170                 fprintf(output, "\n *** invalid interpolation %4i", Fourie.interp);
 171
 172                 return TRUE;
 173         }
 174
 175         if (Fourie.interp > 0)
 176         {
 177                 fprintf(output, "\n         %2i point lagrange interpolation",
 178                                 Fourie.interp);
 179         }
 180         if (Fourie.interp < 0)
 181         {
 182                 fprintf(output, "\n          cubic spline interpolation");
 183         }
 184         if (Fourie.interp == 0)
 185         {
 186                 fprintf(output, "\n          band limiting (sinc) interpolation");
 187         }
 188
 189         ndelta = 2 * Fourie.interp;
 190         if (Fourie.interp < 0)
 191                 ndelta = 6;
 192         gpnum = GeoPar.snrays + ndelta;
 193
 194         // G IS THE BASE OF PROJECTION DATA
 195
 196         g = new REAL[GeoPar.usrays];
 197
 198         // Q IS THE BASE OF THE CONVOLUTING FUNCTION VALUES
 199
 200         qnum = GeoPar.nrays + ndelta;
 201
 202         q = new REAL[qnum];
 203
 204         qinit(q, qnum, filter, Fourie.cutoff);
 205
 206         // GP IS THE BASE OF THE CONVOLUTED PROJECTIONS
 207
 208         gp = new REAL[gpnum];
 209
 210         // BEGIN THE RECONSTRUCTION
 211
 212         Anglst.getang(GeoPar.prjnum - 1, &theta, &sinth, &costh);
 213
 214         theta0 = theta - thmod;
 215         for (np = 0; np < GeoPar.prjnum; np++)
 216         {
 217
 218                 ProjFile.ReadProj(np, g, GeoPar.usrays);
 219
 220                 Anglst.getang(np, &theta1, &sinth, &costh);
 221
 222                 // COMPUTE THE WEIGHT FOR THE PROJECTION
 223
 224                 np2 = np + 1;
 225
 226                 if (np == GeoPar.prjnum - 1)
 227                         np2 = 0;
 228
 229                 Anglst.getang(np2, &theta2, &sinth2, &costh2);
 230
 231                 if (np == (GeoPar.prjnum - 1))
 232                         theta2 += thmod;
 233                 w = (theta2 - theta0) / thfac;
 234                 if (w <= 0.0)
 235                 {
 236                         fprintf(output, "\n *** angles not in increasing order");
 237
 238                         if (g != NULL)
 239                                 delete[] g;
 240                         if (q != NULL)
 241                                 delete[] q;
 242                         if (gp != NULL)
 243                                 delete[] gp;         //wei 3/2005
 244
 245                         return TRUE;
 246                 }
 247
 248                 Fourie.rinc = GeoPar.pinc;
 249
 250                 if (GeoPar.vri)
 251                         Fourie.rinc = GeoPar.pinc * (REAL) MAX0(fabs(sinth), fabs(costh));
 252                 if (GeoPar.div)
 253                         Fourie.rinc = GeoPar.pinc * GeoPar.radius / GeoPar.stod;
 254
 255                 w /= Fourie.rinc;
 256
 257                 if (GeoPar.strip)
 258                         w /= Fourie.rinc;
 259
 260                 theta0 = theta1;
 261
 262                 // COMPUTE CONVOLUTED PROJECTIONS
 263
 264                 izro = (GeoPar.usrays - gpnum) / 2;
 265
 266                 gpi = gp;
 267
 268                 for (nr = 0; nr < gpnum; nr++)
 269                 {
 270                         gi = g;
 271                         qi = q + abs(izro);
 272                         sum = 0.0;
 273
 274                         iz = MIN0(izro, GeoPar.usrays);
 275                         if (iz > 0)
 276                         {
 277                                 for (i = 0; i < iz; i++)
 278                                 {
 279                                         sum += *gi * *qi;
 280                                         gi++;
 281                                         qi--;
 282                                 }
 283                         }
 284
 285                         izro++;
 286                         iz = MAX0(izro, 1);
 287                         if (iz <= GeoPar.usrays)
 288                         {
 289                                 for (i = iz; i <= GeoPar.usrays; i++)
 290                                 {
 291                                         sum += *gi * *qi;
 292                                         gi++;
 293                                         qi++;
 294                                 }
 295                         }
 296
 297                         *gpi++ = w * sum;
 298                 }
 299
 300                 // BACK PROJECT THE DATA
 301                 bckprj(recon, GeoPar.nelem, gp, gpnum, sinth, costh,
 302                                 GeoPar.pixsiz / Fourie.rinc, Fourie.interp);
 303
 304         }
 305
 306         if (g != NULL)
 307                 delete[] g;
 308         if (q != NULL)
 309                 delete[] q;
 310         if (gp != NULL)
 311                 delete[] gp;         //wei 3/2005
 312
 313         return FALSE;
 314 }
 315