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_vcase2.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 ------------------------------------------------------------------
12 VCASE2: VARIABLE SPACING CASE 2
13 PURPOSE: FOR VARIABLE SPACING, ICASE=2 (THE TWO LINES STADLE
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::vcase2(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
70 //tan(z)=sin(z)/cos(z)
71 ///cot(z)=cos(z)/sin(z)
73 uleft = (REAL) cot(ang1)
75 uhigh = (REAL) tan(ang2);
76 kh = Fourie.nsize2 / 2 - 1;
78 for (j = 1; j <= kh; j++)
79 { // bug79 - swr 2/11/05
84 indl = (INTEGER) ceil(aleft - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
86 a = j * Fourie.delp / Fourie.deld;
88 indmin = 2 * low; // bug79 - swr 2/11/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 yy[2] = (REAL) -low - 1;
103 xx[1] = yy[1] * uleft;
104 xx[2] = yy[2] * uleft;
106 xx[4] = (REAL) low + 1;
107 yy[3] = xx[3] * uhigh;
108 yy[4] = xx[4] * uhigh;
112 for (ix = indl; ix <= indh; ix++)
113 { // bug79 - swr 2/11/05
114 xg = ix * Fourie.delp / Fourie.deld;
115 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
116 store(ix, iy, gr, gi);
121 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
124 amid = aleft + distan / (REAL) 2.0;
125 ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
128 for (ix = indl; ix <= indh; ix++)
129 { // bug79 - swr 2/11/05
132 { // bug79 - swr 2/15/05
141 store(ix, iy, gr, gi);
146 // INTERP = 3 LINEAR INTERPOLATION
149 if (t <= Consts.zero)
151 if ((1. - t) <= Consts.zero)
153 pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
154 qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
155 px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
156 qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
159 for (ix = indl; ix <= indh; ix++)
160 { // bug79 - swr 2/11/05
161 s = (ix - aleft) / distan;
162 if (s <= Consts.zero)
164 if ((1. - s) <= Consts.zero)
166 gr = pr * ((REAL) 1.0 - s) + qr * s;
167 gi = px * ((REAL) 1.0 - s) + qi * s;
168 store(ix, iy, gr, gi);
177 indh = (INTEGER) floor(ahigh + Fourie.epsln); // bug82 - changed to floor - swr - 2/20/05
182 switch (Fourie.interp)
185 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
186 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
190 for (iy = indl; iy <= indh; iy++)
191 { // bug79 - swr 2/11/05
192 yg = iy * Fourie.delp / Fourie.deld;
193 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
194 store(ix, iy, gr, gi);
198 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
200 amid = ahigh - distan / (REAL) 2.0;
201 for (iy = indl; iy <= indh; iy++)
202 { // bug79 - swr 2/11/05
205 { // bug79 - swr 2/15/05
214 store(ix, iy, gr, gi);
218 // INTERP = 3 LINEAR INTERPOLATION
220 for (iy = indl; iy <= indh; iy++)
221 { // bug79 - swr 2/11/05
222 s = (hdis + iy - (-j)) / distan;
223 if (s <= Consts.zero)
225 if ((1. - s) <= Consts.zero)
227 gr = pr * ((REAL) 1.0 - s) + qr * s;
228 gi = px * ((REAL) 1.0 - s) + qi * s;
229 store(ix, iy, gr, gi);