src/snark/bckprj.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/bckprj.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  SUBROUTINE TO BACK-PROJECT THE DATA GIVEN THE THE ARRAY 'G'
  11  ONTO THE SQUARE ARRAY 'A'.  G IS A VECTOR OF LENGTH M; A IS AN
  12  N*N MATRIX.  THE BACK PROJECTION IS PERFORMED USING THE
  13  INTERPOLATION METHOD GIVEN BY 'INTERP' PARALLEL TO THE UNIT
  14  VECTOR GIVEN BY (CTH,STH).
  15
  16
  17  This subroutine is modified so that it is not necessary to
  18  call qintp for linear iteration (two-point Lagrange).
  19  The exacution is done within this subroutine.
  20  July 23, 1989
  21  */
  22
  23 #include <cstdio>
  24
  25 #include "blkdta.h"
  26 #include "qintp.h"
  27
  28 #include "bckprj.h"
  29
  30 void bckprj(REAL* a, INTEGER n, REAL* g, INTEGER m, REAL sth, REAL cth,
  31                 REAL spac, INTEGER interp)
  32 {
  33         INTEGER epos;
  34         REAL nmid, mmid;
  35
  36         INTEGER ind;
  37         REAL yinc;
  38         REAL xinc;
  39         REAL ypos;
  40         REAL pos;
  41         INTEGER j;
  42         INTEGER i;
  43         REAL er;
  44         REAL qintp1;
  45
  46         nmid = (n - 1) / 2;
  47         mmid = (m + 1) / 2;
  48         ind = 0;
  49         yinc = cth * spac;
  50         xinc = -sth * spac;
  51         ypos = nmid * yinc + mmid;
  52
  53         if (interp == 2)
  54         {
  55                 for (i = 0; i < n; i++)
  56                 {
  57                         pos = -nmid * xinc + ypos;
  58                         for (j = 0; j < n; j++)
  59                         {
  60                                 epos = (INTEGER) pos;
  61                                 er = pos - epos;
  62                                 qintp1 = g[epos - 1] + er * (g[epos] - g[epos - 1]);
  63                                 a[ind++] += qintp1;
  64                                 pos += xinc;
  65                         }
  66                         ypos -= yinc;
  67                 }
  68         }
  69         else
  70         {
  71                 for (i = 0; i < n; i++)
  72                 {
  73                         pos = -nmid * xinc + ypos;
  74                         for (j = 0; j < n; j++)
  75                         {
  76                                 a[ind++] += qintp(pos, g, m, interp);
  77                                 pos += xinc;
  78                         }
  79                         ypos -= yinc;
  80                 }
  81         }
  82         return;
  83 }