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_uint.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
11 PURPOSE: GIVEN THE COORDINATES OF THE DESIRED POINT (IX,IY)
12 AND THE TWO TRANSFORMED PROJECTIONS PJ1 AND PJ2 FIND THE
13 INTERPOLATED VALUE (GR,GI) THEN CALL THE SUBROUTINE STORE
14 TO STORE THE VALUE IN APPROPRIATE PLACE IN THE ARRAY
16 THE INTERPOLATION SCHEME TO BE USED IS DETERMINED BY
17 INTERP (SEE ALSO COMMENTS IN THE MAIN PROGRAM )
30 #define disn(b,c,the) ((REAL) sqrt((b)*(b)+(c)*(c)-2.*(b)*(c) * (REAL) cos(the)))
32 void foru_class::uint1(INTEGER icase, INTEGER ix, INTEGER iy, REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
65 r = (REAL) sqrt((REAL) x * x + y * y) * Fourie.delp;
66 ang = (REAL) atan2((REAL) y, (REAL) x);
68 if (a >= Fourie.radmax)
73 s = (REAL) fabs((ang - ang1) / (ang2 - ang1));
76 switch (Fourie.interp)
79 // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
84 if ((a * cos(ang - ang1) - low) <= 0.5)
92 if ((a * cos(ang2 - ang) - low) <= 0.5)
99 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
117 // INTERP = 3 LINEAR INTERPOLATION
119 if (t <= Consts.zero)
121 if ((1.0 - t) <= Consts.zero)
123 if (s <= Consts.zero)
125 if ((1.0 - s) <= Consts.zero)
127 pr = t * pj1[indmax] + ((REAL) 1.0 - t) * pj1[indmin];
128 qr = t * pj2[indmax] + ((REAL) 1.0 - t) * pj2[indmin];
129 gr = pr * ((REAL) 1.0 - s) + qr * s;
130 px = t * pj1[indmax + 1] + ((REAL) 1.0 - t) * pj1[indmin + 1];
131 qi = t * pj2[indmax + 1] + ((REAL) 1.0 - t) * pj2[indmin + 1];
132 gi = px * ((REAL) 1.0 - s) + qi * s;
135 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
138 d1 = disn(a, ak, ang - ang1);
139 d2 = disn(a, ak + 1, ang - ang1);
140 d3 = disn(a, ak, ang2 - ang);
141 d4 = disn(a, ak + 1, ang2 - ang);
148 if (e1 <= Consts.zero)
150 if (e2 <= Consts.zero)
152 if (e3 <= Consts.zero)
154 if (e4 <= Consts.zero)
157 den = e1 + e2 + e3 + e4;
160 gr = (pj1[indmin] * e1 + pj1[indmax] * e2 + pj2[indmin] * e3
161 + pj2[indmax] * e4) / den;
163 gi = (pj1[indmin + 1] * e1 + pj1[indmax + 1] * e2 + pj2[indmin + 1] * e3
164 + pj2[indmax + 1] * e4) / den;
167 // CALL THE SUBROUTINE STORE TO STORE THE VALUE (GR,GI) AT POINT
168 // (IX,IY) IN APPROPRIATE PLACE IN THE ARRAY STARTING AT NFRPLN+1
170 store(ix, iy, gr, gi);