1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
/* hst_calculate_limits.F -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/*#include "f2c.h"
*/
#define integer int
#define real float
#define max(a,b) ((a)>(b))? (a):(b)
#define min(a,b) ((a)<(b))? (a):(b)
/* Subroutine */ int hst_calculate_limits__(num_bins__, start_x__, start_y__,
num_x__, num_y__, axis_position__, y_start__, y_end__, x_starts__,
x_ends__, status)
integer *num_bins__, *start_x__, *start_y__, *num_x__, *num_y__;
real *axis_position__;
integer *y_start__, *y_end__, *x_starts__, *x_ends__, *status;
{
/* System generated locals */
integer i__1, i__2, i__3;
integer trim;
/* Builtin functions */
double sqrt();
/* Local variables */
static real minx, maxx;
static integer y;
static real x_diff__, radiusleft, radiussmall, radiusright;
/* Number of bins in each sinogram */
/* First X-pixel in region required by user */
/* First Y-pixel in region required by user */
/* Number of X-pixels in reconstruction region */
/* Number of Y-pixels in reconstruction region */
/* Position of rotation axis */
/* First Y-row in region to reconstruct */
/* Last Y-row in region to reconstruct */
/* First pixel to */
/* Last pixel to */
/* Loop variable for Y-direction */
/* X-difference for points on circle */
/* Parameter adjustments */
--x_ends__;
--x_starts__;
/* Function Body */
trim=1;
radiusleft = *axis_position__;
radiusright = (real) (*num_bins__) - *axis_position__;
if (radiusleft < radiusright) {
radiussmall = radiusleft;
} else {
radiussmall = radiusright;
}
/* Calculate start and end pixels of projection */
*y_start__ = *axis_position__ - radiusright + 1;
*y_end__ = *axis_position__ + radiusleft;
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* we dont go outside users rectangle */
if (*y_start__ < *start_y__) {
*y_start__ = *start_y__;
}
if (*y_end__ > *start_y__ + *num_y__ -1) {
*y_end__ = *start_y__ + *num_y__ -1;
}
/* !!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
maxx = *axis_position__ + radiussmall;
if (maxx > (real) (*num_bins__)) {
maxx = (real) (*num_bins__);
}
minx = *axis_position__ - radiussmall + 1;
if (minx < (float)1.) {
minx = (float)1.;
}
i__1 = *y_end__;
for (y = *y_start__; y <= i__1; ++y) {
if ((real) y < *axis_position__) {
x_diff__ = sqrt(radiusright * radiusright - (y - *axis_position__)
* (y - *axis_position__));
} else {
x_diff__ = sqrt(radiusleft * radiusleft - (y - *axis_position__) *
(y - *axis_position__));
}
if (x_diff__ > radiussmall) {
x_diff__ = radiussmall;
}
/* Computing MAX */
i__2 = 1, i__3 = (integer) (*axis_position__ - x_diff__) - *start_x__
+ 1;
x_starts__[y - *start_y__ + 1] = max(i__2,i__3)+trim;
/* Computing MIN */
i__2 = *num_x__, i__3 = (integer) (*axis_position__ + x_diff__) - *
start_x__ + 1;
x_ends__[y - *start_y__ + 1] = min(i__2,i__3) - trim;
}
*y_start__ = *y_start__ - *start_y__ + 1;
*y_end__ = *y_end__ - *start_y__ + 1;
} /* hst_calculate_limits__ */
|