sddspeakfind.c File Reference

Program to find peaks in a column of SDDS data as a function of row index. More...

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"

Go to the source code of this file.

Macros
#define	THRESHOLD   0x0001UL
#define	FIVEPOINT   0x0002UL
#define	CHANGETHRES   0x0004UL
#define	EZONEFRAC   0x0008UL
#define	PAR_THRESHOLD   0x0010UL
#define	SEVENPOINT   0x0020UL
#define	CURVATURELIMIT   0x0040UL

Enumerations
enum	option_type {   CLO_THRESHOLD , CLO_FIVEPOINT , CLO_CHANGETHRESHOLD , CLO_PIPE ,   CLO_COLUMN , CLO_EXCLUSIONZONE , CLO_MAJOR_ORDER , CLO_SEVENPOINT ,   CLO_CURVATURE_LIMIT , N_OPTIONS }

Functions
void	mark_peaks (double *data, int32_t *row_flag, int64_t rows, long n_points)
void	unmark_flat_peaks (double *data, int32_t *row_flag, int64_t rows, double change_threshold, long five_point, long seven_point, double *x_data, double curvature_limit)
void	unmark_excluded_peaks (double *data, int32_t *row_flag, int64_t rows, double ezone_fraction)
int	main (int argc, char **argv)

Variables
static char *	option [N_OPTIONS]
static char *	USAGE

Detailed Description

Program to find peaks in a column of SDDS data as a function of row index.

This program processes an SDDS (Self Describing Data Set) file and identifies peaks in a specified column based on certain criteria such as threshold, exclusion zones, and curvature limits. The program supports both input and output from files or pipes.

Features:

• Detects peaks using 3-point, 5-point, or 7-point criteria.
• Applies optional thresholds, curvature limits, and exclusion zones to filter peaks.
• Outputs the processed data to an SDDS file.
• Supports row-major and column-major orders for data storage.

Options:

• -threshold: Set a threshold for peak detection.
• -fivePoint or -sevenPoint: Use 5-point or 7-point methods for peak detection.
• -changeThreshold: Apply a fractional change threshold to flat peaks.
• -curvatureLimit: Apply curvature filtering to remove shallow peaks.
• -exclusionZone: Exclude peaks too close to each other.
• -column: Specify the column to analyze.
• -pipe: Use input/output pipes for SDDS data.
• -majorOrder: Specify row or column major order for the output.

Usage:

sddspeakfind [<inputfile>] [<outputfile>] [-pipe=[input][,output]] \
   -column=<columnName> \
   [-threshold=<value>|@<parametername>] \
   [{-fivePoint|-sevenPoint}] \
   [-exclusionZone=<fraction>] \
   [-changeThreshold=<fractionalChange>] \
   [-curvatureLimit=<xColumn>,<minValue>] \
   [-majorOrder=row|column]

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, R. Soliday, H. Shang

Definition in file sddspeakfind.c.

Macro Definition Documentation

CHANGETHRES

#define CHANGETHRES   0x0004UL

Definition at line 54 of file sddspeakfind.c.

CURVATURELIMIT

#define CURVATURELIMIT   0x0040UL

Definition at line 58 of file sddspeakfind.c.

EZONEFRAC

#define EZONEFRAC   0x0008UL

Definition at line 55 of file sddspeakfind.c.

FIVEPOINT

#define FIVEPOINT   0x0002UL

Definition at line 53 of file sddspeakfind.c.

PAR_THRESHOLD

#define PAR_THRESHOLD   0x0010UL

Definition at line 56 of file sddspeakfind.c.

SEVENPOINT

#define SEVENPOINT   0x0020UL

Definition at line 57 of file sddspeakfind.c.

THRESHOLD

#define THRESHOLD   0x0001UL

Definition at line 52 of file sddspeakfind.c.

Enumeration Type Documentation

option_type

enum option_type

