src/snark/qintp.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/qintp.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  INITIALIZE CONVOLUTING FUNCTION FOR CONVOLUTION ALGORITHM.
  11  */
  12
  13 #include <cstdio>
  14 #include <cmath>
  15
  16 #include "blkdta.h"
  17 #include "consts.h"
  18
  19 #include "qintp.h"
  20
  21 REAL qintp(REAL pos, REAL* table, INTEGER n, INTEGER interp)
  22 {
  23         INTEGER epos;
  24         INTEGER opos;
  25         REAL er;
  26         REAL or1;
  27         INTEGER i;
  28
  29         REAL a0;
  30         REAL a1;
  31         REAL a2;
  32         REAL a3;
  33         REAL sum;
  34         REAL sgn;
  35         INTEGER k;
  36         REAL a1t2;
  37         REAL a2t2;
  38         REAL a4;
  39         REAL a0t2;
  40         REAL a1t4;
  41         REAL a2t4;
  42         REAL a3t4;
  43         REAL a4t2;
  44         REAL p1pm2;
  45         REAL p2pm3;
  46         REAL ppm1;
  47         REAL pa23;
  48         REAL pa14;
  49         REAL pa05;
  50
  51         epos = (INTEGER) pos;
  52         opos = (INTEGER) (pos + 0.5);
  53         er = pos - epos;
  54         or1 = pos - opos;
  55         i = interp + 2;
  56
  57         /// fix for indexes in C++
  58         epos--;
  59         opos--;
  60
  61         switch (i)
  62         {
  63
  64         // MODIFIED CUBIC SPLINE
  65
  66         case 1:
  67                 a0 = table[epos - 1];
  68                 a1 = table[epos];
  69                 a2 = table[epos + 1];
  70                 a3 = table[epos + 2];
  71                 return a1 + (REAL) 0.5 * er * (a2 - a0 + er     * ((REAL) 2.0 * a0 - (REAL) 5.0 * a1 + (REAL) 4.0 * a2 - a3     + er * (a3 - a0 + (REAL) 3.0 * (a1 - a2))));
  72
  73                 // BAND LIMITING (SINC) INTERPOLATION
  74
  75         case 2:
  76                 if (fabs(or1) < Consts.zero)
  77                 {
  78                         return table[opos];
  79                 }
  80                 sum = 0.0;
  81                 sgn = -1.0;
  82                 for (k = 0; k < n; k++)
  83                 {
  84                         sum += sgn * table[k] / (pos - k);
  85                         sgn = -sgn;
  86                 }
  87
  88                 if (((opos / 2) * 2) == opos)
  89                 {
  90                         return sum * (REAL) sin(Consts.pi * or1) / Consts.pi;
  91                 }
  92                 else
  93                 {
  94                         return -sum * (REAL) sin(Consts.pi * or1) / Consts.pi;
  95                 }
  96
  97                 // NEAREST POINT
  98
  99         case 3:
 100                 return table[opos];
 101
 102                 // LINEAR INTERPOLATION
 103
 104         case 4:
 105                 return table[epos] + er * (table[epos + 1] - table[epos]);
 106
 107                 // 3 POINT LAGRANGE
 108
 109         case 5:
 110                 a0 = table[opos - 1] / (REAL) 2.0;
 111                 a1 = table[opos];
 112                 a2 = table[opos + 1] / (REAL) 2.0;
 113                 return a1 + or1 * (-a0 + a2 + or1 * (a0 - a1 + a2));
 114
 115                 // 4 POINT LAGRANGE
 116
 117         case 6:
 118                 a0 = table[epos - 1] / (REAL) 6.0;
 119                 a1 = table[epos] / (REAL) 2.0;
 120                 a2 = table[epos + 1] / (REAL) 2.0;
 121                 a3 = table[epos + 2] / (REAL) 6.0;
 122                 a1t2 = table[epos];
 123                 a2t2 = table[epos + 1];
 124                 return a1t2 + er * ((REAL) -2.0 * a0 - a1 + a2t2 - a3 + er * ((REAL) +3.0 * a0 - a1t2 + a2 + er * (-a0 + a1 - a2 + a3)));
 125
 126                 // 5 POINT LAGRANGE
 127
 128         case 7:
 129                 a0 = table[opos - 2] / (REAL) 24.0;
 130                 a1 = table[opos - 1] / (REAL) 6.0;
 131                 a2 = table[opos] / (REAL) 4.0;
 132                 a3 = table[opos + 1] / (REAL) 6.0;
 133                 a4 = table[opos + 2] / (REAL) 24.0;
 134                 a0t2 = (REAL) 2.0 * a0;
 135                 a1t4 = (REAL) 4.0 * a1;
 136                 a2t4 = table[opos];
 137                 a3t4 = (REAL) 4.0 * a3;
 138                 a4t2 = (REAL) 2.0 * a4;
 139                 return a2t4 + or1 * (a0t2 - a1t4 + a3t4 - a4t2 + or1 * (-a0 + a1t4 - (REAL) 5.0 * a2 + a3t4 - a4 + or1 * (-a0t2 + a1 - a3 + a4t2 + or1 * (a0 - a1 + a2 - a3 + a4))));
 140
 141                 // 6 POINT LAGRANGE
 142
 143         case 8:
 144                 p1pm2 = (er + (REAL) 1.0) * (er - (REAL) 2.0);
 145                 p2pm3 = (er + (REAL) 2.0) * (er - (REAL) 3.0);
 146                 ppm1 = er * (er - (REAL) 1.0);
 147                 pa23 = p1pm2 * p2pm3 / (REAL) 12.0;
 148                 pa14 = p2pm3 * ppm1 / (REAL) 24.0;
 149                 pa05 = ppm1 * p1pm2 / (REAL) 120.0;
 150                 return -pa05 * (er - (REAL) 3.0) * table[epos - 2]
 151                                 + pa14 * (er - (REAL) 2.0) * table[epos - 1]
 152                                 - pa23 * (er - (REAL) 1.0) * table[epos]
 153                                 + pa23 * er * table[epos + 1]
 154                                 - pa14 * (er + (REAL) 1.0) * table[epos + 2]
 155                                 + pa05 * (er + (REAL) 2.0) * table[epos + 3];
 156
 157         default:
 158                 // error invalid parameter interp
 159                 return 0;
 160         }
 161 }