src/snark/foru_vcase2.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_vcase2.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  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
  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 #define  cot(z)  (cos(z)/sin(z));
  32
  33 void foru_class::vcase2(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
  34 {
  35         REAL xx[4];
  36         REAL yy[4];
  37
  38         REAL uleft;
  39         REAL uhigh;
  40         INTEGER kh;
  41         INTEGER j;
  42         INTEGER iy;
  43         REAL aleft;
  44         REAL ahigh;
  45         INTEGER indl;
  46         INTEGER indh;
  47         REAL a;
  48         INTEGER low;
  49         INTEGER indmin;
  50         INTEGER indmax;
  51         REAL hdis;
  52         REAL vdis;
  53         REAL distan;
  54         REAL yg;
  55         INTEGER ix;
  56         REAL xg;
  57         REAL gr;
  58         REAL gi;
  59         REAL amid;
  60         INTEGER ind;
  61         INTEGER x;
  62         REAL t;
  63         REAL pr;
  64         REAL qr;
  65         REAL px;
  66         REAL qi;
  67         REAL s;
  68         INTEGER y;
  69
  70         //tan(z)=sin(z)/cos(z)
  71         ///cot(z)=cos(z)/sin(z)
  72
  73         uleft = (REAL) cot(ang1)
  74         ;
  75         uhigh = (REAL) tan(ang2);
  76         kh = Fourie.nsize2 / 2 - 1;
  77
  78         for (j = 1; j <= kh; j++)
  79         {     // bug79 - swr 2/11/05
  80
  81                 iy = -j;
  82                 aleft = uleft * iy;
  83                 ahigh = uhigh * j;
  84                 indl = (INTEGER) ceil(aleft - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
  85                 indh = j - 1;
  86                 a = j * Fourie.delp / Fourie.deld;
  87                 low = (INTEGER) a;
  88                 indmin = 2 * low;     // bug79 - swr 2/11/05
  89                 indmax = indmin + 2;
  90                 hdis = j - aleft;
  91                 vdis = ahigh - (-j);
  92                 distan = hdis + vdis;
  93
  94                 switch (Fourie.interp)
  95                 {
  96
  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)
  99                 case 1:
 100                 case 4:
 101                         yy[1] = (REAL) -low;
 102                         yy[2] = (REAL) -low - 1;
 103                         xx[1] = yy[1] * uleft;
 104                         xx[2] = yy[2] * uleft;
 105                         xx[3] = (REAL) low;
 106                         xx[4] = (REAL) low + 1;
 107                         yy[3] = xx[3] * uhigh;
 108                         yy[4] = xx[4] * uhigh;
 109                         yg = -a;
 110                         if (indl <= indh)
 111                         {
 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);
 117                                 }
 118                         }
 119                         break;
 120
 121                         // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 122                 case 2:
 123
 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
 126                         if (indl <= indh)
 127                         {
 128                                 for (ix = indl; ix <= indh; ix++)
 129                                 {     // bug79 - swr 2/11/05
 130                                         x = ix;
 131                                         if (x <= amid)
 132                                         {     // bug79 - swr 2/15/05
 133                                                 gr = pj1[ind - 1];
 134                                                 gi = pj1[ind];
 135                                         }
 136                                         else
 137                                         {
 138                                                 gr = pj2[ind - 1];
 139                                                 gi = pj2[ind];
 140                                         }
 141                                         store(ix, iy, gr, gi);
 142                                 }
 143                         }
 144                         break;
 145
 146                         // INTERP = 3 LINEAR INTERPOLATION
 147                 case 3:
 148                         t = a - low;
 149                         if (t <= Consts.zero)
 150                                 t = 0.0;
 151                         if ((1. - t) <= Consts.zero)
 152                                 t = 1.0;
 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;
 157                         if (indl <= indh)
 158                         {
 159                                 for (ix = indl; ix <= indh; ix++)
 160                                 {     // bug79 - swr 2/11/05
 161                                         s = (ix - aleft) / distan;
 162                                         if (s <= Consts.zero)
 163                                                 s = 0.0;
 164                                         if ((1. - s) <= Consts.zero)
 165                                                 s = 1.0;
 166                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 167                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 168                                         store(ix, iy, gr, gi);
 169                                 }
 170                         }
 171                         break;
 172
 173                 }
 174
 175                 ix = j;
 176                 indl = -j;
 177                 indh = (INTEGER) floor(ahigh + Fourie.epsln); // bug82 - changed to floor - swr - 2/20/05
 178
 179                 if (indl <= indh)
 180                 {
 181
 182                         switch (Fourie.interp)
 183                         {
 184
 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)
 187                         case 1:
 188                         case 4:
 189                                 xg = a;
 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);
 195                                 }
 196                                 break;
 197
 198                                 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 199                         case 2:
 200                                 amid = ahigh - distan / (REAL) 2.0;
 201                                 for (iy = indl; iy <= indh; iy++)
 202                                 {     // bug79 - swr 2/11/05
 203                                         y = iy;
 204                                         if (y <= amid)
 205                                         {     // bug79 - swr 2/15/05
 206                                                 gr = pj1[ind - 1];
 207                                                 gi = pj1[ind];
 208                                         }
 209                                         else
 210                                         {
 211                                                 gr = pj2[ind - 1];
 212                                                 gi = pj2[ind];
 213                                         }
 214                                         store(ix, iy, gr, gi);
 215                                 }
 216                                 break;
 217
 218                                 // INTERP = 3 LINEAR INTERPOLATION
 219                         case 3:
 220                                 for (iy = indl; iy <= indh; iy++)
 221                                 {     // bug79 - swr 2/11/05
 222                                         s = (hdis + iy - (-j)) / distan;
 223                                         if (s <= Consts.zero)
 224                                                 s = 0.0;
 225                                         if ((1. - s) <= Consts.zero)
 226                                                 s = 1.0;
 227                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 228                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 229                                         store(ix, iy, gr, gi);
 230                                 }
 231                                 break;
 232                         }
 233                 }
 234         }
 235 }