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) $
8 ***********************************************************
10 PURPOSE: INITIALIZE ALL GLOBAL VARIABLES AND ALLOCATE SPACE
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
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
28 IFLG INITIALLY 0, AND SET TO 1 IF ERROR OCCURS
29 RESERV RESERVE FOR FUTURE USE
44 void foru_class::frinit()
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] =
66 CHAR2INT('b', 'a', 'n', 'd'), CHAR2INT('s', 'i', 'n', 'c'),
67 CHAR2INT('c', 'o', 's', 'i'), CHAR2INT('h', 'a', 'm', 'm')
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');
74 // modified calls to getwrd() to use the 4-parameter version. Lajos, Dec 13, 2004
75 static const INTEGER foru_codes[2] =
77 CHAR2INT('s', 'a', 'm', 'p'), CHAR2INT('p', 'i', 'c', 't') };
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");
84 Fourie.rinc = GeoPar.pinc;
86 Fourie.rinc = GeoPar.pinc * GeoPar.radius / GeoPar.stod;
89 fprintf(output, "\nstrip projections");
94 fprintf(output, "\ndivergent rays, will be treated as parallel");
96 filter = InFile.getwrd(TRUE, &eol, fname, 4);
100 for (Fourie.nfiltr = 0; Fourie.nfiltr < 4; Fourie.nfiltr++)
102 if (filter == fname[Fourie.nfiltr])
107 L25: Fourie.cutoff = InFile.getnum(FALSE, &eol);
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)
116 if (Fourie.cutoff > 1.0)
118 Fourie.interp = InFile.getint(FALSE, &eol);
121 if ((Fourie.interp > 4) || (Fourie.interp < 1))
123 size1 = GeoPar.nrays + 1;
124 size2 = GeoPar.nelem + 1;
125 word = InFile.getwrd(FALSE, &eol, foru_codes, 2);
133 size1 = InFile.getint(FALSE, &eol);
138 word = InFile.getwrd(FALSE, &eol, &(foru_codes[1]), 1);
147 L4: size2 = InFile.getint(FALSE, &eol);
155 if (size1 <= GeoPar.usrays)
156 size1 = GeoPar.usrays + 1;
157 if (size1 != (size1 / 2 * 2))
159 Fourie.nsize1 = size1;
160 if (size2 <= GeoPar.nelem)
161 size2 = GeoPar.nelem + 1;
162 if (size2 != (size2 / 2 * 2))
164 Fourie.nsize2 = size2;
165 if (Fourie.iflg == 1)
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);
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;
180 Fourie.nforw[1] = Fourie.nsize1 / 2;
182 Fourie.nback[0] = Fourie.nsize2;
183 Fourie.nback[1] = Fourie.nsize2 / 2;
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];
196 k3 = Fourie.nsize1 + 2;
197 k4 = Fourie.nsize1 + 4;
199 for (k = k3; k < k4; k++)
201 Fourie.nproj1[k] = 0.0;
202 Fourie.nproj2[k] = 0.0;
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++)
210 Fourie.nfrpln[k] = 0.0;
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)))
222 if (GeoPar.par && GeoPar.vri
223 && ((angspn < (0.75 * Consts.pi)) || (angspn >= Consts.pi)))
226 && ((angspn < (1.5 * Consts.pi)) || (angspn >= Consts.twopi)))