src/snark/foru_uint.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_uint.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10
  11  PURPOSE: GIVEN THE COORDINATES OF THE DESIRED POINT  (IX,IY)
  12  AND THE TWO TRANSFORMED PROJECTIONS  PJ1 AND PJ2 FIND THE
  13  INTERPOLATED VALUE (GR,GI) THEN CALL THE SUBROUTINE STORE
  14  TO STORE THE VALUE IN APPROPRIATE PLACE IN THE ARRAY
  15  STARTING AT NFRPLN+1
  16  THE INTERPOLATION SCHEME TO BE USED IS DETERMINED BY
  17  INTERP (SEE ALSO COMMENTS IN THE MAIN PROGRAM )
  18  */
  19
  20 #include <cstdio>
  21 #include <cmath>
  22
  23 #include "blkdta.h"
  24 #include "fourie.h"
  25 #include "consts.h"
  26 #include "uiod.h"
  27
  28 #include "foru.h"
  29
  30 #define disn(b,c,the)  ((REAL) sqrt((b)*(b)+(c)*(c)-2.*(b)*(c) * (REAL) cos(the)))
  31
  32 void foru_class::uint1(INTEGER icase, INTEGER ix, INTEGER iy, REAL ang1, REAL ang2, REAL* pj1, REAL* pj2)
  33 {
  34
  35         INTEGER x;
  36         INTEGER y;
  37         REAL r;
  38         REAL ang;
  39         REAL a;
  40         INTEGER low;
  41         INTEGER indmin;
  42         INTEGER indmax;
  43         REAL s;
  44         REAL t;
  45         INTEGER ind;
  46         REAL gr;
  47         REAL gi;
  48         REAL pr;
  49         REAL qr;
  50         REAL px;
  51         REAL qi;
  52         REAL ak;
  53         REAL d1;
  54         REAL d2;
  55         REAL d3;
  56         REAL d4;
  57         REAL e1;
  58         REAL e2;
  59         REAL e3;
  60         REAL e4;
  61         REAL den;
  62
  63         x = ix;
  64         y = iy;
  65         r = (REAL) sqrt((REAL) x * x + y * y) * Fourie.delp;
  66         ang = (REAL) atan2((REAL) y, (REAL) x);
  67         a = r / Fourie.deld;
  68         if (a >= Fourie.radmax)
  69                 return;
  70         low = (INTEGER) a;
  71         indmin = 2 * low;
  72         indmax = indmin + 2;
  73         s = (REAL) fabs((ang - ang1) / (ang2 - ang1));
  74         t = a - low;
  75
  76         switch (Fourie.interp)
  77         {
  78
  79         // INTERP = 1 NEAREST NEIGHBOR IN CARTESIAN SENSE
  80         case 1:
  81                 if (s <= 0.5)
  82                 {
  83                         ind = indmax;
  84                         if ((a * cos(ang - ang1) - low) <= 0.5)
  85                                 ind = indmin;
  86                         gr = pj1[ind];
  87                         gi = pj1[ind + 1];
  88                 }
  89                 else
  90                 {
  91                         ind = indmax;
  92                         if ((a * cos(ang2 - ang) - low) <= 0.5)
  93                                 ind = indmin;
  94                         gr = pj2[ind];
  95                         gi = pj2[ind + 1];
  96                 }
  97                 break;
  98
  99                 // INTERP = 2 NEAREST NEIGHBOR IN RADIAL SENSE
 100         case 2:
 101                 if (t <= 0.5)
 102                         ind = indmin;
 103                 if (t > 0.5)
 104                         ind = indmax;
 105                 if (s <= 0.5)
 106                 {
 107                         gr = pj1[ind];
 108                         gi = pj1[ind + 1];
 109                 }
 110                 else
 111                 {
 112                         gr = pj2[ind];
 113                         gi = pj2[ind + 1];
 114                 }
 115                 break;
 116
 117                 // INTERP = 3 LINEAR INTERPOLATION
 118         case 3:
 119                 if (t <= Consts.zero)
 120                         t = 0.0;
 121                 if ((1.0 - t) <= Consts.zero)
 122                         t = 1.0;
 123                 if (s <= Consts.zero)
 124                         s = 0.0;
 125                 if ((1.0 - s) <= Consts.zero)
 126                         s = (REAL) 1.0;
 127                 pr = t * pj1[indmax] + ((REAL) 1.0 - t) * pj1[indmin];
 128                 qr = t * pj2[indmax] + ((REAL) 1.0 - t) * pj2[indmin];
 129                 gr = pr * ((REAL) 1.0 - s) + qr * s;
 130                 px = t * pj1[indmax + 1] + ((REAL) 1.0 - t) * pj1[indmin + 1];
 131                 qi = t * pj2[indmax + 1] + ((REAL) 1.0 - t) * pj2[indmin + 1];
 132                 gi = px * ((REAL) 1.0 - s) + qi * s;
 133                 break;
 134
 135                 // INTERP = 4 F=(F1/D1+F2/D2+F3/D3+F4/D4)/(1/D1+1/D2+1/D3+1/D4)
 136         case 4:
 137                 ak = low;
 138                 d1 = disn(a, ak, ang - ang1);
 139                 d2 = disn(a, ak + 1, ang - ang1);
 140                 d3 = disn(a, ak, ang2 - ang);
 141                 d4 = disn(a, ak + 1, ang2 - ang);
 142                 e1 = d2 * d3 * d4;
 143                 e2 = d1 * d3 * d4;
 144                 e3 = d1 * d2 * d4;
 145                 e4 = d1 * d2 * d3;
 146
 147
 148                 if (e1 <= Consts.zero)
 149                         e1 = 0;
 150                 if (e2 <= Consts.zero)
 151                         e2 = 0;
 152                 if (e3 <= Consts.zero)
 153                         e3 = 0;
 154                 if (e4 <= Consts.zero)
 155                         e4 = 0;
 156
 157                 den = e1 + e2 + e3 + e4;
 158
 159
 160                 gr = (pj1[indmin] * e1 + pj1[indmax] * e2 + pj2[indmin] * e3
 161                                 + pj2[indmax] * e4) / den;
 162
 163                 gi = (pj1[indmin + 1] * e1 + pj1[indmax + 1] * e2 + pj2[indmin + 1] * e3
 164                                 + pj2[indmax + 1] * e4) / den;
 165         }
 166
 167         // CALL THE SUBROUTINE STORE TO STORE THE VALUE (GR,GI) AT POINT
 168         // (IX,IY) IN APPROPRIATE PLACE IN THE ARRAY STARTING AT NFRPLN+1
 169
 170         store(ix, iy, gr, gi);
 171         return;
 172 }