# PEAKS -- The following procedures are general numerical functions
# dealing with finding peaks in a data array.
#
# FIND_PEAKS		Find the NMAX peaks in the data array.
# FIND_UPEAKS		Find the uniformly distrib. peaks in the data array.
# FIND_IPEAKS		Find all the isolated peaks in the data array.
# FIND_LOCAL_MAXIMA	Find the local maxima in the data array.
# IS_LOCAL_MAX		Test a point to determine if it is a local maximum.
# FIND_THRESHOLD	Find the peaks with positions satisfying threshold
#			and contrast constraints.
# FIND_ISOLATED		Flag peaks which are within separation of a peak
#			with a higher peak value.
# FIND_NMAX		Select up to the nmax highest ranked peaks.
# FIND_UNMAX		Select up to the nmax ranked peaks in bins.
# COMPARE		Compare procedure for sort used in FIND_PEAKS.


# FIND_PEAKS -- Find the NMAX peaks in the data array.
#
# The peaks are found using the following algorithm:
#
# 1.  Find the local maxima.
# 2.  Reject peaks below the threshold.
# 3.  Determine the ranks of the remaining peaks.
# 4.  Flag weaker peaks within separation of a stronger peak.
# 5.  Accept at most the nmax strongest peaks.
#
# Indefinite points are ignored.  The peak positions are returned in the
# array x.


int procedure find_peaks (data, x, npoints, contrast, separation, edge, nmax,
    threshold, debug)

# Procedure parameters:
real	data[npoints]		# Input data array
real	x[npoints]		# Output peak position array
int	npoints			# Number of data points
real	contrast		# Maximum contrast between strongest and weakest
int	separation		# Minimum separation between peaks
int	edge			# Minimum distance from the edge
int	nmax			# Maximum number of peaks to be returned
real	threshold		# Minimum threshold level for peaks
bool	debug			# Print diagnostic information?

int	nrank, npeaks, find_nmax()
pointer	rank

begin
	# Find all isolated peaks and their rank.
	call find_ipeaks (data, x, npoints, contrast, separation, edge,
    	    threshold, rank, nrank, debug)

	# Select the strongest nmax peaks.
	npeaks = find_nmax (data, x, Memi[rank], nrank, nmax, debug)

	call mfree (rank, TY_INT)
	return (npeaks)
end


# FIND_UPEAKS -- Find the uniformly distrib. peaks in the data array.
#
# The peaks are found using the following algorithm:
#
# 1.  Find the local maxima.
# 2.  Reject peaks below the threshold.
# 3.  Determine the ranks of the remaining peaks.
# 4.  Flag weaker peaks within separation of a stronger peak.
# 5.  Accept at most the nmax uniformly distributed peaks.
#
# Indefinite points are ignored.  The peak positions are returned in the
# array x.

int procedure find_upeaks (data, x, npoints, contrast, separation, edge,
	nmax, nbins, threshold, debug)

# Procedure parameters:
real	data[npoints]		# Input data array
real	x[npoints]		# Output peak position array
int	npoints			# Number of data points
real	contrast		# Maximum contrast between strongest and weakest
int	separation		# Minimum separation between peaks
int	edge			# Minimum distance from the edge
int	nmax			# Maximum number of peaks to be returned
int	nbins			# Number of bins across the data array
real	threshold		# Minimum threshold level for peaks
bool	debug			# Print diagnostic information?

int	npts, nrank, npeaks, find_unmax()
pointer	rank

begin
	npts = npoints

	# Find all isolated peaks and their rank.
	call find_ipeaks (data, x, npoints, contrast, separation, edge,
    	    threshold, rank, nrank, debug)

	# Select the peaks.
	npeaks = find_unmax (data, npts, x, Memi[rank], nrank, nmax, nbins,
	    debug)

	call mfree (rank, TY_INT)
	return (npeaks)
end


# FIND_IPEAKS -- Find the all the isolated peaks in the data array.
#
# The peaks are found using the following algorithm:
#
# 1.  Find the local maxima.
# 2.  Reject peaks below the threshold.
# 3.  Determine the ranks of the remaining peaks.
# 4.  Flag weaker peaks within separation of a stronger peak.
# 5.  Return a rank array
#
# Indefinite points are ignored.  The peak positions are returned in the
# array x.

procedure find_ipeaks (data, x, npoints, contrast, separation, edge, threshold,
	rank, nrank, debug)

