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_vcase4.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ------------------------------------------------------------------
12 VCASE4: VARIABLE SPACING CASE 4
13 PURPOSE: FOR VARIABLE SPACING, ICASE=4 ( THE TWO LINES STRADLE
14 THE 45 DEGREE LINE) 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 #define cot(z) (cos(z)/sin(z));
33 void foru_class::vcase4(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
70 ///tan(z)=sin(z)/cos(z)
71 ///cot(z)=cos(z)/sin(z)
73 ulow = (REAL) tan(ang1);
74 uleft = (REAL) cot(ang2)
76 kh = Fourie.nsize2 / 2 - 1;
78 for (j = 1; j <= kh; j++)
79 { // bug79 - swr 2/11/05
82 a = j * Fourie.delp / Fourie.deld;
84 indmin = 2 * low; // bug79 - swr 2/11/05
91 indl = (INTEGER) ceil(alow - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
94 switch (Fourie.interp)
97 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
98 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
102 xx[2] = (REAL) low + 1;
103 yy[1] = xx[1] * ulow;
104 yy[2] = xx[2] * ulow;
106 yy[4] = (REAL) low + 1;
107 xx[3] = yy[3] * uleft;
108 xx[4] = yy[4] * uleft;
113 for (iy = indl; iy <= indh; iy++)
114 { // bug79 - swr 2/11/05
115 yg = iy * Fourie.delp / Fourie.deld;
116 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
117 store(ix, iy, gr, gi);
122 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
125 amid = alow + distan / (REAL) 2.0;
126 ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
130 for (iy = indl; iy <= indh; iy++)
131 { // bug79 - swr 2/11/05
134 { // bug79 - swr 2/15/05
143 store(ix, iy, gr, gi);
148 // INTERP = 3 LINEAR INTERPOLATION
151 if (t <= Consts.zero)
153 if ((1. - t) <= Consts.zero)
155 pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
156 qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
157 px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
158 qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
163 for (iy = indl; iy <= indh; iy++)
164 { // bug79 - swr 2/11/05
165 s = (iy - alow) / distan;
166 if (s <= Consts.zero)
168 if (((REAL) 1.0 - s) <= Consts.zero)
170 gr = pr * ((REAL) 1.0 - s) + qr * s;
171 gi = px * ((REAL) 1.0 - s) + qi * s;
172 store(ix, iy, gr, gi);
179 indl = (INTEGER) ceil(aleft - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
185 switch (Fourie.interp)
188 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
189 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
193 for (ix = indl; ix <= indh; ix++)
194 { // bug79 - swr 2/11/05
195 xg = ix * Fourie.delp / Fourie.deld;
196 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
197 store(ix, iy, gr, gi);
201 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
203 amid = aleft + distan / (REAL) 2.0;
204 for (ix = indl; ix <= indh; ix++)
205 { // bug79 - swr 2/11/05
208 { // bug79 - swr 2/15/05
217 store(ix, iy, gr, gi);
221 // INTERP = 3 LINEAR INTERPOLATION
223 for (ix = indl; ix <= indh; ix++)
224 { // bug79 - swr 2/11/05
225 s = (vdis + j - ix) / distan;
226 if (s <= Consts.zero)
228 if ((1. - s) <= Consts.zero)
230 gr = pr * ((REAL) 1.0 - s) + qr * s;
231 gi = px * ((REAL) 1.0 - s) + qi * s;
232 store(ix, iy, gr, gi);