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) $
8 **********************************************************************
10 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
12 THIS SUBROUTINE TRANSFORMS ONE DIVERGENT BEAM OF
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
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
33 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
43 void congeo(REAL* projec, INTEGER nrays, REAL* nuproj, INTEGER nurays, REAL* nupnts)
62 // THE FIRST (TANGENT) RAY
68 nuproj[nur2] = projec[nra2];
70 // FILL IN POINTS OUTSIDE ORIGINAL TABLE
77 if (nupnts[j - 1] < float(indist))
80 indnu = nurays - j + 2;
84 // THE POINTS INSIDE THE TABLE
87 indpro = nra2 - indist + 1;
89 for (i = j; i < nur2; i++)
93 space = nupnts[i - 1] - ((REAL) (indist));
99 nuproj[i] = space * projec[indpro - 1] + ((REAL) 1.0 - space) * projec[indpro];
101 // THE CORRESPONDING POINTS ON THE NEGATIVE AXIS
103 irefl = nurays - i + 2;
104 indref = nrays - indpro + 1;
105 nuproj[irefl] = ((REAL) 1.0 - space) * projec[indref]
106 + space * projec[indref + 1];
113 nr2 = nurays / 2 - 1;
115 fprintf(output, "\n interpolation performed at points:");
116 for (l = 0; l < nr2; l++)
120 fprintf(output, "\n ");
122 fprintf(output, "%13.6e", nupnts[l]);
125 fprintf(output, "\n projection data");
126 for (l = 0; l < nrays; l++)
130 fprintf(output, "\n ");
132 fprintf(output, "%13.6e", projec[l]);
135 fprintf(output, "\n interpolated projection data");
136 for (l = 0; l < nurays; l++)
140 fprintf(output, "\n ");
142 fprintf(output, " %13.6e", nuproj[l]);