src/snark/foru_vcase3.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_vcase3.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ------------------------------------------------------------------
  11
  12  VCASE3: VARIABLE SPACING CASE 3
  13  PURPOSE: FOR VARIABLE SPACING, ICASE=3 (THE TWO LINES LIE BETWEEN
  14  -45 AND 45 DEGREE LINES) 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
  18
  19  ------------------------------------------------------------------
  20  */
  21 // bug79 - problem with fourier algorithm
  22 #include <cstdio>
  23 #include <cmath>
  24
  25 #include "blkdta.h"
  26 #include "fourie.h"
  27 #include "consts.h"
  28
  29 #include "foru.h"
  30
  31 void foru_class::vcase3(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
  32 {
  33         REAL xx[4];
  34         REAL yy[4];
  35
  36         REAL ulow;
  37         REAL uhigh;
  38         INTEGER kh;
  39         INTEGER j;
  40         INTEGER ix;
  41         REAL alow;
  42         REAL ahigh;
  43         INTEGER indl;
  44         INTEGER indh;
  45         REAL a;
  46         INTEGER low;
  47         INTEGER indmin;
  48         INTEGER indmax;
  49         REAL distan;
  50         REAL xg;
  51         INTEGER iy;
  52         REAL yg;
  53         REAL gr;
  54         REAL gi;
  55         REAL amid;
  56         INTEGER y;
  57         INTEGER ind;
  58
  59         REAL t;
  60         REAL pr;
  61         REAL qr;
  62         REAL px;
  63         REAL qi;
  64         REAL s;
  65
  66         //tan(z)=sin(z)/cos(z)
  67
  68         ulow = (REAL) tan(ang1);
  69         uhigh = (REAL) tan(ang2);
  70         kh = Fourie.nsize2 / 2 - 1;
  71
  72         for (j = 1; j <= kh; j++)
  73         {     // bug79 - swr 2/11/05
  74
  75                 ix = j;
  76                 alow = ulow * ix;
  77                 ahigh = uhigh * ix;
  78                 indl = (INTEGER) ceil(alow - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
  79                 indh = (INTEGER) floor(ahigh + Fourie.epsln); // bug82 - changed to floor - swr - 2/20/05
  80                 if (indl <= indh)
  81                 {
  82                         a = j * Fourie.delp / Fourie.deld;
  83                         low = (INTEGER) a;
  84                         indmin = 2 * low;     // bug79 - swr 2/11/05
  85                         indmax = indmin + 2;
  86                         distan = ahigh - alow;
  87
  88                         switch (Fourie.interp)
  89                         {
  90
  91                         // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
  92                         // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
  93
  94                         case 1:
  95                         case 4:
  96                                 xx[1] = (REAL) low;
  97                                 xx[2] = (REAL) low + 1;
  98                                 xx[3] = (REAL) low;
  99                                 xx[4] = (REAL) low + 1;
 100                                 yy[1] = xx[1] * ulow;
 101                                 yy[2] = xx[2] * ulow;
 102                                 yy[3] = xx[3] * uhigh;
 103                                 yy[4] = xx[4] * uhigh;
 104                                 xg = a;
 105                                 for (iy = indl; iy <= indh; iy++)
 106                                 {     // bug79 - swr 2/11/05
 107                                         yg = iy * Fourie.delp / Fourie.deld;
 108                                         vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
 109                                         store(ix, iy, gr, gi);
 110                                 }
 111
 112                                 break;
 113
 114                                 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 115
 116                         case 2:
 117                                 amid = (alow + ahigh) / (REAL) 2.0;
 118                                 for (iy = indl; iy <= indh; iy++)
 119                                 {     // bug79 - swr 2/11/05
 120                                         y = iy;
 121                                         ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 122
 123                                         if (y <= amid)
 124                                         {     // bug79 - swr 2/15/05
 125                                                 gr = pj1[ind - 1];
 126                                                 gi = pj1[ind];
 127                                         }
 128                                         else
 129                                         {
 130                                                 gr = pj2[ind - 1];
 131                                                 gi = pj2[ind];
 132                                         }
 133                                         store(ix, iy, gr, gi);
 134                                 }
 135                                 break;
 136
 137                                 // INTERP = 3 LINEAR INTERPOLATION
 138
 139                         case 3:
 140                                 t = a - low;
 141                                 if (t <= Consts.zero)
 142                                         t = 0.0;
 143                                 if (((REAL) 1.0 - t) <= Consts.zero)
 144                                         t = (REAL) 1.0;
 145                                 pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
 146                                 qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
 147                                 px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
 148                                 qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
 149
 150                                 for (iy = indl; iy <= indh; iy++)
 151                                 {     // bug79 - swr 2/11/05
 152                                         s = (iy - alow) / distan;
 153                                         if (s <= Consts.zero)
 154                                                 s = 0.0;
 155                                         if (((REAL) 1.0 - s) <= Consts.zero)
 156                                                 s = (REAL) 1.0;
 157                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 158                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 159                                         store(ix, iy, gr, gi);
 160
 161                                 }
 162                                 break;
 163
 164                         }
 165                 }
 166         }
 167         return;
 168 }