Baseline subtraction tool for SDDS datasets. More...

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

Macros
#define	SELECT_ENDPOINTS 0x0001U

#define	SELECT_OUTSIDEFWHA 0x0002U

#define	SELECT_ANTIOUTLIER 0x0004U

#define	METHOD_FIT 0x0001U

#define	METHOD_AVERAGE 0x0002U

Enumerations
enum	option_type { CLO_PIPE , CLO_COLUMNS , CLO_METHOD , CLO_SELECT , CLO_NONNEGATIVE , CLO_REPEATS , CLO_DESPIKE , CLO_MAJOR_ORDER , N_OPTIONS }

Functions
long	resolveColumnNames (SDDS_DATASET SDDSin, char *column, int32_t columns)

void	selectEndpoints (short *selected, int64_t rows, long endPoints)

void	selectOutsideFWHA (double data, double indepData, short *selected, int64_t rows, double fwhaLimit)

void	selectAntiOutlier (double data, short selected, int64_t rows, long passes)

void	fitAndRemoveBaseline (double data, double indepData, short *selected, int64_t rows, long terms)

void	averageAndRemoveBaseline (double data, short selected, int64_t rows)

void	despikeProfile (double *data, int64_t rows, long widthLimit, long passes)

int	main (int argc, char **argv)

Variables
static char *	option [N_OPTIONS]

char *	USAGE

Detailed Description

Baseline subtraction tool for SDDS datasets.

This program processes SDDS datasets to subtract a baseline from specified columns. The baseline can be computed using various methods and selection criteria. Options allow for nonnegative constraints, despiking, and multiple repeats.

Usage: sddsbaseline [<input>] [<output>] [-pipe=<in>[,<out>]] [-columns=<listOfNames>] [-nonnegative [-despike=passes=<number>,widthlimit=

] [-repeats=<count>]] [-select={endpoints=<number> | outsideFWHA=<multiplier> | antioutlier=<passes>}] [-method={average|fit[,terms=<number>]}] [-majorOrder=row|column]

Options: -pipe Specify input and/or output pipes. -columns List of columns to process. -nonnegative Forces all values to be nonnegative after baseline subtraction. This is accomplished by setting all negative values to 0. -despike Specify that positive features narrower than widthLimit shall be set to zero. Parameters: passes=<number> Number of despike passes. widthlimit=

Width limit for despiking. -repeats Specify how many times to apply the baseline removal algorithm. Meaningful only if used in combination with -nonnegative. -select Specify how to select points to include in baseline determination. Options: endpoints=<number> outsideFWHA=<multiplier> antioutlier=<passes> -method Specify how to process selected points in order to compute baseline. Options: average fit[,terms=<number>] -majorOrder Specify write output in row or column major order.

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, R. Soliday, H. Shang

Definition in file sddsbaseline.c.

Macro Definition Documentation

◆ METHOD_AVERAGE

#define METHOD_AVERAGE 0x0002U

Definition at line 115 of file sddsbaseline.c.

◆ METHOD_FIT

#define METHOD_FIT 0x0001U

Definition at line 114 of file sddsbaseline.c.

◆ SELECT_ANTIOUTLIER

#define SELECT_ANTIOUTLIER 0x0004U

Definition at line 112 of file sddsbaseline.c.

◆ SELECT_ENDPOINTS

#define SELECT_ENDPOINTS 0x0001U

Definition at line 110 of file sddsbaseline.c.

◆ SELECT_OUTSIDEFWHA

#define SELECT_OUTSIDEFWHA 0x0002U

Definition at line 111 of file sddsbaseline.c.

Enumeration Type Documentation

◆ option_type

enum option_type

Definition at line 56 of file sddsbaseline.c.

                 {
  CLO_PIPE,
  CLO_COLUMNS,
  CLO_METHOD,
  CLO_SELECT,
  CLO_NONNEGATIVE,
  CLO_REPEATS,
  CLO_DESPIKE,
  CLO_MAJOR_ORDER,
  N_OPTIONS
};

Function Documentation

◆ averageAndRemoveBaseline()

void averageAndRemoveBaseline	(	double *	data,
		short *	selected,
		int64_t	rows )

Definition at line 472 of file sddsbaseline.c.

                                                                           {
  int64_t i, count;
  double ave;
  for (i = count = ave = 0; i < rows; i++)
    if (selected[i]) {
      count++;
      ave += data[i];
    }
 
  if (count) {
    ave /= count;
    for (i = 0; i < rows; i++)
      data[i] -= ave;
  }
}

◆ despikeProfile()

void despikeProfile	(	double *	data,
		int64_t	rows,
		long	widthLimit,
		long	passes )

