src/snark/foru_prjtrn.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_prjtrn.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  PRJTRN: PROJECTION TRANSFORM
  11  PURPOSE: READ NP'TH PROJECTION INTO VECTOR, PUT THE ANGLE IN
  12  PHI (IN RADIAN ) AND DPHI ( IN DEGREE )  WHERE THE ANGLE
  13  IS ALWAYS BETWEEN -90 AND 90
  14  IF THE ANGLE IS NOT IN THAT RANGE ADD N*PI OR SUBSTRACT
  15  N*PI TO MAKE IT IN THE RANGE AND REVERSE THE PROJECTION
  16  DATA IF NECESSARY
  17  THEN TRANSFORM & FILTER THE DATA
  18  ***NOTE THE UNCONVENTIONAL DEFINITION OF ANGLE AND ORDER OF PROJ
  19  ***DATA IN SNARK
  20  */
  21
  22 #include <cstdio>
  23 #include <cmath>
  24
  25 #include "blkdta.h"
  26 #include "geom.h"
  27 #include "spctrm.h"
  28 #include "fourie.h"
  29 #include "consts.h"
  30 #include "uiod.h"
  31 #include "anglst.h"
  32 #include "projfile.h"
  33
  34 #include "foru.h"
  35
  36 void foru_class::prjtrn(INTEGER np, REAL* vector, REAL* phi, REAL* dphi)
  37 {
  38         REAL theta;
  39         REAL sinth;
  40         REAL costh;
  41         static REAL oldth;
  42
  43         REAL fact;
  44         INTEGER k;
  45         INTEGER invdir;
  46
  47         BOOLEAN revers;
  48
  49         Anglst.getang(np, &theta, &sinth, &costh);
  50         if (np != 0)
  51         {
  52                 if (oldth > theta)
  53                         error(3);
  54                 if (GeoPar.uni && ((theta - oldth) > Consts.pid2))
  55                         error(4);
  56                 if (GeoPar.vri && ((theta - oldth) > (Consts.pi / 4.0)))
  57                         error(4);
  58                 if (Fourie.iflg == 1)
  59                         return;
  60         }
  61         oldth = theta;
  62         *phi = theta - Consts.pid2;
  63
  64         while (*phi < -Consts.pid2)
  65         {
  66                 *phi += Consts.twopi;
  67         }
  68
  69         while (*phi > (Consts.pi + Consts.pid2))
  70         {
  71                 *phi -= Consts.twopi;
  72         }
  73
  74         revers = (*phi > Consts.pid2);
  75         if (revers)
  76                 *phi -= Consts.pi;
  77         *dphi = *phi * (REAL) 180.0 / Consts.pi;
  78
  79         ProjFile.ReadProj(np, Fourie.ntemp, Fourie.nsize1);
  80
  81         fact = Fourie.rinc;
  82         if (GeoPar.vri)
  83                 fact = Fourie.rinc * (REAL) MAX0(fabs(cos(*phi)), fabs(sin(*phi)));
  84         if (GeoPar.strip)
  85         {
  86                 for (k = 0; k < Fourie.nsize1; k++)
  87                 {
  88                         Fourie.ntemp[k] /= fact;
  89                 }
  90         }
  91
  92         ftmap(Fourie.ntemp, vector, revers);
  93
  94
  95         invdir = 1;
  96         fft(vector, Fourie.nforw, invdir, fact);
  97
  98 // Filtering
  99         frfilt(vector, Fourie.nsize1, Fourie.nfiltr, Fourie.cutoff, *phi);
 100         return;
 101 }