src/snark/foru_unif.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_unif.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  PURPOSE: FOR EACH CASE (ICASE=1,...,6) FIND ALL POINTS LIE BETWEEN
  11  THE TWO RAYS AND CALL UINT TO CALAULATE INTERPOLATED
  12  VALUE
  13  */
  14
  15 #include <cstdio>
  16 #include <cmath>
  17
  18 #include "blkdta.h"
  19 #include "fourie.h"
  20 #include "consts.h"
  21
  22 #include "foru.h"
  23
  24 #define  cot(z)  (cos(z)/sin(z));
  25
  26 void foru_class::unif(INTEGER icase, REAL ang1, REAL ang2, REAL* nold, REAL* new1)
  27 {
  28         INTEGER kh, ix, iy, j, khd, kk;
  29         REAL uleft, uright, uhigh, ulow, xhigh, xlow, yhigh, ylow;
  30
  31         kh = Fourie.nsize2 / 2 - 1;
  32
  33         switch (icase)
  34         {
  35
  36 //.... CASE 1
  37         case 1:
  38                 uleft = (REAL) cot(ang1)
  39                 ;
  40                 uright = (REAL) cot(ang2)
  41                 ;
  42                 for (j = 1; j <= kh; j++)
  43                 {
  44                         iy = -j;
  45                         xlow = uleft * iy - Fourie.epsln;
  46                         xhigh = uright * iy + Fourie.epsln;
  47                         ix = lround(xlow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
  48                         while (ix <= xhigh)
  49                         {
  50                                 if (ix >= xlow)
  51                                 {
  52                                         uint1(icase, ix, iy, ang1, ang2, nold, new1);
  53                                 }
  54                                 ix++;
  55                         }
  56                 }
  57                 break;
  58
  59 //.... CASE 2
  60         case 2:
  61                 uleft = (REAL) cot(ang1)
  62                 ;
  63                 for (j = 1; j <= kh; j++)
  64                 {
  65                         iy = -j;
  66                         xlow = uleft * iy - Fourie.epsln;
  67                         ix = lround(xlow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
  68                         while (ix <= j)
  69                         {
  70                                 if (ix >= xlow)
  71                                 {
  72                                         uint1(icase, ix, iy, ang1, ang2, nold, new1);
  73                                 }
  74                                 ix++;
  75                         }
  76                 }
  77
  78                 uhigh = (REAL) tan(ang2);
  79                 for (ix = 1; ix <= kh; ix++)
  80                 {
  81                         yhigh = ix * uhigh + Fourie.epsln;
  82                         iy = -ix + 1;
  83                         while (iy <= yhigh)
  84                         {
  85                                 uint1(icase, ix, iy, ang1, ang2, nold, new1);
  86                                 iy++;
  87                         }
  88                 }
  89                 break;
  90
  91                 // CASE 3
  92         case 3:
  93                 ulow = (REAL) tan(ang1);
  94                 uhigh = (REAL) tan(ang2);
  95                 for (ix = 1; ix <= kh; ix++)
  96                 {
  97                         ylow = ulow * ix - Fourie.epsln;
  98                         yhigh = uhigh * ix + Fourie.epsln;
  99                         iy = lround(ylow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 100                         while (iy <= yhigh)
 101                         {
 102                                 if (iy >= ylow)
 103                                 {
 104                                         uint1(icase, ix, iy, ang1, ang2, nold, new1);
 105                                 }
 106                                 iy++;
 107                         }
 108                 }
 109                 break;
 110
 111                 // CASE 4
 112
 113         case 4:
 114                 ulow = (REAL) tan(ang1);
 115                 for (ix = 1; ix <= kh; ix++)
 116                 {  // !!! ix = 0 does not work mk
 117                         ylow = ulow * ix - Fourie.epsln;
 118                         iy = lround(ylow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 119                         while (iy <= ix)
 120                         { // changed ">=" to "<=" to match Fortran code. Lajos, Jan 13, 2005
 121                                 if (iy >= ylow)
 122                                 {
 123                                         uint1(icase, ix, iy, ang1, ang2, nold, new1);
 124                                 }
 125                                 iy++;
 126                         }
 127                 }
 128
 129                 uleft = (REAL) cot(ang2)
 130                 ;
 131
 132                 for (iy = 1; iy <= kh; iy++)
 133                 {
 134                         xlow = uleft * iy - Fourie.epsln;
 135                         ix = iy - 1;
 136                         while (ix >= xlow)
 137                         {
 138                                 uint1(icase, ix, iy, ang1, ang2, nold, new1);
 139                                 ix--;
 140                         }
 141                 }
 142
 143                 break;
 144
 145                 // CASE 5
 146
 147         case 5:
 148                 uleft = (REAL) cot(ang2)
 149                 ;
 150                 uright = (REAL) cot(ang1)
 151                 ;
 152                 for (iy = 1; iy <= kh; iy++)
 153                 {
 154                         xlow = uleft * iy - Fourie.epsln;
 155                         xhigh = uright * iy + Fourie.epsln;
 156                         ix = lround(xlow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 157                         while (ix <= xhigh)
 158                         {
 159                                 if (ix >= xlow)
 160                                 {
 161                                         uint1(icase, ix, iy, ang1, ang2, nold, new1);
 162                                 }
 163                                 ix++;
 164                         }
 165                 }
 166                 break;
 167
 168                 // CASE 6
 169
 170                 // CONJUGATE THE SECOND PROJECTION
 171         case 6:
 172                 khd = Fourie.nsize1 + 2;
 173                 for (kk = 1; kk < khd; kk += 2)
 174                 {
 175                         new1[kk] = -new1[kk];
 176                 }
 177
 178                 uleft = (REAL) cot(ang2)
 179                 ;
 180                 uright = (REAL) cot(ang1)
 181                 ;
 182                 for (iy = 1; iy <= kh; iy++)
 183                 {
 184                         xlow = uleft * iy - Fourie.epsln;
 185                         xhigh = uright * iy + Fourie.epsln;
 186                         ix = lround(xlow); // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 187                         while (ix <= xhigh)
 188                         {
 189                                 if (ix >= xlow)
 190                                 {
 191                                         uint1(icase, ix, iy, ang1, ang2 + Consts.pi, nold, new1);
 192                                 }
 193                                 ix++;
 194                         }
 195                 }
 196
 197                 // CHANGE BACK
 198
 199                 khd = Fourie.nsize1 + 2;
 200                 for (kk = 1; kk < khd; kk += 2)
 201                 {
 202                         new1[kk] = -new1[kk];
 203                 }
 204                 break;
 205         }
 206
 207         return;
 208 }