aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/ccdred/src/sigma.gx
blob: 8b59f1f69a67ed97a0d126443edaeded96b6d757 (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
$for (sr)
# SIGMA -- Compute sigma line from image lines with rejection.

procedure sigma$t (data, nimages, mean, sigma, npts)

pointer	data[nimages]		# Data vectors
int	nimages			# Number of data vectors
$if (datatype == sil)
real	mean[npts]		# Mean vector
real	sigma[npts]		# Sigma vector (returned)
$else
PIXEL	mean[npts]		# Mean vector
PIXEL	sigma[npts]		# Sigma vector (returned)
$endif
int	npts			# Number of points in each vector

$if (datatype == sil)
real	val, sig, pixval
$else
PIXEL	val, sig, pixval
$endif
int	i, j, n, n1

begin
	n = nimages - 1
	do i = 1, npts {
	    val = mean[i]
	    sig = 0.
	    n1 = n
	    do j = 1, nimages {
		pixval = Mem$t[data[j]+i-1]
		if (IS_INDEF (pixval))
		    n1 = n1 - 1
		else
		    sig = sig + (pixval - val) ** 2
	    }
	    if (n1 > 0)
	        sigma[i] = sqrt (sig / n1)
	    else
	        sigma[i] = 0.
	}
end


# WTSIGMA -- Compute scaled and weighted sigma line from image lines with
# rejection.

procedure wtsigma$t (data, scales, zeros, wts, nimages, mean, sigma, npts)

pointer	data[nimages]		# Data vectors
real	scales[nimages]		# Scale factors
real	zeros[nimages]		# Zero levels
real	wts[nimages]		# Weights
int	nimages			# Number of data vectors
$if (datatype == sil)
real	mean[npts]		# Mean vector
real	sigma[npts]		# Sigma vector (returned)
real	val, sig, pixval
$else
PIXEL	mean[npts]		# Mean vector
PIXEL	sigma[npts]		# Sigma vector (returned)
PIXEL	val, sig, pixval
$endif
int	npts			# Number of points in each vector

int	i, j, n
real	sumwts

begin
	do i = 1, npts {
	    val = mean[i]
	    n = 0
	    sig = 0.
	    sumwts = 0.
	    do j = 1, nimages {
		pixval = Mem$t[data[j]+i-1]
		if (!IS_INDEF (pixval)) {
		    n = n + 1
		    sig = sig + wts[j]*(pixval/scales[j]-zeros[j]-val) ** 2
		    sumwts = sumwts + wts[j]
		}
	    }
	    if (n > 1)
	        sigma[i] = sqrt (sig / sumwts * n / (n - 1))
	    else
	        sigma[i] = 0.
	}
end
$endfor