src/snark/congeo.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/congeo.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  **********************************************************************
   9
  10  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  11
  12  THIS SUBROUTINE TRANSFORMS ONE DIVERGENT BEAM OF
  13  PROJECTION DATA
  14  PROJEC(NR),NR=1,2,...,NRAYS,
  15  WITH DETECTOR SPACING PINC ( AND EITHER ARC OR
  16  TANGENT ) INTO A COMPLETE FAN-BEAM OF PROJECTION
  17  DATA
  18  NUPROJ(L),L=1,NURAYS,
  19  AT DETECTOR SPACING NUPINC ( AS IF WE HAD THE ARC
  20  CASE ), WHERE NURAYS IS EVEN, AND NURAYS*NUPINC
  21  EQUALS PI*STOD, SUCH THAT NUPROJ(L) IS THE LINE INTEGRAL
  22  FOR THE RAY THROUGH THE SOURCE, THAT MAKES ANGLE
  23  L*(NUPINC/RADIUS) RADIANS
  24  WITH THE LINE ORTHOGONAL TO THE CENTRAL RAY IN THE
  25  DIVERGENT BEAM.THE CENTRAL RAYS ARE PROJEC(NRAYS/2+1)
  26  AND PROJEC(NURAYS/2+1).
  27  NUPNTS HAS TO BE INITIALIZED AS FOLLOWS:
  28  NUPNTS(NURAYS/2-I)=I*NUPINC/PINC, I=1,...,NURAYS/2-1,
  29  SO NUPNTS(NURAYS/2-I) IS THE SCALED POSITION OF THE
  30  I-TH RAY ON THE CIRCLE OF DETECTORS, MEASURED RELATIVE
  31  TO THE CENTRAL RAY.
  32
  33  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  34  */
  35
  36 #include <cstdio>
  37
  38 #include "blkdta.h"
  39 #include "uiod.h"
  40
  41 #include "congeo.h"
  42
  43 void congeo(REAL* projec, INTEGER nrays, REAL* nuproj, INTEGER nurays, REAL* nupnts)
  44 {
  45
  46         INTEGER nur2;
  47         INTEGER nra2;
  48         INTEGER indist;
  49         INTEGER j;
  50         INTEGER indnu;
  51         INTEGER indpro;
  52         INTEGER i;
  53         REAL space;
  54         INTEGER irefl;
  55         INTEGER indref;
  56         INTEGER nr2;
  57         INTEGER l;
  58
  59         nur2 = nurays / 2;
  60         nra2 = nrays / 2;
  61
  62         // THE FIRST (TANGENT) RAY
  63
  64         nuproj[0] = 0.0;
  65
  66         // THE CENTER RAY
  67
  68         nuproj[nur2] = projec[nra2];
  69
  70         // FILL IN POINTS OUTSIDE ORIGINAL TABLE
  71
  72         indist = nra2;
  73         j = 1;
  74         for (;;)
  75         {
  76                 j++;
  77                 if (nupnts[j - 1] < float(indist))
  78                         break;
  79                 nuproj[j] = 0.0;
  80                 indnu = nurays - j + 2;
  81                 nuproj[indnu] = 0.0;
  82         }
  83
  84 //      THE POINTS INSIDE THE TABLE
  85
  86         indist = indist - 1;
  87         indpro = nra2 - indist + 1;
  88
  89         for (i = j; i < nur2; i++)
  90         {
  91                 for (;;)
  92                 {
  93                         space = nupnts[i - 1] - ((REAL) (indist));
  94                         if (space >= 0)
  95                                 break;
  96                         indist--;
  97                         indpro++;
  98                 }
  99                 nuproj[i] = space * projec[indpro - 1] + ((REAL) 1.0 - space) * projec[indpro];
 100
 101                 // THE CORRESPONDING POINTS ON THE NEGATIVE AXIS
 102
 103                 irefl = nurays - i + 2;
 104                 indref = nrays - indpro + 1;
 105                 nuproj[irefl] = ((REAL) 1.0 - space) * projec[indref]
 106                                 + space * projec[indref + 1];
 107
 108         }
 109
 110         if (trace < 5)
 111                 return;
 112
 113         nr2 = nurays / 2 - 1;
 114
 115         fprintf(output, "\n         interpolation performed at points:");
 116         for (l = 0; l < nr2; l++)
 117         {
 118                 if ((l % 5) == 0)
 119                 {
 120                         fprintf(output, "\n          ");
 121                 }
 122                 fprintf(output, "%13.6e", nupnts[l]);
 123         }
 124
 125         fprintf(output, "\n         projection data");
 126         for (l = 0; l < nrays; l++)
 127         {
 128                 if ((l % 5) == 0)
 129                 {
 130                         fprintf(output, "\n          ");
 131                 }
 132                 fprintf(output, "%13.6e", projec[l]);
 133         }
 134
 135         fprintf(output, "\n         interpolated projection data");
 136         for (l = 0; l < nurays; l++)
 137         {
 138                 if ((l % 5) == 0)
 139                 {
 140                         fprintf(output, "\n          ");
 141                 }
 142                 fprintf(output, " %13.6e", nuproj[l]);
 143         }
 144
 145         trace--;
 146
 147         return;
 148 }