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_vc56.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 VC56: VARIABLE SPACING CASE 5 AND 6
11 PURPOSE: FOR VARIABLE SPACING, ICASE=5 OR 6 (THE TWO LINES LIE
12 BETWEEN 45 AND 90 DEGREE LINES OR ANG1 IS IN 1ST QUADRANT
13 AND ANG2 IS IN 4TH QUADRANT) FIND ALL POINTS THAT LIE
14 BETWEEN THE TWO LINES AND INTERPOLATE
15 PJ1 AND PJ2 ARE THE ARRAYS HOLDING THE TRANSFORMED PROJS
16 ANG1 AND ANG2 ARE THE CORRESPONDING ANGLES
17 FOR CASE 6 THE SECOND ANGLE, ANG2, WILL BE ADDED PI TO
18 BECOME AN ANGLE IN SECOND QUADRANT TO PERFORM
22 // bug79 - problem with fourier algorithm
32 #define cot(z) (cos(z)/sin(z));
34 void foru_class::vc56(INTEGER icase, REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
70 ///cot(z)=cos(z)/sin(z)
75 khd = Fourie.nsize1 + 2;
77 for (kk = 1; kk < khd; kk += 2)
78 { // bug79 - swr 2/11/05
85 uleft = (REAL) cot(ang2)
87 uright = (REAL) cot(ang1)
89 kh = Fourie.nsize2 / 2 - 1;
91 for (j = 1; j <= kh; j++)
92 { // bug79 - swr 2/11/05
97 indl = (INTEGER) ceil(aleft - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
98 indh = (INTEGER) floor(aright + Fourie.epsln); // bug82 - changed to floor - swr - 2/20/05
102 a = j * Fourie.delp / Fourie.deld;
104 indmin = 2 * low; // bug79 - swr 2/11/05
106 distan = aright - aleft;
108 switch (Fourie.interp)
110 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
111 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
116 yy[2] = (REAL) low + 1;
117 yy[4] = (REAL) low + 1;
118 xx[1] = yy[1] * uright;
119 xx[2] = yy[2] * uright;
120 xx[3] = yy[3] * uleft;
121 xx[4] = yy[4] * uleft;
124 for (ix = indl; ix <= indh; ix++)
125 { // bug79 - swr 2/11/05
126 xg = ix * Fourie.delp / Fourie.deld;
127 vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
128 store(ix, iy, gr, gi);
132 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
134 amid = (aleft + aright) / (REAL) 2.0;
135 ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
136 for (ix = indl; ix <= indh; ix++)
137 { // bug79 - swr 2/11/05
140 { // bug79 - swr 2/15/05
149 store(ix, iy, gr, gi);
153 // INTERP = 3 LINEAR INTERPOLATION
156 if (t <= Consts.zero)
158 if ((1. - t) <= Consts.zero)
160 pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
161 qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
162 px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
163 qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
165 for (ix = indl; ix <= indh; ix++)
166 { // bug79 - swr 2/11/05
167 s = (aright - ix) / distan;
168 if (s <= Consts.zero)
170 if ((1. - s) <= Consts.zero)
172 gr = pr * ((REAL) 1.0 - s) + qr * s;
173 gi = px * ((REAL) 1.0 - s) + qi * s;
174 store(ix, iy, gr, gi);
184 khd = Fourie.nsize1 + 2;
185 for (kk = 1; kk < khd; kk += 2)
186 { // bug79 - swr 2/11/05