# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	"../icombine.h"

define	MINCLIP		3	# Minimum number of images for this algorithm

$for (sird)
# IC_AAVSIGCLIP -- Reject pixels using an average sigma about the average
# The average sigma is normalized by the expected poisson sigma.

procedure ic_aavsigclip$t (d, m, n, scales, zeros, nimages, npts, average)

pointer	d[nimages]		# Data pointers
pointer	m[nimages]		# Image id pointers
int	n[npts]			# Number of good pixels
real	scales[nimages]		# Scales
real	zeros[nimages]		# Zeros
int	nimages			# Number of images
int	npts			# Number of output points per line
$if (datatype == sil)
real	average[npts]		# Average
$else
PIXEL	average[npts]		# Average
$endif

int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
$if (datatype == sil)
real	d1, low, high, sum, a, s, s1, r, one
data	one /1.0/
$else
PIXEL	d1, low, high, sum, a, s, s1, r, one
data	one /1$f/
$endif
pointer	sp, sums, resid, dp1, dp2, mp1, mp2

include	"../icombine.com"

begin
	# If there are insufficient pixels go on to the combining.
	if (nkeep < 0)
	    maxkeep = max (0, nimages + nkeep)
	else
	    maxkeep = min (nimages, nkeep)
	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
	    docombine = true
	    return
	}

	call smark (sp)
	call salloc (sums, npts, TY_REAL)
	call salloc (resid, nimages+1, TY_REAL)

	# Since the unweighted average is computed here possibly skip combining
	if (dowts || combine != AVERAGE)
	    docombine = true
	else
	    docombine = false

	# Compute the unweighted average with the high and low rejected and
	# the poisson scaled average sigma.  There must be at least three
	# pixels at each point to define the average and contributions to
	# the mean sigma.  Corrections for differences in the image
	# scale factors are selected by the doscale1 flag.

	nin = n[1]
	s = 0.
	n2 = 0
	do i = 1, npts {
	    k = i - 1
	    n1 = n[i]
	    if (n1 < 3)
		next

	    # Unweighted average with the high and low rejected
	    low = Mem$t[d[1]+k]
	    high = Mem$t[d[2]+k]
	    if (low > high) {
		d1 = low
		low = high
		high = d1
	    }
	    sum = 0.
	    do j = 3, n1 {
	       d1 = Mem$t[d[j]+k]
	       if (d1 < low) {
		   sum = sum + low
		   low = d1
		} else if (d1 > high) {
		    sum = sum + high
		    high = d1
		} else
		    sum = sum + d1
	    }
	    a = sum / (n1 - 2)
	    sum = sum + low + high

	    # Poisson scaled sigma accumulation
	    if (doscale1) {
		do j = 1, n1 {
		    dp1 = d[j] + k
		    mp1 = m[j] + k

		    d1 = Mem$t[dp1]
		    l = Memi[mp1]
		    s1 = max (one, (a + zeros[l]) /  scales[l])
		    s = s + (d1 - a) ** 2 / s1
		}
	    } else {
		s1 = max (one, a)
		do j = 1, n1
		    s = s + (Mem$t[d[j]+k] - a) ** 2 / s1
	    }
	    n2 = n2 + n1

	    # Save the average and sum for later.
	    average[i] = a
	    Memr[sums+k] = sum
	}

	# Here is the final sigma.
	if (n2 > 1)
	    s = sqrt (s / (n2 - 1))

	# Reject pixels and compute the final average (if needed).
	# There must be at least three pixels at each point for rejection.
	# Iteratively scale the mean sigma and reject pixels
	# Compact the data and keep track of the image IDs if needed.

	do i = 1, npts {
	    k = i - 1
	    n1 = n[i]
	    if (nkeep < 0)
		maxkeep = max (0, n1 + nkeep)
	    else
		maxkeep = min (n1, nkeep)
	    if (n1 <= max (2, maxkeep)) {
		if (!docombine) {
		    if (n1 == 0)
			average[i] = blank
		    else {
			sum = Mem$t[d[1]+k]
			do j = 2, n1
			    sum = sum + Mem$t[d[j]+k]
			average[i] = sum / n1
		    }
		}
		next
	    }

	    a = average[i]
	    sum = Memr[sums+k]

	    repeat {
		n2 = n1
		if (s > 0.) {
		    if (doscale1) {
			for (j=1; j<=n1; j=j+1) {
			    dp1 = d[j] + k
			    mp1 = m[j] + k

			    d1 = Mem$t[dp1]
			    l = Memi[mp1]
			    s1 = s * sqrt (max (one, (a+zeros[l]) /  scales[l]))
			    r = (d1 - a) / s1
			    if (r < -lsigma || r > hsigma) {
				Memr[resid+n1] = abs(r)
				if (j < n1) {
				    dp2 = d[n1] + k
				    Mem$t[dp1] = Mem$t[dp2]
				    Mem$t[dp2] = d1
				    mp2 = m[n1] + k
				    Memi[mp1] = Memi[mp2]
				    Memi[mp2] = l
				    j = j - 1
				}
				sum = sum - d1
				n1 = n1 - 1
			    }
			}
		    } else {
			s1 = s * sqrt (max (one, a))
			for (j=1; j<=n1; j=j+1) {
			    dp1 = d[j] + k
			    d1 = Mem$t[dp1]
			    r = (d1 - a) / s1
			    if (r < -lsigma || r > hsigma) {
				Memr[resid+n1] = abs(r)
				if (j < n1) {
				    dp2 = d[n1] + k
				    Mem$t[dp1] = Mem$t[dp2]
				    Mem$t[dp2] = d1
				    if (keepids) {
					mp1 = m[j] + k
					mp2 = m[n1] + k
					l = Memi[mp1]
					Memi[mp1] = Memi[mp2]
					Memi[mp2] = l
				    }
				    j = j - 1
				}
				sum = sum - d1
				n1 = n1 - 1
			    }
			}
		    }
		}
		if (n1 > 1)
		    a = sum / n1
	    } until (n1 == n2 || n1 <= max (2, maxkeep))

	    # If too many are rejected add some back in.
	    # Pixels with equal residuals are added together.
	    if (n1 < maxkeep) {
		nk = maxkeep
		if (doscale1) {
		    for (j=n1+1; j<=nk; j=j+1) {
			dp1 = d[j] + k
			mp1 = m[j] + k
			r = Memr[resid+j]
			jj = 0
			do l = j+1, n2 {
			    s = Memr[resid+l]
			    if (s < r + TOL) {
				if (s > r - TOL)
				    jj = jj + 1
				else {
				    jj = 0
				    Memr[resid+l] = r
				    r = s
				    dp2 = d[l] + k
				    d1 = Mem$t[dp1]
				    Mem$t[dp1] = Mem$t[dp2]
				    Mem$t[dp2] = d1
				    mp2 = m[l] + k
				    s = Memi[mp1]
				    Memi[mp1] = Memi[mp2]
				    Memi[mp2] = s
				}
			    }
			}
			sum = sum + Mem$t[dp1]
			n1 = n1 + 1
			nk = max (nk, j+jj)
		    }
		} else {
		    for (j=n1+1; j<=nk; j=j+1) {
			dp1 = d[j] + k
			r = Memr[resid+j]
			jj = 0
			do l = j+1, n2 {
			    s = Memr[resid+l]
			    if (s < r + TOL) {
				if (s > r - TOL)
				    jj = jj + 1
				else {
				    jj = 0
				    Memr[resid+l] = r
				    r = s
				    dp2 = d[l] + k
				    d1 = Mem$t[dp1]
				    Mem$t[dp1] = Mem$t[dp2]
				    Mem$t[dp2] = d1
				    if (keepids) {
					mp1 = m[j] + k
					mp2 = m[l] + k
					s = Memi[mp1]
					Memi[mp1] = Memi[mp2]
					Memi[mp2] = s
				    }
				}
			    }
			}
			sum = sum + Mem$t[dp1]
			n1 = n1 + 1
			nk = max (nk, j+jj)
		    }
		}
		if (n1 > 1)
		    a = sum / n1
	    }

	    # Save the average if needed.
	    n[i] = n1
	    if (!docombine) {
		if (n1 > 0)
		    average[i] = a
		else
		    average[i] = blank
	    }
	}

	# Check if the data flag has to be reset for rejected pixels
	if (dflag == D_ALL) {
	    do i = 1, npts {
		if (n[i] != nin) {
		    dflag = D_MIX
		    break
		}
	    }
	}

	call sfree (sp)
