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_vcase3.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ------------------------------------------------------------------
12 VCASE3: VARIABLE SPACING CASE 3
13 PURPOSE: FOR VARIABLE SPACING, ICASE=3 (THE TWO LINES LIE BETWEEN
14 -45 AND 45 DEGREE LINES) FIND ALL POINTS THAT LIE
15 BETWEEN THE TWO LINES AND INTERPOLATE
16 PJ1 AND PJ2 ARE THE ARRAYS HOLDING THE TRANSFORMED PROJS
17 ANG1 AND ANG2 ARE THE CORRESPONDING ANGLES
19 ------------------------------------------------------------------
21 // bug79 - problem with fourier algorithm
31 void foru_class::vcase3(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
66 //tan(z)=sin(z)/cos(z)
68 ulow = (REAL) tan(ang1);
69 uhigh = (REAL) tan(ang2);
70 kh = Fourie.nsize2 / 2 - 1;
72 for (j = 1; j <= kh; j++)
73 { // bug79 - swr 2/11/05
78 indl = (INTEGER) ceil(alow - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
79 indh = (INTEGER) floor(ahigh + Fourie.epsln); // bug82 - changed to floor - swr - 2/20/05
82 a = j * Fourie.delp / Fourie.deld;
84 indmin = 2 * low; // bug79 - swr 2/11/05
86 distan = ahigh - alow;
88 switch (Fourie.interp)
91 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
92 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
97 xx[2] = (REAL) low + 1;
99 xx[4] = (REAL) low + 1;
100 yy[1] = xx[1] * ulow;
101 yy[2] = xx[2] * ulow;
102 yy[3] = xx[3] * uhigh;
103 yy[4] = xx[4] * uhigh;
105 for (iy = indl; iy <= indh; iy++)
106 { // bug79 - swr 2/11/05
107 yg = iy * Fourie.delp / Fourie.deld;
108 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
109 store(ix, iy, gr, gi);
114 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
117 amid = (alow + ahigh) / (REAL) 2.0;
118 for (iy = indl; iy <= indh; iy++)
119 { // bug79 - swr 2/11/05
121 ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
124 { // bug79 - swr 2/15/05
133 store(ix, iy, gr, gi);
137 // INTERP = 3 LINEAR INTERPOLATION
141 if (t <= Consts.zero)
143 if (((REAL) 1.0 - t) <= Consts.zero)
145 pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
146 qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
147 px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
148 qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
150 for (iy = indl; iy <= indh; iy++)
151 { // bug79 - swr 2/11/05
152 s = (iy - alow) / distan;
153 if (s <= Consts.zero)
155 if (((REAL) 1.0 - s) <= Consts.zero)
157 gr = pr * ((REAL) 1.0 - s) + qr * s;
158 gi = px * ((REAL) 1.0 - s) + qi * s;
159 store(ix, iy, gr, gi);