# Procedure parameters:
real	data[npoints]		# Input data array
real	x[npoints]		# Output peak position array
int	npoints			# Number of data points
real	contrast		# Maximum contrast between strongest and weakest
int	separation		# Minimum separation between peaks
int	edge			# Minimum distance from the edge
real	threshold		# Minimum threshold level for peaks
pointer	rank			# Rank array
int	nrank			# Size of rank array
bool	debug			# Print diagnostic information?

int	i, j
int	nlmax, nisolated
pointer	sp, y

int	find_local_maxima(), find_threshold(), find_isolated()
int	compare()

extern	compare()

common	/sort/ y

begin
	# Find the local maxima in data and put column positions in x..
	nlmax = find_local_maxima (data, x, npoints, debug)

	# Reject local maxima near the edge.
	if (edge > 0) {
	    j = 0
	    do i = 1, nlmax {
		if ((x[i] > edge) && (x[i] <= npoints - edge)) {
		    j = j + 1
		    x[j] = x[i]
		}
	    }
	    nlmax = j
	}

	# Allocate a working array y.
	call smark (sp)
	call salloc (y, npoints, TY_REAL)

	# Reject the local maxima which do not satisfy the thresholds.
	# The array y is set to the peak values of the remaining peaks.
	nrank = find_threshold (data, x, Memr[y], nlmax,
	    contrast, threshold, debug)

	# Rank the peaks by peak value.
	call malloc (rank, nrank, TY_INT)
	do i = 1, nrank
	    Memi[rank + i - 1] = i
	call qsort (Memi[rank], nrank, compare)

	# Reject the weaker peaks within sep of a stronger peak.
	nisolated = find_isolated (x, Memi[rank], nrank, separation, debug)

	call sfree (sp)
end


# FIND_LOCAL_MAXIMA -- Find the local maxima in the data array.
#
# A data array is input and the local maxima positions array is output.
# The number of local maxima found is returned.

int procedure find_local_maxima (data, x, npoints, debug)

real	data[npoints]			# Input data array
real	x[npoints]			# Output local maxima positions array
int	npoints				# Number of input points
bool	debug				# Print debugging information?

int	i, nlmax

bool	is_local_max()

begin
	nlmax = 0
	do i = 1, npoints {
	    if (is_local_max (i, data, npoints)) {
		nlmax = nlmax + 1
		x[nlmax] = i
	    }
	}

	if (debug) {
	    call printf ("  Number of local maxima found = %d.\n")
		call pargi (nlmax)
	}

	return (nlmax)
end


# IS_LOCAL_MAX -- Test a point to determine if it is a local maximum.
#
# Indefinite points are ignored.

bool procedure is_local_max (index, data, npoints)

# Procedure parameters:
int	index			# Index to test for local maximum
real	data[npoints]		# Data values
int	npoints			# Number of points in the data vector

int	i, j, nright, nleft

begin
	# INDEF points cannot be local maxima.
	if (IS_INDEFR (data[index]))
	    return (FALSE)

	# Find the left and right indices where data values change and the
	# number of points with the same value.  Ignore INDEF points.
	nleft = 0
	for (i = index - 1; i >= 1; i = i - 1) {
	    if (!IS_INDEFR (data[i])) {
		if (data[i] != data[index])
		    break
		nleft = nleft + 1
	    }
	}
	nright = 0
	for (j = index + 1; j <= npoints; j = j + 1) {
	    if (!IS_INDEFR (data[j])) {
		if (data[j] != data[index])
		    break
		nright = nright + 1
	    }
	}

	# Test for failure to be a local maxima
	if ((i == 0) && (j == npoints+1)) {
	    return (FALSE)		# Data is constant
	} else if (i == 0) {
	    if (data[j] > data[index])
	        return (FALSE)		# Data increases to right
	} else if (j == npoints+1) {
	    if (data[i] > data[index])	# Data increase to left
		return (FALSE)
	} else if ((data[i] > data[index]) || (data[j] > data[index])) {
	    return (FALSE)		# Not a local maximum
	} else if (!((nleft - nright == 0) || (nleft - nright == 1))) {
	    return (FALSE)		# Not center of plateau
	}

	# Point is a local maxima
	return (TRUE)
end


# FIND_THRESHOLD -- Find the peaks with positions satisfying threshold
# and contrast constraints.
#
# The input is the data array, data, and the peak positions array, x.
# The x array is resorted to the nthreshold peaks satisfying the constraints.
# The corresponding nthreshold data values are returned the y array.
# The number of peaks satisfying the constraints (nthreshold) is returned.