Definition at line 488 of file sddsbaseline.c.

                                                                              {
  int64_t i, i0, i1, j;
  while (--passes >= 0) {
    for (i = 0; i < rows; i++) {
      if (i == 0) {
        if (data[i] != 0)
          continue;
      } else if (!(data[i - 1] == 0 && data[i] != 0))
        continue;
      i0 = i;
      for (i1 = i + 1; i1 < rows; i1++)
        if (!data[i1])
          break;
      if ((i1 - i0) <= widthLimit)
        for (j = i0; j < i1; j++)
          data[j] = 0;
    }
  }
}

◆ fitAndRemoveBaseline()

void fitAndRemoveBaseline	(	double *	data,
		double *	indepData,
		short *	selected,
		int64_t	rows,
		long	terms )

Definition at line 431 of file sddsbaseline.c.

                                                                                                           {
  double *data, *indepData = NULL;
  int64_t count, i, j;
  double chi;
  double *coef, *sCoef;
 
  if (!(coef = malloc(sizeof(*coef) * fitTerms)) || !(sCoef = malloc(sizeof(*sCoef) * fitTerms)))
    SDDS_Bomb("allocation failure (fitAndRemoveBaseline)");
  for (i = count = 0; i < rows; i++)
    if (selected[i])
      count++;
  if (count < 3)
    return;
  if (!(data = malloc(sizeof(*data) * count)) || !(indepData = malloc(sizeof(*indepData) * count)))
    SDDS_Bomb("memory allocation failure");
  for (i = j = 0; i < rows; i++) {
    if (selected[i]) {
      data[j] = data0[i];
      indepData[j] = indepData0[i];
      j++;
    }
  }
 
  if (!lsfn(indepData, data, NULL, (long)count, fitTerms - 1, coef, sCoef, &chi, NULL))
    return;
 
  for (i = 0; i < rows; i++) {
    double term;
    term = 1;
    for (j = 0; j < fitTerms; j++) {
      data0[i] -= term * coef[j];
      term *= indepData0[i];
    }
  }
 
  free(data);
  free(indepData);
  free(coef);
  free(sCoef);
}

◆ main()

int main	(	int	argc,
		char **	argv )