end


# IC_MAVSIGCLIP -- Reject pixels using an average sigma about the median
# The average sigma is normalized by the expected poisson sigma.

procedure ic_mavsigclip$t (d, m, n, scales, zeros, nimages, npts, median)

pointer	d[nimages]		# Data pointers
pointer	m[nimages]		# Image id pointers
int	n[npts]			# Number of good pixels
real	scales[nimages]		# Scales
real	zeros[nimages]		# Zeros
int	nimages			# Number of images
int	npts			# Number of output points per line
$if (datatype == sil)
real	median[npts]		# Median
$else
PIXEL	median[npts]		# Median
$endif

int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
pointer	sp, resid, mp1, mp2
$if (datatype == sil)
real	med, low, high, r, s, s1, one
data	one /1.0/
$else
PIXEL	med, low, high, r, s, s1, one
data	one /1$f/
$endif

include	"../icombine.com"

begin
	# If there are insufficient pixels go on to the combining.
	if (nkeep < 0)
	    maxkeep = max (0, nimages + nkeep)
	else
	    maxkeep = min (nimages, nkeep)
	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
	    docombine = true
	    return
	}

	call smark (sp)
	call salloc (resid, nimages+1, TY_REAL)

	# Compute the poisson scaled average sigma about the median.
	# There must be at least three pixels at each point to define
	# the mean sigma.  Corrections for differences in the image
	# scale factors are selected by the doscale1 flag.

	s = 0.
	n2 = 0
	nin = n[1]
	do i = 1, npts {
	    k = i - 1
	    n1 = n[i]
	    if (n1 < 3) {
		if (n1 == 0)
		    median[i] = blank
		else if (n1 == 1)
		    median[i] = Mem$t[d[1]+k]
		else {
		    low = Mem$t[d[1]+k]
		    high = Mem$t[d[2]+k]
		    median[i] = (low + high) / 2.
		}
		next
	    }

	    # Median
	    n3 = 1 + n1 / 2
	    if (mod (n1, 2) == 0) {
		low = Mem$t[d[n3-1]+k]
		high = Mem$t[d[n3]+k]
		med = (low + high) / 2.
	    } else
		med = Mem$t[d[n3]+k]

	    # Poisson scaled sigma accumulation
	    if (doscale1) {
		do j = 1, n1 {
		    l = Memi[m[j]+k]
		    s1 = max (one, (med + zeros[l]) /  scales[l])
		    s = s + (Mem$t[d[j]+k] - med) ** 2 / s1
		}
	    } else {
		s1 = max (one, med)
		do j = 1, n1
		    s = s + (Mem$t[d[j]+k] - med) ** 2 / s1
	    }
	    n2 = n2 + n1

	    # Save the median for later.
	    median[i] = med
	}

	# Here is the final sigma.
	if (n2 > 1)
	    s = sqrt (s / (n2 - 1))
	else {
	    call sfree (sp)
	    return
	}

	# Compute individual sigmas and iteratively clip.
	do i = 1, npts {
	    k = i - 1
	    n1 = n[i]
	    if (nkeep < 0)
		maxkeep = max (0, n1 + nkeep)
	    else
		maxkeep = min (n1, nkeep)
	    if (n1 < max (3, maxkeep+1))
		next
	    nl = 1
	    nh = n1
	    med = median[i]

	    repeat {
		n2 = n1
		n3 = nl + n1 / 2

		if (n1 >= max (MINCLIP, maxkeep+1) && s > 0.) {
		    if (doscale1) {
			for (; nl <= n2; nl = nl + 1) {
			    l = Memi[m[nl]+k]
			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
			    r = (med - Mem$t[d[nl]+k]) / s1
			    if (r <= lsigma)
				break
			    Memr[resid+nl] = r
			    n1 = n1 - 1
			}
			for (; nh >= nl; nh = nh - 1) {
			    l = Memi[m[nh]+k]
			    s1 = s * sqrt (max (one, (med+zeros[l])/scales[l]))
			    r = (Mem$t[d[nh]+k] - med) / s1
			    if (r <= hsigma)
				break
			    Memr[resid+nh] = r
			    n1 = n1 - 1
			}
		    } else {
			s1 = s * sqrt (max (one, med))
			for (; nl <= n2; nl = nl + 1) {
			    r = (med - Mem$t[d[nl]+k]) / s1
			    if (r <= lsigma)
				break
			    Memr[resid+nl] = r
			    n1 = n1 - 1
			}
			for (; nh >= nl; nh = nh - 1) {
			    r = (Mem$t[d[nh]+k] - med) / s1
			    if (r <= hsigma)
				break
			    Memr[resid+nh] = r
			    n1 = n1 - 1
			}
		    }

		    # Recompute median
		    if (n1 < n2) {
			if (n1 > 0) {
			    n3 = nl + n1 / 2
			    if (mod (n1, 2) == 0) {
				low = Mem$t[d[n3-1]+k]
				high = Mem$t[d[n3]+k]
				med = (low + high) / 2.
			    } else
				med = Mem$t[d[n3]+k]
			} else
			    med = blank
		    }
		}
	    } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1))

	    # If too many are rejected add some back in.
	    # Pixels with equal residuals are added together.
	    while (n1 < maxkeep) {
		if (nl == 1)
		    nh = nh + 1
		else if (nh == n[i])
		    nl = nl - 1
		else {
		    r = Memr[resid+nl-1]
		    s = Memr[resid+nh+1]
		    if (r < s) {
			nl = nl - 1
			r = r + TOL
			if (s <= r)
			    nh = nh + 1
			if (nl > 1) {
			    if (Memr[resid+nl-1] <= r)
				nl = nl - 1
			}
		    } else {
			nh = nh + 1
			s = s + TOL
			if (r <= s)
			    nl = nl - 1
			if (nh < n2) {
			    if (Memr[resid+nh+1] <= s)
				nh = nh + 1
			}
		    }
		}
		n1 = nh - nl + 1

		# Recompute median
		if (n1 < n2) {
		    if (n1 > 0) {
			n3 = nl + n1 / 2
			if (mod (n1, 2) == 0) {
			    low = Mem$t[d[n3-1]+k]
			    high = Mem$t[d[n3]+k]
			    med = (low + high) / 2.
			} else
			    med = Mem$t[d[n3]+k]
		    } else
			med = blank
		}
	    }

	    # Only set median and reorder if needed
	    n[i] = n1
	    if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow >= 1.)) {
		j = max (nl, n1 + 1)
		if (keepids) {
		    do l = 1, min (n1, nl-1) {
			Mem$t[d[l]+k] = Mem$t[d[j]+k]
			if (grow >= 1.) {
			    mp1 = m[l] + k
			    mp2 = m[j] + k
			    id = Memi[mp1]
			    Memi[mp1] = Memi[mp2]
			    Memi[mp2] = id
			} else
			    Memi[m[l]+k] = Memi[m[j]+k]
			j = j + 1
		    }
		} else {
		    do l = 1, min (n1, nl - 1) {
			Mem$t[d[l]+k] = Mem$t[d[j]+k]
			j = j + 1
		    }
		}
	    }

	    if (combine == MEDIAN)
		median[i] = med
	}

	# Check if data flag needs to be reset for rejected pixels
	if (dflag == D_ALL) {
	    do i = 1, npts {
		if (n[i] != nin) {
		    dflag = D_MIX
		    break
		}
	    }
	}

	# Flag that the median is computed.
	if (combine == MEDIAN)
	    docombine = false
	else
	    docombine = true

	call sfree (sp)
end
$endfor