Extracting and Filtering Minima and Maxima of 1D Functions

Persistence1D is a class for finding local extrema and their persistence in one-dimensional data. Local minima and local maxima are extracted, paired, and sorted according to their persistence.

The code runs in O(n log n) time, where n is the number of input points.

Additionally, the function can be reconstructed such that it only contains a filtered set of local extrema. This is useful for smoothing the data.

Persistence1D has been written by Yeara Kozlov and Tino Weinkauf, Max Planck Institute for Informatics, Saarbrücken, Germany. You may use it as you wish, it is in the public domain. If you find it useful, it would be nice to hear from you. Just drop us a line.

example of the input and output of the Persistence1D code

Persistence1D detects all local minima and maxima. It pairs them according to persistence. Subsequent filtering reveals the most dominant local extrema of the function.

example of reconstructing a function based on a filtered set of extrema

Reconstructing a smoothed version of the input function that contains only the filtered set of extrema. This works only from within Matlab.

Overview

All relevant code is found in a single header file for simple inclusion in other projects. The only dependency is the C++ standard library.

Besides C++, the code can also be used from within Matlab via mex. There is also a command line program to quickly process text files with data.

Input

  • One-dimensional vector of float values. It is assumed that this represents the function values of a one-dimensional function.

Output

  • Indices of minima and maxima points.
  • List of all paired extrema and their persistence.
  • The global minimum is never paired and is returned separately.

If desired, the output can be filtered by a persistence threshold and adjusted to use the 1-indexing convention (e.g., for Matlab).

Example for C++

The following example shows how simple it is to use this class. It is practically one function call.

#include "persistence1d.hpp"

using namespace std;
using namespace p1d;

int main()
{
    //Create some data
    vector< float > data;
    data.push_back(2.0);   data.push_back(5.0);   data.push_back(7.0);
    data.push_back(-12.0); data.push_back(-13.0); data.push_back(-7.0);
    data.push_back(10.0);  data.push_back(18.0);  data.push_back(6.0);
    data.push_back(8.0);   data.push_back(7.0);   data.push_back(4.0);

    //Run persistence on data - this is the main call.
    Persistence1D p;
    p.RunPersistence(data);

    //Get all extrema with a persistence larger than 10.
    vector< TPairedExtrema > Extrema;
    p.GetPairedExtrema(Extrema, 10);

    //Print all found pairs - pairs are sorted ascending wrt. persistence.
    for(vector< TPairedExtrema >::iterator it = Extrema.begin(); it != Extrema.end(); it++)
    {
        cout << "Persistence: " << (*it).Persistence
             << " minimum index: " << (*it).MinIndex
             << " maximum index: " << (*it).MaxIndex
             << std::endl;
    }

    //Also, print the global minimum.
    cout << "Global minimum index: " << p.GetGlobalMinimumIndex()
         << " Global minimum value: " << p.GetGlobalMinimumValue() << endl;

    /*
        Note that filtering and printing can also be done with one single function call:
        p.PrintResults(10);
    */

    return 0;
}

Running the code from above yields this output:

Persistence: 14 minimum index: 11 maximum index: 7
Global minimum index: 4 Global minimum value: -13

Example for Matlab

The matlab interface is just as convenient and the results are easily plotted. The result of this script is one of the plots from above.

% Compile the mex file.
mex run_persistence1d.cpp

% Generate some data
x = 1:600;
SineLowFreq = sin(x * 0.01 * pi);
SineMedFreq = 0.25 * sin(x * 0.01 * pi * 4.9);
SineHighFreq = 0.15 * sin(x * 0.01 * pi * 12.1);
data = SineLowFreq + SineMedFreq + SineHighFreq;

% Input should be in SINGLE precision
single_precision_data = single(data);

% Call persistence1d - this is the main call
[minIndices maxIndices persistence globalMinIndex globalMinValue] = run_persistence1d(single_precision_data);

% Set a threshold for filtering the result
threshold = 0.5;

% Use filter_features_by_persistence to filter the pairs
persistent_features = filter_features_by_persistence(minIndices, maxIndices, persistence, threshold);
filteredMinima = [persistent_features(:,1) ; globalMinIndex];
filteredMaxima = persistent_features(:,2);

% Plot only features with persistence > threshold.
figure;
plot(data, '-k', 'LineWidth', 2);
title(strcat('extrema with persistence > ', num2str(threshold)));
hold on;

% Add a scatter plot for the filtered minima
plot(filteredMinima, data(filteredMinima), 'o', 'MarkerSize', 9, 'MarkerFaceColor', [0.3 0.3 1], 'MarkerEdgeColor', [0 0 1]);

% Add a scatter plot for the filtered maxima
plot(filteredMaxima, data(filteredMaxima), 'o', 'MarkerSize', 9, 'MarkerFaceColor', [1 0.2 0.2], 'MarkerEdgeColor', [1 0 0]);

hold off;

Example for Smoothing Data

The Matlab scripts for smoothing the function are located in the reconstruct1d folder. Smoothing is only possible from within Matlab.

% Add Reconstruct1D folder to Matlab's path
addpath('..');

% Setup Persistence1D and MOSEK
setup_persistence1d();
turn_on_mosek();

% Load the data set
load '..\datasets\test_data.mat';

% Choose smoothness for the reconstructed function.
% 'biharmonic' smoothness guarantees that the reconstructed function is C1 smooth
% 'triharmonic' smoothness guarantees that the reconstructed function is C2 smooth
smoothness = 'biharmonic';

% Choose a threshold for persistence features
threshold = 0.2;

% The data term weight affects how closely the reconstructed function
% adheres to the data.
data_weight = 0.0000001;

x = reconstruct1d(data, threshold, smoothness, data_weight);
plot_reconstructed_data(data, x, smoothness, threshold, data_weight);
turn_off_mosek();

Documentation

The download package comes with extensive documentation and examples.

Download

Persistence1D has been written by Yeara Kozlov and Tino Weinkauf, Max Planck Institute for Informatics, Saarbrücken, Germany. You may use it as you wish, it is in the public domain. If you find it useful, it would be nice to hear from you. Just drop us a line.

Contents