src/snark/foru_frinit.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/foru_frinit.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  PURPOSE: INITIALIZE ALL GLOBAL VARIABLES AND ALLOCATE SPACE
  11  IMPORTANT VARIABLES:
  12  RINC  MODIFIED PINC,FOR UNI RINC=PINC
  13  FOR VRI RINC=PINC*RADIUS/STOD
  14  DELD  SAMPLE SPACING ON TRANSFORMED PROJECTIONS
  15  DELP  SAMPLE SPACING ON FOURIER PLAN
  16  RADMAX  CUTOFF RADIUS
  17  MIDR  NSIZE1/2
  18  MIDP  NSIZE2/2
  19  EPSLN  A SMALL NUMBER TO HANDLE ROUNDING ERROR
  20  NFORW(3)  3 ELEMENTS ARRAY TO INDICATE THE DIMENSION
  21  OF THE ARRAY TO BE FORWARD TRANSFORMED
  22  NBACK(3)  3 ELEMENTS ARRAY TO INDICATE THE DIMENSION
  23  OF THE ARRAY TO BE BACKTRANSFORMED
  24  NPROJ1,NPROJ2  STARTING ADDRESSES OF ARRAYS TO HOLD
  25  TRANSFORMED PROJECTIONS
  26  NFRPLN  STARTING ADDRESS OF THE ARRAY TO REPRESENT
  27  FOURIER PLANE
  28  IFLG  INITIALLY 0, AND SET TO 1 IF ERROR OCCURS
  29  RESERV  RESERVE FOR FUTURE USE
  30  */
  31
  32 #include <cstdio>
  33
  34 #include "blkdta.h"
  35 #include "geom.h"
  36 #include "fourie.h"
  37 #include "consts.h"
  38 #include "uiod.h"
  39 #include "anglst.h"
  40 #include "infile.h"
  41
  42 #include "foru.h"
  43
  44 void foru_class::frinit()
  45 {
  46         INTEGER filter;
  47         BOOLEAN eol;
  48         INTEGER size1;
  49         INTEGER size2;
  50         INTEGER word;
  51         INTEGER k3;
  52         INTEGER k4;
  53         INTEGER k;
  54         INTEGER kh;
  55         REAL angfst;
  56         REAL sinth;
  57         REAL costh;
  58         REAL anglst;
  59         REAL angspn;
  60
  61
  62         //data fname/4hband,4hsinc,4hcosi,4hhamm/
  63         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
  64         static const INTEGER fname[4] =
  65         {
  66                         CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'),
  67                         CHAR2INT('c', 'o', 's', 'i'), CHAR2INT('h', 'a', 'm', 'm')
  68         };
  69
  70         //data samp/4hsamp/,pict/4hpict/
  71         static const INTEGER samp = CHAR2INT('s', 'a', 'm', 'p');
  72         static const INTEGER pict = CHAR2INT('p', 'i', 'c', 't');
  73
  74         // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
  75         static const INTEGER foru_codes[2] =
  76         {
  77         CHAR2INT('s', 'a', 'm', 'p'), CHAR2INT('p', 'i', 'c', 't') };
  78
  79         Fourie.iflg = 0;
  80
  81         fprintf(output, "\nf o u r i e r  r e c o n s t r u c t i o n      \n");
  82
  83         if (GeoPar.par)
  84                 Fourie.rinc = GeoPar.pinc;
  85         if (GeoPar.div)
  86                 Fourie.rinc = GeoPar.pinc * GeoPar.radius / GeoPar.stod;
  87         if (GeoPar.strip)
  88         {
  89                 fprintf(output, "\nstrip projections");
  90
  91         }
  92         if (GeoPar.div)
  93         {
  94                 fprintf(output, "\ndivergent rays, will be treated as parallel");
  95         }
  96         filter = InFile.getwrd(TRUE, &eol, fname, 4);
  97         if (eol)
  98                 error(5);
  99
 100         for (Fourie.nfiltr = 0; Fourie.nfiltr < 4; Fourie.nfiltr++)
 101         {
 102                 if (filter == fname[Fourie.nfiltr])
 103                         goto L25;
 104         }
 105         error(6);
 106
 107         L25: Fourie.cutoff = InFile.getnum(FALSE, &eol);
 108         if (eol)
 109                 error(5);
 110         // bug 151 - swr - 9/24/05
 111         if (Fourie.cutoff < 0.0)
 112                 Fourie.cutoff = (REAL) 2.0 * GeoPar.prjnum * Fourie.rinc
 113                                 / ((REAL) Consts.pi * GeoPar.nelem * GeoPar.pixsiz);
 114         if (Fourie.cutoff <= Consts.zero)
 115                 Fourie.cutoff = 1.0;
 116         if (Fourie.cutoff > 1.0)
 117                 Fourie.cutoff = 1.0;
 118         Fourie.interp = InFile.getint(FALSE, &eol);
 119         if (eol)
 120                 error(5);
 121         if ((Fourie.interp > 4) || (Fourie.interp < 1))
 122                 error(8);
 123         size1 = GeoPar.nrays + 1;
 124         size2 = GeoPar.nelem + 1;
 125         word = InFile.getwrd(FALSE, &eol, foru_codes, 2);
 126         if (!eol)
 127         {
 128
 129                 if (word == pict)
 130                         goto L4;
 131                 if (word != samp)
 132                         error(6);
 133                 size1 = InFile.getint(FALSE, &eol);
 134                 if (eol)
 135                         error(5);
 136                 if (eol)
 137                         return;
 138                 word = InFile.getwrd(FALSE, &eol, &(foru_codes[1]), 1);
 139
 140                 if (!eol)
 141                 {
 142                         if (word == pict)
 143                                 goto L4;
 144                         error(6);
 145                         return;
 146
 147                         L4: size2 = InFile.getint(FALSE, &eol);
 148                         if (eol)
 149                                 error(5);
 150                         if (eol)
 151                                 return;
 152                 }
 153         }
 154
 155         if (size1 <= GeoPar.usrays)
 156                 size1 = GeoPar.usrays + 1;
 157         if (size1 != (size1 / 2 * 2))
 158                 size1++;
 159         Fourie.nsize1 = size1;
 160         if (size2 <= GeoPar.nelem)
 161                 size2 = GeoPar.nelem + 1;
 162         if (size2 != (size2 / 2 * 2))
 163                 size2++;
 164         Fourie.nsize2 = size2;
 165         if (Fourie.iflg == 1)
 166                 return;
 167
 168         fprintf(output, "\nfilter=   %2i", Fourie.nfiltr);
 169         fprintf(output, "\ncutoff=    %5.2f", Fourie.cutoff);
 170         fprintf(output, "\ninterp=     %2i", Fourie.interp);
 171         fprintf(output, "\nnsize1=    %4i", Fourie.nsize1);
 172         fprintf(output, "\nnsize2=   %4i", Fourie.nsize2);
 173
 174         Fourie.radmax = (Fourie.nsize1 / 2) * Fourie.cutoff;
 175         if (Fourie.nfiltr == 4)
 176                 Fourie.radmax = (REAL) Fourie.nsize1 / 2;
 177         Fourie.midr = Fourie.nsize1 / 2;
 178         Fourie.midp = Fourie.nsize2 / 2;
 179         Fourie.nforw[0] = 1;
 180         Fourie.nforw[1] = Fourie.nsize1 / 2;
 181         Fourie.nforw[2] = 1;
 182         Fourie.nback[0] = Fourie.nsize2;
 183         Fourie.nback[1] = Fourie.nsize2 / 2;
 184         Fourie.nback[2] = 1;
 185
 186         if (Fourie.ntemp != NULL)
 187                 delete[] Fourie.ntemp;
 188         if (Fourie.nproj1 != NULL)
 189                 delete[] Fourie.nproj1;
 190         if (Fourie.nproj2 != NULL)
 191                 delete[] Fourie.nproj2;
 192         Fourie.ntemp = new REAL[Fourie.nsize1];
 193         Fourie.nproj1 = new REAL[Fourie.nsize1 + 4];
 194         Fourie.nproj2 = new REAL[Fourie.nsize1 + 4];
 195
 196         k3 = Fourie.nsize1 + 2;
 197         k4 = Fourie.nsize1 + 4;
 198
 199         for (k = k3; k < k4; k++)
 200         {
 201                 Fourie.nproj1[k] = 0.0;
 202                 Fourie.nproj2[k] = 0.0;
 203         }
 204         if (Fourie.nfrpln != NULL)
 205                 delete[] Fourie.nfrpln;
 206         Fourie.nfrpln = new REAL[Fourie.nsize2 * (Fourie.nsize2 + 2)];
 207         kh = Fourie.nsize2 * (Fourie.nsize2 + 2);
 208         for (k = 0; k < kh; k++)
 209         {
 210                 Fourie.nfrpln[k] = 0.0;
 211         }
 212
 213         Fourie.deld = (REAL) 1.0 / (Fourie.nsize1 * Fourie.rinc);
 214         Fourie.delp = (REAL) 1.0 / (Fourie.nsize2 * GeoPar.pixsiz);
 215         Anglst.getang(0, &angfst, &sinth, &costh);
 216         Anglst.getang(GeoPar.prjnum - 1, &anglst, &sinth, &costh);
 217         angspn = anglst - angfst;
 218         Fourie.epsln = (REAL) 1.0e-3;
 219         if (GeoPar.par && GeoPar.uni
 220                         && ((angspn < Consts.pid2) || (angspn >= Consts.pi)))
 221                 error(1);
 222         if (GeoPar.par && GeoPar.vri
 223                         && ((angspn < (0.75 * Consts.pi)) || (angspn >= Consts.pi)))
 224                 error(1);
 225         if (GeoPar.div
 226                         && ((angspn < (1.5 * Consts.pi)) || (angspn >= Consts.twopi)))
 227                 error(2);
 228
 229         return;
 230 }