int procedure find_threshold (data, x, y, npoints, contrast, threshold, debug)

real	data[ARB]			# Input data values
real	x[npoints]			# Input/Output peak positions
real	y[npoints]			# Output peak data values
int	npoints				# Number of peaks input
real	contrast			# Contrast constraint
real	threshold			# Threshold constraint
bool	debug				# Print debugging information?

int	i, j, nthreshold
real	minval, maxval, lcut

begin
	# Set the y array to be the values at the peak positions.
	do i = 1, npoints {
	    j = x[i]
	    y[i] = data[j]
	}

	# Determine the min and max values of the peaks.
	call alimr (y, npoints, minval, maxval)

	# Set the threshold based on the max of the absolute threshold and the
	# contrast.  Use arlt to set peaks below threshold to INDEF.
	if (!IS_INDEFR(threshold) || !IS_INDEFR(contrast)) {
	    if (IS_INDEFR(threshold))
		lcut = real (contrast * maxval)
	    else if (IS_INDEFR(contrast))
		lcut = real (threshold)
	    else
		lcut = max (real (threshold), real (contrast * maxval))
	    call arltr (y, npoints, lcut, INDEFR)
	}

	if (debug) {
	    call printf ("  Highest peak value = %g.\n")
		call pargr (maxval)
	    call printf ("  Peak cutoff threshold = %g.\n")
		call pargr (lcut)
	    do i = 1, npoints {
		if (IS_INDEFR (y[i])) {
		    j = x[i]
		    call printf (
			"  Peak at column %d with value %g below threshold.\n")
			call pargi (j)
			call pargr (data[j])
		}
	    }
	}

	# Determine the number of acceptable peaks & resort the x and y arrays.
	nthreshold = 0
	do i = 1, npoints {
	    if (IS_INDEFR (y[i]))
		next
	    nthreshold = nthreshold + 1
	    x[nthreshold] = x[i]
	    y[nthreshold] = y[i]
	}

	if (debug) {
	    call printf ("  Number of peaks above the threshold = %d.\n")
		call pargi (nthreshold)
	}

	return (nthreshold)
end

# FIND_ISOLATED -- Flag peaks which are within separation of a peak
# with a higher peak value.
#
# The peak positions, x, and their ranks, rank, are input.
# The rank array contains the indices of the peak positions in order from
# the highest peak value to the lowest peak value.  Starting with
# highest rank (rank[1]) all peaks of lower rank within separation
# are marked by setting their positions to INDEF.  The number of
# unflaged peaks is returned.

int procedure find_isolated (x, rank, npoints, separation, debug)

# Procedure parameters:
real	x[npoints]		# Positions of points
int	rank[npoints]		# Rank of peaks
int	npoints			# Number of peaks
int	separation		# Minimum allowed separation
bool	debug			# Print diagnostic information

int	i, j
int	nisolated

begin
	# Eliminate close neighbors.  The eliminated
	# peaks are marked by setting their positions to INDEF.
	nisolated = 0
	do i = 1, npoints {
	    if (IS_INDEFR (x[rank[i]]))
		next
	    nisolated = nisolated + 1
	    do j = i + 1, npoints {
		if (IS_INDEFR (x[rank[j]]))
		    next
		if (abs (x[rank[i]] - x[rank[j]]) < separation) {
		    if (debug) {
			call printf (
			    "  Peak at column %d too near peak at column %d.\n")
			    call pargi (int (x[rank[j]]))
			    call pargi (int (x[rank[i]]))
		    }
		    x[rank[j]] = INDEFR
		}
	    }
	}

	if (debug) {
	    call printf ("  Number of peaks separated by %d pixels = %d.\n")
		call pargi (separation)
		call pargi (nisolated)
	}

	# Return number of isolated peaks.
	return (nisolated)
end


# FIND_NMAX -- Select up to the nmax highest ranked peaks.
#
# The data values, data, peak positions, x, and their ranks, rank, are input.
# The data values are used only in printing debugging information.
# Peak positions previously eliminated are flaged by the value INDEF.
# The rank array contains the indices to the peak positions in order from
# the highest peak value to the lowest peak value.
# First all but the nmax highest ranked peaks (which have not been previously
# eliminated) are eliminated by marking their positions with the value INDEF.
# Then the remaining peaks are resorted to contain only the unflaged
# peaks and the number of such peaks is returned.

int procedure find_nmax (data, x, rank, npoints, nmax, debug)

