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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
include <imhdr.h>
include <imio.h>
include <mach.h>
include "../lib/apphot.h"
include "../lib/noise.h"
include "../lib/find.h"
# AP_BFDFIND -- Find stars in an image using a pattern matching
# technique.
procedure ap_bfdfind (im, cnv, sky, out, ap, boundary, constant, verbose)
pointer im # pointer to the input image
pointer cnv # pointer to the convolved image
pointer sky # pointer to the sky image
int out # the output file descriptor
pointer ap # pointer to the apphot structure
int boundary # type of boundary extension
real constant # constant for constant boundary extension
int verbose # verbose switch
int norm, nxk, nyk, nstars, stid
pointer sp, gker2d, ngker2d, dker2d, skip
real a, b, c, f, gsums[LEN_GAUSS], skysigma, skymode, threshold, relerr
real dmax, dmin, xsigsq, ysigsq
int apstati(), ap_find()
real ap_egkernel(), apstatr()
begin
# Compute the parameters of the Gaussian kernel.
call ap_egparams (FWHM_TO_SIGMA * apstatr (ap, FWHMPSF) *
apstatr (ap, SCALE), apstatr (ap, RATIO), apstatr (ap, THETA),
apstatr (ap, NSIGMA), a, b, c, f, nxk, nyk)
# Allocate working space.
call smark (sp)
call salloc (gker2d, nxk * nyk, TY_REAL)
call salloc (ngker2d, nxk * nyk, TY_REAL)
call salloc (dker2d, nxk * nyk, TY_REAL)
call salloc (skip, nxk * nyk, TY_INT)
# Compute the 1 and 2 D kernels.
if (IS_INDEFR(apstatr(ap, DATAMIN)) && IS_INDEFR(apstatr(ap,
DATAMAX))) {
norm = YES
dmin = -MAX_REAL
dmax = MAX_REAL
} else {
norm = NO
if (IS_INDEFR(apstatr (ap, DATAMIN)))
dmin = -MAX_REAL
else
dmin = apstatr (ap, DATAMIN)
if (IS_INDEFR(apstatr (ap, DATAMAX)))
dmax = MAX_REAL
else
dmax = apstatr (ap, DATAMAX)
}
relerr = ap_egkernel (Memr[gker2d], Memr[ngker2d], Memr[dker2d],
Memi[skip], nxk, nyk, gsums, a, b, c, f)
# Set up the image boundary extension characteristics.
call ap_imset (im, boundary, max (1 + nxk / 2, 1 + nyk / 2),
constant)
call ap_imset (cnv, boundary, max (1 + nxk / 2, 1 + nyk / 2),
constant)
# Convolve the input image with the Gaussian kernel. The resultant
# picture constains in each pixel the height of the Gaussian
# function centered in the subarray which best represents the data
# within a circle of nsigma * sigma of the Gaussian.
if (IM_ACMODE(cnv) != READ_ONLY) {
if (norm == YES)
call ap_fconvolve (im, cnv, sky, Memr[ngker2d], Memr[dker2d],
Memi[skip], nxk, nyk, gsums[GAUSS_SGOP])
else
call ap_gconvolve (im, cnv, sky, Memr[gker2d], Memi[skip],
nxk, nyk, gsums, dmin, dmax)
}
# Save the task parameters in the database file if the savepars
# switch is enabled, otherwise a simple list of detected objects
# is written to the data base file.
call ap_wfdparam (out, ap)
# Find all the objects in the input image with the specified image
# characteristics.
if (verbose == YES) {
call printf ("\nImage: %s ")
call pargstr (IM_HDRFILE(im))
call printf ("fwhmpsf: %g ratio: %g theta: %g nsigma: %g\n\n")
call pargr (apstatr (ap, FWHMPSF))
call pargr (apstatr (ap, RATIO))
call pargr (apstatr (ap, THETA))
call pargr (apstatr (ap, NSIGMA))
}
# Get the skymode and threshold.
skysigma = apstatr (ap, SKYSIGMA)
if (IS_INDEFR(skysigma)) {
skymode = 0.0
threshold = 0.0
} else {
skymode = apstatr (ap, EPADU) * (skysigma ** 2 -
(apstatr (ap, READNOISE) / apstatr (ap, EPADU)) ** 2)
skymode = max (0.0, skymode)
threshold = apstatr (ap, THRESHOLD) * skysigma
}
# Compute the x and y sigma.
xsigsq = (apstatr (ap, SCALE) * apstatr (ap, FWHMPSF) / 2.35482) ** 2
ysigsq = (apstatr (ap, SCALE) * apstatr (ap, RATIO) *
apstatr (ap, FWHMPSF) / 2.35482) ** 2
# Find the stars.
stid = 1
nstars = ap_find (ap, im, cnv, out, NULL, Memr[gker2d],
Memi[skip], nxk, nyk, skymode, threshold, relerr,
apstati (ap,POSITIVE), xsigsq, ysigsq, dmin, dmax,
apstatr (ap, SHARPLO), apstatr (ap, SHARPHI), apstatr (ap,
ROUNDLO), apstatr (ap, ROUNDHI), verbose, stid, NO)
if (verbose == YES) {
call printf (
"\nthreshold: %g relerr: %5.3f %g <= sharp <= %g ")
call pargr (threshold)
call pargr (relerr)
call pargr (apstatr (ap, SHARPLO))
call pargr (apstatr (ap, SHARPHI))
call printf ("%g <= round <= %g \n\n")
call pargr (apstatr (ap, ROUNDLO))
call pargr (apstatr (ap, ROUNDHI))
}
call sfree (sp)
end
|
