src/snark/foru_vc56.cpp

   1 /*
   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) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  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
  19  INTERPOLATION
  20  */
  21
  22 // bug79 - problem with fourier algorithm
  23 #include <cstdio>
  24 #include <cmath>
  25
  26 #include "blkdta.h"
  27 #include "fourie.h"
  28 #include "consts.h"
  29
  30 #include "foru.h"
  31
  32 #define  cot(z)  (cos(z)/sin(z));
  33
  34 void foru_class::vc56(INTEGER icase, REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
  35 {
  36         REAL xx[4];
  37         REAL yy[4];
  38
  39         INTEGER khd;
  40         INTEGER kk;
  41         REAL uleft;
  42         REAL uright;
  43         INTEGER kh;
  44         INTEGER j;
  45         INTEGER iy;
  46         REAL aleft;
  47         REAL aright;
  48         INTEGER indl;
  49         INTEGER indh;
  50         REAL a;
  51         INTEGER low;
  52         INTEGER indmin;
  53         INTEGER indmax;
  54         REAL distan;
  55         REAL yg;
  56         INTEGER ix;
  57         REAL xg;
  58         REAL gr;
  59         REAL gi;
  60         REAL amid;
  61         INTEGER ind;
  62         INTEGER x;
  63         REAL t;
  64         REAL pr;
  65         REAL qr;
  66         REAL px;
  67         REAL qi;
  68         REAL s;
  69
  70         ///cot(z)=cos(z)/sin(z)
  71
  72         if (icase != 5)
  73         {
  74
  75                 khd = Fourie.nsize1 + 2;
  76
  77                 for (kk = 1; kk < khd; kk += 2)
  78                 {     // bug79 - swr 2/11/05
  79                         pj2[kk] = -pj2[kk];
  80                 }
  81
  82                 ang2 += Consts.pi;
  83         }
  84
  85         uleft = (REAL) cot(ang2)
  86         ;
  87         uright = (REAL) cot(ang1)
  88         ;
  89         kh = Fourie.nsize2 / 2 - 1;
  90
  91         for (j = 1; j <= kh; j++)
  92         {     // bug79 - swr 2/11/05
  93
  94                 iy = j;
  95                 aleft = uleft * iy;
  96                 aright = uright * iy;
  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
  99
 100                 if (indl <= indh)
 101                 {
 102                         a = j * Fourie.delp / Fourie.deld;
 103                         low = (INTEGER) a;
 104                         indmin = 2 * low;     // bug79 - swr 2/11/05
 105                         indmax = indmin + 2;
 106                         distan = aright - aleft;
 107
 108                         switch (Fourie.interp)
 109                         {
 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)
 112                         case 1:
 113                         case 4:
 114                                 yy[1] = (REAL) low;
 115                                 yy[3] = (REAL) low;
 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;
 122                                 yg = a;
 123
 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);
 129                                 }
 130                                 break;
 131
 132                                 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 133                         case 2:
 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
 138                                         x = ix;
 139                                         if (x <= amid)
 140                                         {     // bug79 - swr 2/15/05
 141                                                 gr = pj2[ind - 1];
 142                                                 gi = pj2[ind];
 143                                         }
 144                                         else
 145                                         {
 146                                                 gr = pj1[ind - 1];
 147                                                 gi = pj1[ind];
 148                                         }
 149                                         store(ix, iy, gr, gi);
 150                                 }
 151                                 break;
 152
 153                                 // INTERP = 3 LINEAR INTERPOLATION
 154                         case 3:
 155                                 t = a - low;
 156                                 if (t <= Consts.zero)
 157                                         t = 0.0;
 158                                 if ((1. - t) <= Consts.zero)
 159                                         t = (REAL) 1.0;
 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;
 164
 165                                 for (ix = indl; ix <= indh; ix++)
 166                                 {     // bug79 - swr 2/11/05
 167                                         s = (aright - ix) / distan;
 168                                         if (s <= Consts.zero)
 169                                                 s = 0.0;
 170                                         if ((1. - s) <= Consts.zero)
 171                                                 s = 1.0;
 172                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 173                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 174                                         store(ix, iy, gr, gi);
 175                                 }
 176                                 break;
 177                         }
 178                 }
 179         }
 180
 181         if (icase == 5)
 182                 return;
 183
 184         khd = Fourie.nsize1 + 2;
 185         for (kk = 1; kk < khd; kk += 2)
 186         {     // bug79 - swr 2/11/05
 187                 pj2[kk] = -pj2[kk];
 188         }
 189         ang2 -= Consts.pi;
 190         return;
 191 }