Demo entry 3870570

Thresholding functions

   

Submitted by Chris Simpson on Mar 04, 2016 at 13:19
Language: Python 3. Code size: 4.2 kB.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#  Copyright (c) 2016 Chris Simpson, University of Manchester

#  Permission is hereby granted, free of charge, to any person obtaining a
#  copy of this software and associated documentation files (the “Software”),
#  to deal in the Software without restriction, including without limitation
#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
#  and/or sell copies of the Software, and to permit persons to whom the
#  Software is furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included
#  in all copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#  DEALINGS IN THE SOFTWARE.

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal

def plot_binned(fname, window = 101, total_peaks = 3):
    """
    Plot the binned data and visualise the peaks.
    Allows for the optimisation of peak for peak finding
    algorithm.

    # fname:       File name (inc path if not local)
    # window:      Window to use for peak finding algorithm
    # total_peaks:
    """
    data = np.loadtxt(fname, delimiter=',')
    peak_locs = signal.argrelmax(data[:,1], order=window)[0]
    plt.figure(figsize = (9, 6))
    plt.semilogy(data[:,0], data[:,1])
    plt.ylabel('Counts')
    plt.xlabel('Greyscale value')
    for i in range(total_peaks):
        plt.semilogy(data[peak_locs[i], 0], data[peak_locs[i], 1], 'ro')


def thresholding(fname, window, method = 'centre', peak = 1,
                 total_peaks = 3, plot = True):
    """
    Simple thresholding function for binned imaging data. Returns
    number of thresholded pixels. Specify the peak of interest
    and total peaks. Works best with clearly defined peaks.

    # fname:       File name (inc path if not local)
    # window:      Window to use for peak finding algorithm
    # method:      Centre of minima thresholding
    # peak:        Peak of interest (1 based indexing)
    # total_peaks: Total peaks (1 based indexing)
    # plot:        Plot data, peaks and thresholding limits
    """
    data = np.loadtxt(fname, delimiter = ',')
    peaks = signal.argrelmax(data[:,1], order=window)[0]

    if method.lower() == 'minima':
        if peak == 1:
            th1 = 0
        else:
            th1 = (signal.argrelmin(data[peaks[peak - 2]:peaks[peak - 1], 1],
                     order=window)[0] + peaks[peak - 2])[0]
        th2 = (signal.argrelmin(data[peaks[peak - 1]:peaks[peak], 1],
                 order=window)[0] + peaks[peak - 1])[0]
    else:
        if peak == 1:
            th1 = 0
        else:
            th1 = (peaks[peak -1] - peaks[peak - 2])//2 + peaks[peak - 2]
        th2 = (peaks[peak] - peaks[peak - 1])//2 + peaks[peak - 1]

    th1 = 0 if peak == 1 else th1
    th2 = data[:, 0].size -1 if total_peaks == peak else th2

    if plot:
        y_lim = np.min(data[:, 1]), np.max(data[:, 1])
        plt.figure(figsize = (9, 6))
        plt.semilogy(data[:,0], data[:,1])
        for i in range(total_peaks):
            plt.semilogy(data[peaks[i], 0], data[peaks[i], 1], 'ro')
        plt.semilogy([data[th1, 0], data[th1, 0]], [y_lim[0], y_lim[1]], 'k--')
        plt.semilogy([data[th2, 0], data[th2, 0]], [y_lim[0], y_lim[1]], 'k--')
        plt.ylabel('Counts')
        plt.xlabel('Greyscale value')

    thresholded_pix = int(np.sum(data[th1:th2, 1]))

    return thresholded_pix

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).