Definition at line 125 of file sddsbaseline.c.

                                {
  int iArg;
  char **inputColumn;
  char *input, *output;
  long readCode, repeats, repeat, fitTerms = 2;
  int64_t rows, i, j;
  int32_t columns;
  unsigned long pipeFlags, methodFlags, selectFlags, dummyFlags, majorOrderFlag;
  SCANNED_ARG *scanned;
  SDDS_DATASET SDDSin, SDDSout;
  double *data, *indepData;
  int32_t endPoints, antiOutlierPasses;
  short *selected;
  double outsideFWHA;
  long nonnegative;
  int32_t despikePasses, despikeWidthLimit;
  short columnMajorOrder = -1;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&scanned, argc, argv);
  if (argc < 2)
    bomb(NULL, USAGE);
 
  output = input = NULL;
  inputColumn = NULL;
  columns = nonnegative = 0;
  repeats = 1;
  pipeFlags = methodFlags = selectFlags = dummyFlags = 0;
  endPoints = antiOutlierPasses = 0;
  outsideFWHA = 0;
  despikePasses = 0;
  despikeWidthLimit = 2;
 
  for (iArg = 1; iArg < argc; iArg++) {
    if (scanned[iArg].arg_type == OPTION) {
      /* process options here */
      switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
      case CLO_MAJOR_ORDER:
        majorOrderFlag = 0;
        scanned[iArg].n_items -= 1;
        if (scanned[iArg].n_items > 0 &&
            (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].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 (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
          columnMajorOrder = 1;
        else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
          columnMajorOrder = 0;
        break;
      case CLO_COLUMNS:
        if (scanned[iArg].n_items < 2)
          SDDS_Bomb("invalid -columns syntax");
        inputColumn = tmalloc(sizeof(*inputColumn) * (columns = scanned[iArg].n_items - 1));
        for (i = 0; i < columns; i++)
          inputColumn[i] = scanned[iArg].list[i + 1];
        break;
      case CLO_PIPE:
        if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
          SDDS_Bomb("invalid -pipe syntax");
        break;
      case CLO_METHOD:
        if (!(scanned[iArg].n_items -= 1))
          SDDS_Bomb("invalid -method syntax");
        if (!scanItemList(&methodFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "average", -1, NULL, 0, METHOD_AVERAGE,
                          "fit", -1, NULL, 0, METHOD_FIT,
                          "terms", SDDS_LONG, &fitTerms, 1, 0, NULL) ||
            bitsSet(methodFlags) != 1 || fitTerms < 2)
          SDDS_Bomb("invalid -method syntax/values");
        break;
      case CLO_SELECT:
        if (!(scanned[iArg].n_items -= 1))
          SDDS_Bomb("invalid -select syntax");
        if (!scanItemList(&selectFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "endpoints", SDDS_LONG, &endPoints, 1, SELECT_ENDPOINTS,
                          "outsidefwha", SDDS_DOUBLE, &outsideFWHA, 1, SELECT_OUTSIDEFWHA,
                          "antioutlier", SDDS_LONG, &antiOutlierPasses, 1, SELECT_ANTIOUTLIER, NULL) ||
            bitsSet(selectFlags) != 1)
          SDDS_Bomb("invalid -select syntax/values");
        break;
      case CLO_NONNEGATIVE:
        nonnegative = 1;
        break;
      case CLO_REPEATS:
        if (scanned[iArg].n_items != 2 ||
            sscanf(scanned[iArg].list[1], "%ld", &repeats) != 1 ||
            repeats <= 0)
          SDDS_Bomb("invalid -repeats syntax");
        break;
      case CLO_DESPIKE:
        despikePasses = 1;
        if (!(scanned[iArg].n_items -= 1))
          SDDS_Bomb("invalid -despike syntax");
        if (!scanItemList(&dummyFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "passes", SDDS_LONG, &despikePasses, 1, 0,
                          "widthlimit", SDDS_LONG, &despikeWidthLimit, 1, 0, NULL) ||
            despikePasses < 0 ||
            despikeWidthLimit < 1)
          SDDS_Bomb("invalid -despike syntax/values");
        break;
      default:
        fprintf(stderr, "error: unknown/ambiguous option: %s\n", scanned[iArg].list[0]);
        exit(EXIT_FAILURE);
        break;
      }
    } else {
      if (!input)
        input = scanned[iArg].list[0];
      else if (!output)
        output = scanned[iArg].list[0];
      else
        SDDS_Bomb("too many filenames seen");
    }
  }
 
  if (!bitsSet(selectFlags))
    SDDS_Bomb("no -select option given");
  if (!bitsSet(methodFlags))
    SDDS_Bomb("no -method option given");
 
  if (!nonnegative && despikePasses)
    SDDS_Bomb("not meaningful to despike without setting -nonnegative");
  if (!nonnegative && repeats > 1)
    SDDS_Bomb("not meaningful to repeat without setting -nonnegative");
 
  processFilenames("sddsbaseline", &input, &output, pipeFlags, 0, NULL);
 
  if (!columns)
    SDDS_Bomb("supply the names of columns to process with the -columns option");
 
  if (!SDDS_InitializeInput(&SDDSin, input))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  if (!resolveColumnNames(&SDDSin, &inputColumn, &columns))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!columns)
    SDDS_Bomb("no columns selected for processing");
 
  if (!SDDS_InitializeCopy(&SDDSout, &SDDSin, output, "w"))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!SDDS_DefineSimpleColumn(&SDDSout, "SelectedForBaselineDetermination", NULL, SDDS_SHORT))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (columnMajorOrder != -1)
    SDDSout.layout.data_mode.column_major = columnMajorOrder;
  else
    SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
  if (!SDDS_WriteLayout(&SDDSout))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  indepData = NULL;
  selected = NULL;
  while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
    if (!SDDS_CopyPage(&SDDSout, &SDDSin))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    if ((rows = SDDS_CountRowsOfInterest(&SDDSin))) {
      if (!(indepData = SDDS_Realloc(indepData, sizeof(*indepData) * rows)) ||
          !(selected = SDDS_Realloc(selected, sizeof(*selected) * rows)))
        SDDS_Bomb("memory allocation failure");
      for (i = 0; i < rows; i++)
        indepData[i] = i;
      for (i = 0; i < columns; i++) {
        if (!(data = SDDS_GetColumnInDoubles(&SDDSin, inputColumn[i])))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
        for (repeat = 0; repeat < repeats; repeat++) {
          for (j = 0; j < rows; j++)
            selected[j] = 0;
          switch (selectFlags) {
          case SELECT_ENDPOINTS:
            selectEndpoints(selected, rows, endPoints);
            break;
          case SELECT_OUTSIDEFWHA:
            selectOutsideFWHA(data, indepData, selected, rows, outsideFWHA);
            break;
          case SELECT_ANTIOUTLIER:
            selectAntiOutlier(data, selected, rows, antiOutlierPasses);
            break;
          default:
            SDDS_Bomb("invalid select flag");
            break;
          }
          switch (methodFlags) {
          case METHOD_FIT:
            fitAndRemoveBaseline(data, indepData, selected, rows, fitTerms);
            break;
          case METHOD_AVERAGE:
            averageAndRemoveBaseline(data, selected, rows);
            break;
          default:
            break;
          }
          if (nonnegative) {
            for (j = 0; j < rows; j++)
              if (data[j] < 0)
                data[j] = 0;
            if (despikePasses)
              despikeProfile(data, rows, despikeWidthLimit, despikePasses);
          }
        }
        if (!SDDS_SetColumnFromDoubles(&SDDSout, SDDS_SET_BY_NAME, data, rows, inputColumn[i]) ||
            !SDDS_SetColumn(&SDDSout, SDDS_SET_BY_NAME, selected, rows, "SelectedForBaselineDetermination"))
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
        free(data);
      }
    }
    if (!SDDS_WritePage(&SDDSout))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (!SDDS_Terminate(&SDDSin)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    return EXIT_FAILURE;
  }
  if (!SDDS_Terminate(&SDDSout)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    return EXIT_FAILURE;
  }
  if (indepData)
    free(indepData);
  if (selected)
    free(selected);
  return EXIT_SUCCESS;
}

