src/snark/foru_vcase4.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_vcase4.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  VCASE4: VARIABLE SPACING CASE 4
  13  PURPOSE: FOR VARIABLE SPACING, ICASE=4 ( THE TWO LINES STRADLE
  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::vcase4(REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
  34 {
  35         REAL xx[4];
  36         REAL yy[4];
  37
  38         REAL ulow;
  39         REAL uleft;
  40         INTEGER kh;
  41         INTEGER j;
  42         INTEGER ix;
  43         REAL a;
  44         INTEGER low;
  45         INTEGER indmin;
  46         INTEGER indmax;
  47         REAL alow;
  48         REAL aleft;
  49         REAL vdis;
  50         REAL hdis;
  51         REAL distan;
  52         INTEGER indl;
  53         INTEGER indh;
  54         REAL xg;
  55         INTEGER iy;
  56         REAL yg;
  57         REAL gr;
  58         REAL gi;
  59         REAL amid;
  60         INTEGER ind;
  61         INTEGER y;
  62         REAL t;
  63         REAL pr;
  64         REAL qr;
  65         REAL px;
  66         REAL qi;
  67         REAL s;
  68         INTEGER x;
  69
  70         ///tan(z)=sin(z)/cos(z)
  71         ///cot(z)=cos(z)/sin(z)
  72
  73         ulow = (REAL) tan(ang1);
  74         uleft = (REAL) cot(ang2)
  75         ;
  76         kh = Fourie.nsize2 / 2 - 1;
  77
  78         for (j = 1; j <= kh; j++)
  79         {      // bug79 - swr 2/11/05
  80
  81                 ix = j;
  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                 alow = ulow * ix;
  87                 aleft = uleft * j;
  88                 vdis = j - alow;
  89                 hdis = j - aleft;
  90                 distan = hdis + vdis;
  91                 indl = (INTEGER) ceil(alow - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
  92                 indh = j - 1;
  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                         xx[1] = (REAL) low;
 102                         xx[2] = (REAL) low + 1;
 103                         yy[1] = xx[1] * ulow;
 104                         yy[2] = xx[2] * ulow;
 105                         yy[3] = (REAL) low;
 106                         yy[4] = (REAL) low + 1;
 107                         xx[3] = yy[3] * uleft;
 108                         xx[4] = yy[4] * uleft;
 109                         xg = a;
 110
 111                         if (indl <= indh)
 112                         {
 113                                 for (iy = indl; iy <= indh; iy++)
 114                                 {     // bug79 - swr 2/11/05
 115                                         yg = iy * Fourie.delp / Fourie.deld;
 116                                         vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
 117                                         store(ix, iy, gr, gi);
 118                                 }
 119                         }
 120                         break;
 121
 122                         // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 123                 case 2:
 124
 125                         amid = alow + distan / (REAL) 2.0;
 126                         ind = 2 * (lround(a)) + 1; // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 127
 128                         if (indl <= indh)
 129                         {
 130                                 for (iy = indl; iy <= indh; iy++)
 131                                 {     // bug79 - swr 2/11/05
 132                                         y = iy;
 133                                         if (y <= amid)
 134                                         {     // bug79 - swr 2/15/05
 135                                                 gr = pj1[ind - 1];
 136                                                 gi = pj1[ind];
 137                                         }
 138                                         else
 139                                         {
 140                                                 gr = pj2[ind - 1];
 141                                                 gi = pj2[ind];
 142                                         }
 143                                         store(ix, iy, gr, gi);
 144                                 }
 145                         }
 146                         break;
 147
 148                         // INTERP = 3 LINEAR INTERPOLATION
 149                 case 3:
 150                         t = a - low;
 151                         if (t <= Consts.zero)
 152                                 t = 0.0;
 153                         if ((1. - t) <= Consts.zero)
 154                                 t = 1.0;
 155                         pr = pj1[indmin] * ((REAL) 1.0 - t) + pj1[indmax] * t;
 156                         qr = pj2[indmin] * ((REAL) 1.0 - t) + pj2[indmax] * t;
 157                         px = pj1[indmin + 1] * ((REAL) 1.0 - t) + pj1[indmax + 1] * t;
 158                         qi = pj2[indmin + 1] * ((REAL) 1.0 - t) + pj2[indmax + 1] * t;
 159
 160                         if (indl <= indh)
 161                         {
 162
 163                                 for (iy = indl; iy <= indh; iy++)
 164                                 {     // bug79 - swr 2/11/05
 165                                         s = (iy - alow) / distan;
 166                                         if (s <= Consts.zero)
 167                                                 s = 0.0;
 168                                         if (((REAL) 1.0 - s) <= Consts.zero)
 169                                                 s = (REAL) 1.0;
 170                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 171                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 172                                         store(ix, iy, gr, gi);
 173                                 }
 174                         }
 175                         break;
 176                 }
 177
 178                 iy = j;
 179                 indl = (INTEGER) ceil(aleft - Fourie.epsln); // bug82 - changed to ceil - swr - 2/20/05
 180                 indh = j;
 181
 182                 if (indl <= indh)
 183                 {
 184
 185                         switch (Fourie.interp)
 186                         {
 187
 188                         // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
 189                         // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
 190                         case 1:
 191                         case 4:
 192                                 yg = a;
 193                                 for (ix = indl; ix <= indh; ix++)
 194                                 {     // bug79 - swr 2/11/05
 195                                         xg = ix * Fourie.delp / Fourie.deld;
 196                                         vint14(xx, yy, xg, yg, pj1, pj2, indmin, indmax, &gr, &gi); // changed to pass address of "gr" and "gi". Lajos, Feb 10, 2005
 197                                         store(ix, iy, gr, gi);
 198                                 }
 199                                 break;
 200
 201                                 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 202                         case 2:
 203                                 amid = aleft + distan / (REAL) 2.0;
 204                                 for (ix = indl; ix <= indh; ix++)
 205                                 {     // bug79 - swr 2/11/05
 206                                         x = ix;
 207                                         if (x <= amid)
 208                                         {     // bug79 - swr 2/15/05
 209                                                 gr = pj2[ind - 1];
 210                                                 gi = pj2[ind];
 211                                         }
 212                                         else
 213                                         {
 214                                                 gr = pj1[ind - 1];
 215                                                 gi = pj1[ind];
 216                                         }
 217                                         store(ix, iy, gr, gi);
 218                                 }
 219                                 break;
 220
 221                                 // INTERP = 3 LINEAR INTERPOLATION
 222                         case 3:
 223                                 for (ix = indl; ix <= indh; ix++)
 224                                 {     // bug79 - swr 2/11/05
 225                                         s = (vdis + j - ix) / distan;
 226                                         if (s <= Consts.zero)
 227                                                 s = 0.0;
 228                                         if ((1. - s) <= Consts.zero)
 229                                                 s = 1.0;
 230                                         gr = pr * ((REAL) 1.0 - s) + qr * s;
 231                                         gi = px * ((REAL) 1.0 - s) + qi * s;
 232                                         store(ix, iy, gr, gi);
 233                                 }
 234                                 break;
 235                         }
 236                 }
 237         }
 238         return;
 239 }