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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <mach.h>
# ARAV -- Compute the mean and standard deviation of a sample array by
# iteratively rejecting points further than KSIG from the mean. If the
# value of KSIG is given as 0.0, a cutoff value will be automatically
# calculated from the standard deviation and number of points in the sample.
# The number of pixels remaining in the sample upon termination is returned
# as the function value.
int procedure arav$t (a, npix, mean, sigma, ksig)
PIXEL a[ARB] # input data array
$if (datatype == dl)
double mean, sigma, ksig, deviation, lcut, hcut, lgpx
$else
real mean, sigma, ksig, deviation, lcut, hcut, lgpx
$endif
int npix, ngpix, old_ngpix, awvg$t()
begin
lcut = -MAX_REAL # no rejection to start
hcut = MAX_REAL
ngpix = MAX_INT
# Iteratively compute mean, sigma and reject outliers until no
# more pixels are rejected, or until there are no more pixels.
repeat {
old_ngpix = ngpix
ngpix = awvg$t (a, npix, mean, sigma, lcut, hcut)
$if (datatype == dl)
if (ngpix <= 1 || sigma <= EPSILOND)
$else
if (ngpix <= 1 || sigma <= EPSILONR)
$endif
break
if (ksig == 0.0) { # Chauvenet's relation
lgpx = log10 (real(ngpix))
deviation = (lgpx * (-0.1042 * lgpx + 1.1695) + .8895) * sigma
} else
deviation = sigma * abs(ksig)
lcut = mean - deviation # compute window
hcut = mean + deviation
} until (ngpix >= old_ngpix)
return (ngpix)
end
|