◆ resolveColumnNames()

long resolveColumnNames	(	SDDS_DATASET *	SDDSin,
		char ***	column,
		int32_t *	columns )

Definition at line 348 of file sddsbaseline.c.

                                                                                {
  long i;
 
  SDDS_SetColumnFlags(SDDSin, 0);
  for (i = 0; i < *columns; i++) {
    if (!SDDS_SetColumnsOfInterest(SDDSin, SDDS_MATCH_STRING, (*column)[i], SDDS_OR))
      return 0;
  }
  free(*column);
  if (!(*column = SDDS_GetColumnNames(SDDSin, columns)) || *columns == 0) {
    SDDS_SetError("no columns found");
    return 0;
  }
  return 1;
}

◆ selectAntiOutlier()

void selectAntiOutlier	(	double *	data,
		short *	selected,
		int64_t	rows,
		long	passes )

Definition at line 406 of file sddsbaseline.c.

                                                                                 {
  double ave, sum2, limit;
  int64_t i, count;
 
  for (i = 0; i < rows; i++)
    selected[i] = 1;
  while (--passes >= 0) {
    for (i = ave = count = 0; i < rows; i++)
      if (selected[i]) {
        ave += data[i];
        count++;
      }
    if (!count)
      break;
    ave /= count;
    for (i = sum2 = 0; i < rows; i++)
      if (selected[i])
        sum2 += sqr(data[i] - ave);
    limit = 2 * sqrt(sum2 / count);
    for (i = 0; i < rows; i++)
      if (selected[i] && fabs(data[i] - ave) > limit)
        selected[i] = 0;
  }
}

◆ selectEndpoints()

void selectEndpoints	(	short *	selected,
		int64_t	rows,
		long	endPoints )

Definition at line 364 of file sddsbaseline.c.

                                                                    {
  int64_t i;
  for (i = 0; i < endPoints && i < rows; i++)
    selected[i] = 1;
  for (i = rows - 1; i > rows - endPoints - 1 && i >= 0; i--)
    selected[i] = 1;
}

◆ selectOutsideFWHA()

void selectOutsideFWHA	(	double *	data,
		double *	indepData,
		short *	selected,
		int64_t	rows,
		double	fwhaLimit )

Definition at line 372 of file sddsbaseline.c.

                                                                                                         {
  double top, base, fwha;
  int64_t i, i1, i2, i2save;
  int64_t imin, imax;
  double point1, point2;
 
  if (rows < 3 || fwhaLimit <= 0)
    return;
 
  index_min_max(&imin, &imax, data, rows);
  if (!data[imax])
    return;
 
  if (!findTopBaseLevels(&top, &base, data, rows, 50, 2.0))
    return;
 
  if ((i1 = findCrossingPoint(0, data, rows, (top - base) * 0.5 + base, 1, indepData, &point1)) < 0 ||
      (i2 = i2save = findCrossingPoint(i1, data, rows, (top - base) * 0.9 + base, -1, NULL, NULL)) < 0 ||
      (i2 = findCrossingPoint(i2, data, rows, (top - base) * 0.5 + base, -1, indepData, &point2)) < 0) {
    return;
  }
  fwha = point2 - point1;
 
  for (i = 0; i < rows; i++)
    selected[i] = 1;
  for (i = imax - fwha * fwhaLimit; i <= imax + fwha * fwhaLimit; i++) {
    if (i < 0)
      continue;
    if (i >= rows)
      break;
    selected[i] = 0;
  }
}

Variable Documentation

◆ option

char* option[N_OPTIONS]

static

Initial value:

= {
  "pipe",
  "columns",
  "method",
  "select",
  "nonnegative",
  "repeats",
  "despike",
  "majorOrder",
}

Definition at line 68 of file sddsbaseline.c.

                                 {
  "pipe",
  "columns",
  "method",
  "select",
  "nonnegative",
  "repeats",
  "despike",
  "majorOrder",
};

◆ USAGE

char* USAGE