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