Thesis/canny_shock_finder.py

"""
Canny_shock_finder.py

A little shock finding code which I [Chris] wrote using a cut down one-dimensional Canny edge detector.
Seems to work quite well with the appropriate selection of input parameters.

The original paper is:

Canny, J. (1986). A computational approach to edge detection. 
IEEE Transactions on pattern analysis and machine intelligence, (6), 679-698.

Chris James (c.james4@uq.edu.au) - 04/07/17

Updated in 2024 by Cal Wing to use Python 3 
    This refactor *REMOVED* some features, namely graphing
"""

__version__ = "0.1.0"
__date__ = "15/10/2024"

def canny_shock_finder(time_list, pressure_list, sigma = 4, derivative_threshold = 0.001, auto_derivative = False, post_suppression_threshold = 0.05,
                       post_shock_pressure = None, start_time = None, find_second_peak = None, plot = True, plot_scale_factor = None,
                       amount_of_values_before_and_after = 100, time_unit = 's', y_label = 'Pressure (kPa)', plot_title = 'canny shock finding result',
                       return_processing_lists = False, calculate_automatic_derivative_threshold = False,
                       check_canny_runtime = False, plot_time_unit = 's', plot_time_scale = 1.0, print_diag_info: bool = False):
    """
    This is a basic function with uses a cut down one-dimensional Canny edge detector to find shock arrival.
    It seems to work quite well with the appropriate selection of input parameters 
    (sigma, derivative_threshold, post_suppression_threshold).
    
    The original paper is:

    Canny, J. (1986). A computational approach to edge detection. 
    IEEE Transactions on pattern analysis and machine intelligence, (6), 679-698.
    
    Canny edge detection is a multi-stage algorithm generally used to find edges for image detection.
    It is a widely used and discussed, and it should not be hard for a user to find adequate 
    information about it. (The original paper by Canny is also quite a well presented and easy to read
    resource as well.)
    
    In two dimensions (i.e. for an image) Canny edge detection is a bit more complicated, but
    a simplified one dimensioanl version has been implemented here. The first step is to filter
    the image by convolving it with the derivative of a Gaussian (this is done by using the
    gaussian_filter1d functions from scipy.ndimage.filters with an order of 1), and then find
    gradients in the image in the x and y directions with generally the Sobel method. As our 
    implementation here is one-dimensional, this second step has been performed by convolving the
    original signal with the second derivative of a Gaussian to get the derivative of the original 
    derivative of a Gaussian convolution (this is done using the gaussian_filter1d functions from 
    scipy.ndimage.filters with an order of 2). The standard deviation (sigma) of these Gaussians is 
    defined by the user as an input. 
    At this point, there are two arrays, which are the original data convolved with both the 
    first and second derivatives of a Gaussian. The first array has maximums where step changes occur, 
    and the other crosses zero where step changes occur. Here the derivative_threshold input is used by 
    the code to inform it what to call zero in the second order Gaussian, as due to floating point uncertainty the 
    values don't come out as zero, but very close to it. This value is somewhat dependent on the magnitude of the 
    original signal, and it would be expected to scale with the magnitude of the original signal.
    The next step in the Canny edge detection process is "Non-maximum suppression" where any value 
    which does not have a different gradient on either side of it is set to zero, "thinning" the edges 
    to generally a single point.
    (If the code ends up with two edges next to each other due to issues with the signal, a midpoint
    with a +- uncertainty is taken). 
    Finally, a low and high threshold is used to define less important edges and what is certaintly an
    edge. Here, as we want one strong edge for the shock arrival, any edges below the lower threshold 
    (the post_suppression_threshold) are set to 0, and a high threshold is not used. Similar to the 
    derivative_threshold discussed above, the post_suppression_threshold would also be expected to scale
    with the magnitude of the input signal.
    
    As the selection of input parameters do scale with the magnitude of the signal, they must be modified
    as required, but it is hoped that they can hopefully be setup for one experimental campaign,
    and then used for all analysis after that for each signal. 
    
    The default setting is for the code to find only one peak, but it can also find a second peak 
    for reflected shock tunnel usage if required.
    
    Returns the time and uncertainty for both the first and second peaks in the data,
    which should be the first and second shock arrival times if the code is configured correctly
    for the data. If the second peak is not asked for, its value are returned as None
    
    :param time_list: List or array of time values. Is turned into an array by the code either way.
    :param pressure_list: List or array of pressure (or other quantity) values. Is turned into an array by the code either way.
    :param sigma: standard deviation to be used with the Gaussian filters. Defaults to 4 but should be changed
        to suit the data.
    :param derivative_threshold: threshold below which the non-maximum supression part of the calculation
        will take values to be zero. Used to counter numerical errors causing the second order Gaussian values
        where it crosses zero to not be exactly zero. Defaults to 0.001 but should be changed to suit the data.
    :param post_suppression_threshold: threshold applied to the data after the non-maximum supression
        has been performed. Is basically tuned to ensure that any values left in the noise before
        shock arrival are removed in favour of the larger peak when the shock arrives. Defaults to 
        0.05 but should be modified as required until the first peak found by the code is the first shock arrival time.
    :param post_shock_pressure: A post shock pressure estimate can be entered instead of the post_supression_threshold.
        It will then be used to calculate a new rough post_suppression threshold based on the post-shock pressure value
        and the sigma value. Input should be in the same units as the data
    :param start_time: start_time can be used to make the code start at a certain time. Any value below that will be ignored. 
        Defaults to None (i.e. not used)
    :param find_second_peak: Input flag for reflected shock tunnels where there is expected to be a 
        second shock arrival. Defaults to None and will then be set to the first value in the time_list.
    :param plot: flag to plot the result. Defaults to True.
    :param plot_scale_factor: used to scale the results so they have a similar magnitude to the original data.
        Defaults to None, which means that no scaling will be used.
    :param amount_of_values_before_and_after: amount of values before and after the found shock arrival time
        to use for calculating the y limits for the plot. Defaults to 100.    
    :param time_unit: time unit for the plotting. Defaults to 's (seconds)'.
    :param pressure_unit: time unit for the y axis Defaults to 'kPa' (pressure in kilopascals).
    :param plot_title: string title for the plot. Defaults to 'canny shock finding result'
    :param return_processing_lists: returns a dictionary with the tine and pressure lists, as well as the gaussian filtered results,
        and the non-maximum suppression and final result lists. Useful for making plots of how the process has worked.
        Defaults to False.
    :param calculate_automatic_derivative_threshold: Makes the program calculate an automatic derivative threshold basically using
        an emprirical correlation which I worked out based on how the required derivative threshold changes with the sigma value.
    :param check_canny_runtime: flag to make the program time how long the Canny edge detection procedure
        actually takes. Mainly done for testing purposes.
    :param plot_time_unit: Allows a different time unit to be specified for the results plot, so that the
        original calculation can be performed in the original units, but then plotted in another. Relies on
        the user using the correct plot_time_scale factor to connect the two time units. Defaults to 's' (seconds).
    :param plot_time_scale: See plot_time_unit above. These two inputs allow the results plot to be in a different
        time unit to the original data. The user must specify this scale to connect the original and plotted
        time units. Defaults to 1.0 (so no change if time units are seconds).
    """