real	data[ARB]			# Input data values
real	x[npoints]			# Peak positions
int	rank[npoints]			# Ranks of peaks
int	npoints				# Number of input peaks
int	nmax				# Max number of peaks to be selected
bool	debug				# Print debugging information?

int	i, j, npeaks

begin
	# Only mark peaks to reject if the number peaks is greater than nmax.
	if (nmax < npoints) {
	    npeaks = 0
	    do i = 1, npoints {
	        if (IS_INDEFR (x[rank[i]]))
		    next
	        npeaks = npeaks + 1
	        if (npeaks > nmax) {
		    if (debug) {
			j = x[rank[i]]
			call printf (
		    "  Reject peak at column %d with rank %d and value %g.\n")
			    call pargi (j)
			    call pargi (i)
			    call pargr (data[j])
		    }
	            x[rank[i]] = INDEFR
	        }
	    }
	}

	# Eliminate INDEF points and determine the number of spectra found.
	npeaks = 0
	do i = 1, npoints {
	    if (IS_INDEFR (x[i]))
		next
	    npeaks = npeaks + 1
	    x[npeaks] = x[i]
	}

	return (npeaks)
end


# FIND_UNMAX -- Select up to the nmax highest ranked peaks in bins.
#
# The data values, data, peak positions, x, and their ranks, rank, are input.
# The data values are used only in printing debugging information.
# Peak positions previously eliminated are flaged by the value INDEF.
# The rank array contains the indices to the peak positions in order from
# the highest peak value to the lowest peak value.
# First all but the nmax highest ranked peaks (which have not been previously
# eliminated) are eliminated by marking their positions with the value INDEF.
# Then the remaining peaks are resorted to contain only the unflaged
# peaks and the number of such peaks is returned.

int procedure find_unmax (data, npts, x, rank, npoints, nmax, nbins, debug)

real	data[npts]			# Input data values
int	npts				# Number of input data points
real	x[npoints]			# Peak positions
int	rank[npoints]			# Ranks of peaks
int	npoints				# Number of input peaks
int	nmax				# Max number of peaks to be selected
int	nbins				# Number of sample bins
bool	debug				# Print debugging information?

int	i, j, npeaks, width, x1, x2
real	a

begin
	# Only mark peaks to reject if the number peaks is greater than nmax.
	if (nmax < npoints) {

	    # Set up circular bins and select highest peak in each bin
	    # until the desired number of peaks is selected.

	    width = min (npts-1, nint ((npts-1) / (nbins-.5)))
	    x2 = 1
	    npeaks = 0
	    repeat {
		x1 = x2
		x2 = mod (x1 + width, npts) + 1
		j = 0
		do i = 1, npoints {
		    a = x[rank[i]]
		    if (IS_INDEFR (a) || a < 0) {
			j = j + 1
			next
		    }
		    if (x1 < x2) {
			if (a >= x1 && a <= x2) {
			    x[rank[i]] = -a
			    npeaks = npeaks + 1
			    break
			}
		    } else {
			if (a <= x2 || a >= x1) {
			    x[rank[i]] = -a
			    npeaks = npeaks + 1
			    break
			}
		    }
		}
	    } until (npeaks >= nmax || j == npoints)

	    # Now eliminate all unused peaks and reset the selected peaks.
	    do i = 1, npoints {
		if (!IS_INDEFR (x[i]) && x[i] < 1)
		    x[i] = -x[i]
		else
		    x[i] = INDEFR
	    }
	}

	# Eliminate INDEF points and determine the number of peaks found.
	npeaks = 0
	do i = 1, npoints {
	    if (IS_INDEFR (x[i]))
		next
	    npeaks = npeaks + 1
	    x[npeaks] = x[i]
	}

	return (npeaks)
end


# COMPARE -- Compare procedure for sort used in FIND_PEAKS.
# Larger values are indexed first.  INDEF values are indexed last.

int procedure compare (index1, index2)

# Procedure parameters:
int	index1				# Comparison index
int	index2				# Comparison index

pointer	y

common	/sort/ y

begin
	# INDEF points are considered to be smallest possible values.
	if (IS_INDEFR (Memr[y - 1 + index1]))
	    return (1)
	else if (IS_INDEFR (Memr[y - 1 + index2]))
	    return (-1)
	else if (Memr[y - 1 + index1] < Memr[y - 1 + index2])
	    return (1)
	else if (Memr[y - 1 + index1] > Memr[y - 1 + index2])
	    return (-1)
	else
	    return (0)
end