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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <mach.h>
include "im2interpdef.h"
include <math/iminterp.h>
# MSISQGRL -- Procedure to integrate the 2D interpolant over a rectangular
# region.
real procedure msisqgrl (msi, x1, x2, y1, y2)
pointer msi # pointer to the interpolant descriptor structure
real x1, x2 # x integration limits
real y1, y2 # y integration limits
int i, interp_type, nylmin, nylmax, offset
pointer xintegrl, ptr
real xmin, xmax, ymin, ymax, accum
real ii_1dinteg()
begin
# set up 1D interpolant type to match 2-D interpolant
switch (MSI_TYPE(msi)) {
case II_BINEAREST:
interp_type = II_NEAREST
case II_BILINEAR:
interp_type = II_LINEAR
case II_BIDRIZZLE:
interp_type = II_DRIZZLE
case II_BIPOLY3:
interp_type = II_POLY3
case II_BIPOLY5:
interp_type = II_POLY5
case II_BISPLINE3:
interp_type = II_SPLINE3
case II_BISINC:
interp_type = II_SINC
case II_BILSINC:
interp_type = II_LSINC
}
# set up temporary storage for x integrals
call calloc (xintegrl, MSI_NYCOEFF(msi), TY_REAL)
# switch order of x integration at the end
xmin = x1
xmax = x2
if (x2 < x1) {
xmax = x2
xmin = x1
}
# switch order of y integration at end
ymin = y1
ymax = y2
if (y2 < y1) {
ymax = y2
ymin = y1
}
# find the appropriate range in y in the coeff array
offset = mod (MSI_FSTPNT(msi), MSI_NXCOEFF(msi))
nylmin = max (1, min (MSI_NYCOEFF(msi), int (ymin + 0.5)))
nylmax = min (MSI_NYCOEFF(msi), max (1, int (ymax + 0.5)))
nylmin = nylmin + offset
nylmax = nylmax + offset
# integrate in x
ptr = MSI_COEFF(msi) + offset + (nylmin - 1) * MSI_NXCOEFF(msi)
do i = nylmin, nylmax {
Memr[xintegrl+i-1] = ii_1dinteg (COEFF(ptr), MSI_NXCOEFF(msi),
xmin, xmax, interp_type, MSI_NSINC(msi), DX, MSI_XPIXFRAC(msi))
if (x2 < x1)
Memr[xintegrl+i-1] = - Memr[xintegrl+i-1]
ptr = ptr + MSI_NXCOEFF(msi)
}
# integrate in y
if (interp_type == II_SPLINE3) {
call amulkr (Memr[xintegrl], 6.0, Memr[xintegrl],
MSI_NYCOEFF(msi))
accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi),
ymin, ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
} else
accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi),
ymin, ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
# free space
call mfree (xintegrl, TY_REAL)
# correct for integration error.
if (y2 < y1)
return (-accum)
else
return (accum)
end
|
