# PEAKS -- The following procedures are general numerical functions
# dealing with finding peaks in a data array.
#
# FIND_PEAKS		Find the 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.
# COMPARE		Compare procedure for sort used in FIND_PEAKS.

# FIND_PEAKS -- Find the 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:
double	data[npoints]		# Input data array
double	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	i, j
int	nlmax, nthreshold, nisolated, npeaks
pointer	sp, y, rank

int	find_local_maxima(), find_threshold(), find_isolated(), find_nmax()
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_DOUBLE)

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

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

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

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

	call sfree (sp)
	return (npeaks)
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)

double	data[npoints]			# Input data array
double	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
double	data[npoints]		# Data values
int	npoints			# Number of points in the data vector

int	i, j, nright, nleft

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

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

	# Test for failure to be a local maxima
	if ((i == 0) && (j == npoints)) {
	    return (FALSE)		# Data is constant
	} else if (i == 0) {
	    if (data[j] > data[index])
	        return (FALSE)		# Data increases to right
	} else if (j == npoints) {
	    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)

double	data[ARB]			# Input data values
double	x[npoints]			# Input/Output peak positions
double	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
double	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 alimd (y, npoints, minval, maxval)

	# Set the threshold based on the max of the absolute threshold and the
	# contrast.  Use arltd to set peaks below threshold to INDEFD.
	lcut = max (double (threshold), double (contrast * maxval))
	call arltd (y, npoints, lcut, INDEFD)

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

	# Determine the number of acceptable peaks & resort the x and y arrays.
	nthreshold = 0
	do i = 1, npoints {
	    if (IS_INDEFD (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 INDEFD.  The number of
# unflaged peaks is returned.

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

# Procedure parameters:
double	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 INDEFD.
	nisolated = 0
	do i = 1, npoints {
	    if (IS_INDEFD (x[rank[i]]))
		next
	    nisolated = nisolated + 1
	    do j = i + 1, npoints {
		if (IS_INDEFD (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]] = INDEFD
		}
	    }
	}

	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 INDEFD.
# 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 INDEFD.
# 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)

double	data[ARB]			# Input data values
double	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_INDEFD (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 pargd (data[j])
		    }
	            x[rank[i]] = INDEFD
	        }
	    }
	}

	# Eliminate INDEFD points and determine the number of spectra found.
	npeaks = 0
	do i = 1, npoints {
	    if (IS_INDEFD (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.  INDEFD 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
	# INDEFD points are considered to be smallest possible values.
	if (IS_INDEFD (Memd[y - 1 + index1]))
	    return (1)
	else if (IS_INDEFD (Memd[y - 1 + index2]))
	    return (-1)
	else if (Memd[y - 1 + index1] < Memd[y - 1 + index2])
	    return (1)
	else if (Memd[y - 1 + index1] > Memd[y - 1 + index2])
	    return (-1)
	else
	    return (0)
end