Definition at line 61 of file sddspeakfind.c.

                 {
  CLO_THRESHOLD,
  CLO_FIVEPOINT,
  CLO_CHANGETHRESHOLD,
  CLO_PIPE,
  CLO_COLUMN,
  CLO_EXCLUSIONZONE,
  CLO_MAJOR_ORDER,
  CLO_SEVENPOINT,
  CLO_CURVATURE_LIMIT,
  N_OPTIONS
};

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 105 of file sddspeakfind.c.

                                {
  SDDS_DATASET in_set, out_set;
  SCANNED_ARG *s_arg;
  long i_arg, page_returned;
  int64_t rows, row;
  int32_t *row_flag;
  char *input, *output, *column_name, *par_threshold_name, *x_column_name;
  double *data, *x_data;
  unsigned long pipe_flags, flags, major_order_flag;
  double threshold, ezone_fraction, change_threshold, curvature_limit;
  short column_major_order = -1;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&s_arg, argc, argv);
  if (argc < 2 || argc > (2 + N_OPTIONS))
    bomb(NULL, USAGE);
 
  flags = pipe_flags = 0;
  input = output = NULL;
  column_name = NULL;
  ezone_fraction = change_threshold = curvature_limit = 0;
  row_flag = NULL;
  par_threshold_name = NULL;
  x_column_name = NULL;
  x_data = NULL;
 
  for (i_arg = 1; i_arg < argc; i_arg++) {
    if (s_arg[i_arg].arg_type == OPTION) {
      switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
      case CLO_MAJOR_ORDER:
        major_order_flag = 0;
        s_arg[i_arg].n_items--;
        if (s_arg[i_arg].n_items > 0 &&
            (!scanItemList(&major_order_flag, s_arg[i_arg].list + 1,
                           &s_arg[i_arg].n_items, 0, "row", -1, NULL, 0,
                           SDDS_ROW_MAJOR_ORDER, "column", -1, NULL, 0,
                           SDDS_COLUMN_MAJOR_ORDER, NULL)))
          SDDS_Bomb("invalid -majorOrder syntax/values");
        if (major_order_flag & SDDS_COLUMN_MAJOR_ORDER)
          column_major_order = 1;
        else if (major_order_flag & SDDS_ROW_MAJOR_ORDER)
          column_major_order = 0;
        break;
      case CLO_THRESHOLD:
        if (s_arg[i_arg].n_items == 2) {
          if (s_arg[i_arg].list[1][0] == '@') {
            SDDS_CopyString(&par_threshold_name, s_arg[i_arg].list[1] + 1);
            flags |= PAR_THRESHOLD;
          } else {
            if (sscanf(s_arg[i_arg].list[1], "%lf", &threshold) != 1)
              SDDS_Bomb("incorrect -threshold syntax");
            flags |= THRESHOLD;
          }
        } else
          SDDS_Bomb("incorrect -threshold syntax");
        break;
      case CLO_FIVEPOINT:
        flags |= FIVEPOINT;
        break;
      case CLO_SEVENPOINT:
        flags |= SEVENPOINT;
        break;
      case CLO_CHANGETHRESHOLD:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%lf", &change_threshold) != 1 ||
            change_threshold <= 0)
          SDDS_Bomb("incorrect -changeThreshold syntax or values");
        flags |= CHANGETHRES;
        break;
      case CLO_CURVATURE_LIMIT:
        if (s_arg[i_arg].n_items != 3 ||
            !strlen(x_column_name = s_arg[i_arg].list[1]) ||
            sscanf(s_arg[i_arg].list[2], "%lf", &curvature_limit) != 1 ||
            curvature_limit <= 0)
          SDDS_Bomb("incorrect -curvatureLimit syntax or values");
        flags |= CURVATURELIMIT;
        break;
      case CLO_PIPE:
        if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipe_flags))
          SDDS_Bomb("invalid -pipe syntax");
        break;
      case CLO_COLUMN:
        if (s_arg[i_arg].n_items != 2)
          SDDS_Bomb("invalid -column syntax");
        column_name = s_arg[i_arg].list[1];
        break;
      case CLO_EXCLUSIONZONE:
        if (s_arg[i_arg].n_items != 2 ||
            sscanf(s_arg[i_arg].list[1], "%lf", &ezone_fraction) != 1 ||
            ezone_fraction <= 0)
          SDDS_Bomb("invalid -exclusionZone syntax or value");
        flags |= EZONEFRAC;
        break;
      default:
        fprintf(stderr, "error: unknown/ambiguous option: %s\n",
                s_arg[i_arg].list[0]);
        exit(1);
        break;
      }
    } else {
      if (input == NULL)
        input = s_arg[i_arg].list[0];
      else if (output == NULL)
        output = s_arg[i_arg].list[0];
      else
        SDDS_Bomb("too many filenames");
    }
  }
 
  processFilenames("sddspeakfind", &input, &output, pipe_flags, 0, NULL);
 
  if (!column_name)
    SDDS_Bomb("-column option must be given");
 
  if (!SDDS_InitializeInput(&in_set, input))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  if (!SDDS_FindColumn(&in_set, FIND_NUMERIC_TYPE, column_name, NULL))
    SDDS_Bomb("the given column is nonexistent or nonnumeric");
 
  if (x_column_name &&
      !SDDS_FindColumn(&in_set, FIND_NUMERIC_TYPE, x_column_name, NULL))
    SDDS_Bomb("the given x column is nonexistent or nonnumeric");
 
  if (!SDDS_InitializeCopy(&out_set, &in_set, output, "w"))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  out_set.layout.data_mode.column_major =
    column_major_order != -1 ? column_major_order : in_set.layout.data_mode.column_major;
 
  if (!SDDS_WriteLayout(&out_set))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  /* Main loop for processing pages */
  while ((page_returned = SDDS_ReadPage(&in_set)) > 0) {
    if (!SDDS_CopyPage(&out_set, &in_set)) {
      SDDS_SetError("Problem copying data for output file");
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
 
    if ((rows = SDDS_CountRowsOfInterest(&out_set)) > 1) {
      data = SDDS_GetColumnInDoubles(&in_set, column_name);
      if (!data)
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
      if (x_column_name)
        x_data = SDDS_GetColumnInDoubles(&in_set, x_column_name);
      if (x_column_name && !x_data)
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
      row_flag = SDDS_Realloc(row_flag, sizeof(*row_flag) * rows);
      for (row = 0; row < rows; row++)
        row_flag[row] = 0;
 
      mark_peaks(data, row_flag, rows,
                 flags & SEVENPOINT ? 7 : (flags & FIVEPOINT ? 5 : 3));
 
      if (flags & PAR_THRESHOLD)
        if (!SDDS_GetParameter(&in_set, par_threshold_name, &threshold))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
      if (flags & THRESHOLD || flags & PAR_THRESHOLD)
        for (row = 0; row < rows; row++)
          if (row_flag[row] && data[row] < threshold)
            row_flag[row] = 0;
 
      if (flags & CHANGETHRES || flags & CURVATURELIMIT)
        unmark_flat_peaks(data, row_flag, rows, change_threshold,
                          flags & FIVEPOINT, flags & SEVENPOINT, x_data, curvature_limit);
 
      if (flags & EZONEFRAC)
        unmark_excluded_peaks(data, row_flag, rows, ezone_fraction);
 
      if (!SDDS_AssertRowFlags(&out_set, SDDS_FLAG_ARRAY, row_flag, rows))
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
      free(data);
    }
 
    if (!SDDS_WritePage(&out_set))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  if (!SDDS_Terminate(&in_set) || !SDDS_Terminate(&out_set)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(1);
  }
 
  free_scanargs(&s_arg, argc);
  free(par_threshold_name);
  free(row_flag);
 
  return 0;
}

mark_peaks()

void mark_peaks	(	double *	data,
		int32_t *	row_flag,
		int64_t	rows,
		long	n_points )

Definition at line 300 of file sddspeakfind.c.

                                                                              {
  double y0, y1, y2;
  int64_t i, rows_minus;
  short five_point = 0;
 
  switch (n_points) {
  case 5:
    y0 = data[1];
    y1 = data[2];
    if (rows < 5)
      return;
    five_point = 1;
    break;
  case 7:
    if (rows < 7)
      return;
    for (i = 3; i < (rows - 3); i++) {
      if ((data[i - 3] < data[i - 2] && data[i - 2] < data[i - 1] &&
           data[i - 1] < data[i] && data[i + 1] < data[i] &&
           data[i + 2] < data[i + 1] && data[i + 3] < data[i + 2]))
        row_flag[i] = 1;
    }
    return;
  default:
    if (rows < 3)
      return;
    y0 = data[0];
    y1 = data[1];
  }
 
  rows_minus = rows - (five_point ? 2 : 1);
  for (i = (five_point ? 2 : 1); i < rows_minus; i++) {
    y2 = data[i + 1];
    if ((y1 > y0 && y1 > y2) || (y1 == y0 && y1 > y2) || (y1 > y0 && y1 == y2)) {
      if (!five_point || (data[i - 2] < y0 && data[i + 2] < y2))
        row_flag[i] = 1;
    }
    y0 = y1;
    y1 = y2;
  }
}

unmark_excluded_peaks()

void unmark_excluded_peaks	(	double *	data,
		int32_t *	row_flag,
		int64_t	rows,
		double	ezone_fraction )

Definition at line 377 of file sddspeakfind.c.

                                                  {
  int64_t lower, upper, i, j, rows1;
  rows1 = rows - 1;
 
  for (i = 0; i < rows; i++) {
    if (row_flag[i]) {
      lower = i - ((int64_t)(ezone_fraction / 2.0 * rows + 0.5));
      if (lower < 0)
        lower = 0;
      upper = i + ((int64_t)(ezone_fraction / 2.0 * rows + 0.5));
      if (upper > rows1)
        upper = rows1;
      for (j = lower; j <= upper; j++) {
        if (j != i && row_flag[j] && data[j] <= data[i])
          row_flag[j] = 0;
      }
    }
  }
}

unmark_flat_peaks()

void unmark_flat_peaks	(	double *	data,
		int32_t *	row_flag,
		int64_t	rows,
		double	change_threshold,
		long	five_point,
		long	seven_point,
		double *	x_data,
		double	curvature_limit )

Definition at line 342 of file sddspeakfind.c.

                                                               {
  long delta;
  int64_t i, rows1;
 
  delta = five_point ? 2 : (seven_point ? 3 : 1);
  rows1 = rows - delta;
 
  if (change_threshold > 0) {
    for (i = delta; i < rows1; i++) {
      if (row_flag[i]) {
        if (!data[i] || (data[i] - data[i + delta]) / data[i] > change_threshold)
          continue;
        if (!data[i] || (data[i] - data[i - delta]) / data[i] > change_threshold)
          continue;
        row_flag[i] = 0;
      }
    }
  }
 
  if (curvature_limit > 0) {
    double chi, coef[3], scoef[3], diff[7];
    int32_t order[3] = {0, 1, 2};
    for (i = delta; i < rows1; i++) {
      if (row_flag[i]) {
        lsfg(x_data + i - delta, data + i - delta, NULL, 2 * delta + 1,
             3, order, coef, scoef, &chi, diff, ipower);
        if (fabs(coef[2]) < curvature_limit)
          row_flag[i] = 0;
      }
    }
  }
}

Variable Documentation

option

char* option[N_OPTIONS]

static

Initial value:

= {
  "threshold",
  "fivepoint",
  "changethreshold",
  "pipe",
  "column",
  "exclusionzone",
  "majorOrder",
  "sevenpoint",
  "curvatureLimit"}

Definition at line 74 of file sddspeakfind.c.

                                 {
  "threshold",
  "fivepoint",
  "changethreshold",
  "pipe",
  "column",
  "exclusionzone",
  "majorOrder",
  "sevenpoint",
  "curvatureLimit"};

USAGE

char* USAGE

static

Initial value:

=
    "sddspeakfind [<inputfile>] [<outputfile>] [-pipe=[input][,output]] \n"
    "   -column=<columnName> \n"
    "  [-threshold=<value>|@<parametername>] \n"
    "  [{-fivePoint|-sevenPoint}] \n"
    "  [-exclusionZone=<fraction>] \n"
    "  [-changeThreshold=<fractionalChange>] \n"
    "  [-curvatureLimit=<xColumn>,<minValue>] \n"
    "  [-majorOrder=row|column] \n\n"
    "Finds peaks in a column of data as a function of row index.\n"
    "Program by Michael Borland. ("__DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n"

Definition at line 85 of file sddspeakfind.c.