aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/daophot/psf/dppsubrast.x
blob: 1fafb408547d6aaef445155d60f62394d6a219cc (plain) (blame)
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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
include <mach.h>
include <imhdr.h>
include "../lib/daophotdef.h"
include "../lib/psfdef.h"

# DP_PSUBRAST -- Fetch the prospective PSF star data and check it for bad
# pixel values.

pointer procedure dp_psubrast (dao, im, lowbad, highbad, x1, x2, y1, y2,
	saturated)

pointer	dao		# pointer to the daophot structure
pointer	im		# pointer to the input image
real	lowbad, highbad	# minimum and maximum good data values
int	x1, x2		# output x limits of the extracted subraster
int	y1, y2		# output y limits of the extracted subraster
int	saturated	# is the star saturated ?

pointer	psf, buf
real	psfrad, fitrad
int	dp_chksr()
pointer	dp_subrast()

begin
	# Initialize.
	psf = DP_PSF(dao)
	buf = NULL
	psfrad = DP_PSFRAD(dao)
	fitrad = DP_FITRAD(dao)

	# Get the data.
 	buf = dp_subrast (im, DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), psfrad,
	    x1, x2, y1, y2)
	if (buf == NULL)
	    return (NULL)

	# Check for bad pixels in subraster, compute the min and max.
	if (dp_chksr (DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), Memr[buf],
	    x2 - x1 + 1, y2 - y1 + 1, x1, y1, psfrad, fitrad,
	    lowbad, highbad, saturated, DP_CUR_PSFMIN(psf),
	        DP_CUR_PSFMAX(psf), DP_CUR_PSFGMAX(psf)) == ERR) {
	    call mfree (buf, TY_REAL)
	    return (NULL)
	}
	
	return (buf)
end


# DL_LSUBRAST -- Give a valid PSF star compute the limits of the data to
# be extracted around it.

int procedure dp_lsubrast (im, xcen, ycen, radius, x1, x2, y1, y2)

pointer	im			# input image descriptor
real	xcen, ycen		# center of subraster
real	radius			# radius of the box
int	x1, y1, x2, y2		# boundaries of subraster

begin
	# Calculate start position of extraction box.
	x1 = int (xcen - radius) - 2
	x2 = int (xcen + radius) + 3
	y1 = int (ycen - radius) - 2
	y2 = int (ycen + radius) + 3
	if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1)
	    return (ERR)

	x1 = max (1, x1)
	x2 = min (IM_LEN(im,1), x2)
	y1 = max (1, y1)
	y2 = min (IM_LEN(im,2), y2)
	return (OK)
end


# DP_SUBRAST -- Given a valid PSF star extract the data around it.

pointer procedure dp_subrast (im, xcen, ycen, radius, x1, x2, y1, y2)

pointer	im			# input image descriptor
real	xcen, ycen		# center of subraster
real	radius			# radius of the box
int	x1, y1, x2, y2		# boundaries of subraster

int	j, ncols
pointer	buf, ptr, imbuf
pointer imgs2r()

begin
	# Calculate start position of extraction box.
	x1 = int (xcen - radius) - 2
	x2 = int (xcen + radius) + 3
	y1 = int (ycen - radius) - 2
	y2 = int (ycen + radius) + 3
	if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1)
	    return (NULL)

	x1 = max (1, x1)
	x2 = min (IM_LEN(im,1), x2)
	y1 = max (1, y1)
	y2 = min (IM_LEN(im,2), y2)
	call malloc (buf, (x2 - x1 + 1) * (y2 - y1 + 1), TY_REAL)

	ptr = buf
	ncols = x2 - x1 + 1
	do j = y1, y2 {
	    imbuf = imgs2r (im, x1, x2, j, j)
	    call amovr (Memr[imbuf], Memr[ptr], ncols) 
	    ptr = ptr + ncols
	}

	return (buf)
end


# DP_CHKSR -- Check the input subraster for bad pixels.

int procedure dp_chksr (x, y, sr, xdim, ydim, x1, y1, psfrad, fitrad, lowbad,
        highbad, saturated, dmin, dmax, gmax)

real	x, y			# position of the star
real	sr[xdim,ydim]		# the data subraster
int	xdim, ydim		# the dimensions of the subraster
int	x1, y1			# the lower left hand coordinates of the array
real	psfrad			# the psf radius
real	fitrad			# the fitting radius
real	lowbad, highbad		# the good data limits
int	saturated		# is the star saturated
real	dmin, dmax		# output data limits
real	gmax			# maximum good data limit

int	i,j
real	pradsq, fradsq, dy2, r2

begin
	pradsq = psfrad * psfrad
	fradsq = fitrad * fitrad
	dmin = MAX_REAL
	dmax = -MAX_REAL
	gmax = -MAX_REAL
	saturated = NO

	# Loop through the pixels looking for bad values.
	do j = 1, ydim {
	    dy2 = (y - (y1 + j -1)) ** 2 
	    if (dy2 > pradsq)
		next
	    do i = 1, xdim {
		r2 = (x - (x1 + i - 1)) ** 2 + dy2
	        if (r2 > pradsq)
		    next
		if (sr[i,j] < dmin)
		    dmin = sr[i,j]
		if (sr[i,j] > dmax)
		    dmax = sr[i,j]
		if (sr[i,j] < lowbad || sr[i,j] > highbad) {
		    if (r2 <= fradsq) {
		        if (sr[i,j] < lowbad)
		    	    return (ERR)
			else if (saturated == NO)
			    saturated = YES
		    }
		} else {
		    if (sr[i,j] > gmax)
		        gmax = sr[i,j]
		}
	    }
	}

	return (OK)
end