Source code for maad.features.shape

#!/usr/bin/env python3 
# -*- coding: utf-8 -*- 
""" 
Ensemble of functions to compute acoustic descriptors from 2D spectrograms.
""" 
#
# Authors:  Juan Sebastian ULLOA <lisofomia@gmail.com>
#           Sylvain HAUPERT <sylvain.haupert@mnhn.fr>        
#
# License: New BSD License

#%%
#***************************************************************************
# -------------------       Load modules         ---------------------------
#***************************************************************************
# Import external modules

from __future__ import print_function 
import numpy as np 
import pandas as pd 
from scipy import ndimage as ndi
import itertools as it 
import matplotlib.pyplot as plt 
from skimage import transform, measure 
from scipy import ndimage 

# Import internal modules
from maad.util import format_features, overlay_rois, plot_shape, overlay_centroid
from maad.rois import rois_to_imblobs
from maad.sound import spectrogram 
 
#%%
#============================================================================ 
# PRIVATE FUNCTIONS 
#============================================================================ 
def _gabor_kernel_nodc(frequency, theta=0, bandwidth=1, gamma=1, 
                      n_stds=3, offset=0): 
    """ 
    Computes complex 2D Gabor filter kernel with no DC offset. 
         
    Gabor kernel is a Gaussian kernel modulated by a complex harmonic function. 
    Harmonic function consists of an imaginary sine function and a real 
    cosine function. Spatial frequency is inversely proportional to the 
    wavelength of the harmonic and to the standard deviation of a Gaussian 
    kernel. The bandwidth is also inversely proportional to the standard 
    deviation. 
     
    Parameters 
    ---------- 
    frequency : float 
        Spatial frequency of the harmonic function. Specified in pixels. 
    theta : float, optional 
        Orientation in radians. If 0, the harmonic is in the x-direction. 
    bandwidth : float, optional 
        The bandwidth captured by the filter. For fixed bandwidth, `sigma_x` 
        and `sigma_y` will decrease with increasing frequency. This value is 
        ignored if `sigma_x` and `sigma_y` are set by the user. 
    gamma : float, optional 
        gamma changes the aspect ratio (ellipsoidal) of the gabor filter.  
        By default, gamma=1 which means no aspect ratio (circle) 
         
        if gamma>1, the filter is larger (x-dir) 
         
        if gamma<1, the filter is higher (y-dir) 
         
        This value is ignored if `sigma_x` and `sigma_y` are set by the user. 
    sigma_x, sigma_y : float, optional 
        Standard deviation in x- and y-directions. These directions apply to 
        the kernel before rotation. If `theta = pi/2`, then the kernel is 
        rotated 90 degrees so that `sigma_x` controls the vertical direction.    
    n_stds : scalar, optional 
        The linear size of the kernel is n_stds (3 by default) standard 
        deviations 
    offset : float, optional 
        Phase offset of harmonic function in radians. 
     
    Returns 
    ------- 
    g_nodc : complex 2d array 
        A complex gabor kernel with no DC offset 
     
    Notes 
    ----- 
    This function is a modification of the gabor_kernel function of scikit-image. 
     
    References 
    ---------- 
    .. [1] http://en.wikipedia.org/wiki/Gabor_filter 
    .. [2] http://mplab.ucsd.edu/tutorials/gabor.pdf 
    .. [3] http://www.cs.rug.nl/~imaging/simplecell.html 
     
    """ 
     
     # set gaussian parameters 
    b = bandwidth 
    sigma_pref = 1.0 / np.pi * np.sqrt(np.log(2) / 2.0) * (2.0 ** b + 1) / (2.0 ** b - 1) 
    sigma_y = sigma_pref / frequency 
    sigma_x = sigma_y / gamma 
    # meshgrid 
    x0 = np.ceil(max(np.abs(n_stds * sigma_x * np.cos(theta)), 
                     np.abs(n_stds * sigma_y * np.sin(theta)), 1)) 
    y0 = np.ceil(max(np.abs(n_stds * sigma_y * np.cos(theta)), 
                     np.abs(n_stds * sigma_x * np.sin(theta)), 1)) 
    y, x = np.mgrid[-y0:y0 + 1, -x0:x0 + 1] 
    # rotation matrix 
    rotx = x * np.cos(theta) + y * np.sin(theta) 
    roty = -x * np.sin(theta) + y * np.cos(theta) 
    # combine gambor and  
    g = np.zeros(y.shape, dtype=complex) 
    g[:] = np.exp(-0.5 * (rotx ** 2 / sigma_x ** 2 + roty ** 2 / sigma_y ** 2)) 
    g /= 2 * np.pi * sigma_x * sigma_y # gaussian envelope 
    oscil = np.exp(1j * (2 * np.pi * frequency * rotx + offset)) # harmonic / oscilatory function 
    g_dc = g*oscil 
    # remove dc component by subtracting the envelope weighted by K 
    K = np.sum(g_dc)/np.sum(g) 
    g_nodc = g_dc - K*g 
  
    return g_nodc 
 
#%%
def _plot_filter_bank(kernels, frequency, ntheta): 
    """ 
    Display filter bank 
     
    Parameters 
    ---------- 
    kernels: list 
        List of kernels from filter_bank_2d_nodc() 
         
    frequency: 1d ndarray of scalars 
        Spatial frequencies used to built the Gabor filters. Values should be 
        in [0;1] 
     
    ntheta: int 
        Number of angular steps between 0° to 90° 
     
    Returns 
    ------- 
    fig : Figure 
        The Figure instance  
    ax : Axis 
        The Axis instance 
    """ 
 
    ntheta = ntheta 
    nfreq = len(frequency) 
     
    # get maximum size 
    aux = list() 
    for kernel in kernels: 
        aux.append(max(kernel[0].shape)) 
    max_size = np.max(aux) 
     
    # plot kernels 
    fig, ax = plt.subplots(nfreq, ntheta) 
    ax = ax.transpose() 
    ax = ax.ravel() 
    for idx, k in enumerate(kernels): 
        kernel = np.real(k[0]) 
        ki, kj = kernel.shape 
        ci, cj = int(max_size/2 - ki/2), int(max_size/2 - kj/2) 
        canvas = np.zeros((max_size,max_size)) 
        canvas[ci:ci+ki,cj:cj+kj] = canvas[ci:ci+ki,cj:cj+kj] + kernel 
        ax[idx].imshow(canvas) 
        ax[idx].axis('off') 
    plt.tight_layout() 
    plt.subplots_adjust(wspace=0.02, hspace=0.02) 
    return ax, fig 
 
#%% 
def _plot_filter_results(im_ref, im_list, kernels, params, m, n): 
    """ 
    Display the resulting spectrograms after filtering with gabor filters 
     
    Parameters 
    ---------- 
    im_ref : 2D array 
        Reference image 
    im_list : list 
        List of filtered images 
    kernels: list 
        List of kernels from filter_bank_2d_nodc() 
    m: int 
        number of columns 
    n: int 
        number of rows 
         
    Returns 
    ------- 
    fig : Figure 
        The Figure instance  
    ax : Axis 
        The Axis instance 
    """     
         
    ncols = m  
    nrows = n 
     
    fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(15, 5)) 
    plt.gray() 
    fig.suptitle('Image responses for Gabor filter kernels', fontsize=12) 
    axes[0][0].axis('off') 
     
    # Plot original images 
    axes[0][1].imshow(im_ref, origin='lower') 
    axes[0][1].set_title('spectrogram', fontsize=9) 
    axes[0][1].axis('off') 
    plt.tight_layout 
     
    params_label = [] 
    for param in params: 
        params_label.append('theta=%d,\nf=%.2f' % (param[1] * 180 / np.pi, param[0])) 
     
    ii = 0 
    for ax_row in axes[1:]:  
        plotGabor = True 
        for ax in ax_row: 
            if plotGabor == True: 
                # Plot Gabor kernel 
                print(params_label[ii]) 
                ax.imshow(np.real(kernels[ii]), interpolation='nearest') 
                ax.set_ylabel(params_label[ii], fontsize=7) 
                ax.set_xticks([]) 
                ax.set_yticks([]) 
                plotGabor = False 
            else: 
                im_filtered = im_list[ii] 
                ax.imshow(im_filtered, origin='lower') 
                ax.axis('off') 
                plotGabor = True 
                ii=ii+1      
         
    plt.show() 
    return ax, fig 
 
#%%
def _filter_mag(im, kernel): 
    """ 
    Normalizes the image and computes im and real part of filter response using  
    the complex kernel and the modulus operation 
     
    Parameters 
    ---------- 
    im: 2D array 
        Input image to process  
    kernel: 2D array 
        Complex kernel or filter 
             
    Returns 
    ------- 
    im_out: Modulus operand on filtered image 
 
    """     
     
    im = (im - im.mean()) / im.std() 
    im_out = np.sqrt(ndi.convolve(im, np.real(kernel), mode='reflect')**2 + 
                   ndi.convolve(im, np.imag(kernel), mode='reflect')**2) 
    return im_out 
 
#%%
def _params_to_df(params_filter_bank, npyr): 
    """ 
    Organises parameters used in multiresolution analysis into a dataframe 
     
    Parameters 
    ---------- 
    params_filter_bank: numpy array 
        Output of parameters from function filter_bank_2d_nodc 
    npyr: int 
        Number of pyramids used in multiresolution analysis 
             
    Returns 
    ------- 
    params_multires: pandas DataFrame 
        Ordered parameters for each shape feature used 
    """     
 
    params = np.asarray(params_filter_bank) 
    orient = params[:,0] * 180 / np.pi 
    orient = orient.tolist() * npyr 
    pyr_level = np.sort(np.arange(npyr).tolist()*len(params))+1 
    freq = params[:,1].tolist() * npyr 
    nparams = len(params) * npyr 
    params_multires = np.zeros(nparams, dtype={'names':('theta', 'freq', 'pyr_level','scale'), 
                                               'formats':('f8', 'f8', 'f8','f8')}) 
    params_multires['theta'] = orient 
    params_multires['freq'] = freq 
    params_multires['scale'] = 1 / np.asarray(freq) 
    params_multires['pyr_level'] = pyr_level 
    params_multires = pd.DataFrame(params_multires) 
 
    return params_multires 
 
#%%
#============================================================================ 
# Public functions 
#============================================================================ 
[docs] def filter_multires(Sxx, kernels, npyr=4, rescale=True): """ Computes 2D wavelet coefficients at multiple scales using Gaussian pyramid transformation to downscale the input spectrogram. Parameters ---------- Sxx: list of 2D arrays List of input spectrograms to filter kernels: list of 2D arrays List of 2D kernels or filters npyr: int, optional Number of pyramids to compute. Default is 4. rescale: boolean, optional Indicates if the reduced images should be rescaled. Default is True. Returns ------- Sxx_out: list of 2D arrays List of spectrograms filtered by each 2D kernel Examples -------- >>> import maad >>> from maad.sound import load, spectrogram >>> from maad.features import filter_bank_2d_nodc, filter_multires >>> from maad import util >>> s, fs = load('../data/spinetail.wav') >>> Sxx, dt, df, ext = spectrogram(s, fs) >>> Sxx_db = util.power2dB(Sxx, db_range=80) + 80 >>> ax, fig = util.plot2d(Sxx_db, **{'extent':ext}) >>> params, kernels = filter_bank_2d_nodc(frequency=(0.5, 0.25), ntheta=2,gamma=2) >>> Sxx_out = filter_multires(Sxx, kernels, npyr=2) Plot one of the resulting spectrograms. >>> ax, fig = util.plot2d(Sxx_out[5], **{'extent':ext}) """ # Downscale image using gaussian pyramid if npyr<2: print('Warning: npyr should be int and larger than 2 for multiresolution') im_pyr = tuple(transform.pyramid_gaussian(Sxx, downscale=2, max_layer=1)) else: im_pyr = tuple(transform.pyramid_gaussian(Sxx, downscale=2, max_layer=npyr-1)) # filter 2d array at multiple resolutions using gabor kernels im_filt=[] for im in im_pyr: # for each pyramid for kernel, param in kernels: # for each kernel im_filt.append(_filter_mag(im, kernel)) # magnitude response of filter # Rescale image using gaussian pyramid if rescale: dims_raw = Sxx.shape Sxx_out=[] for im in im_filt: ratio = np.array(dims_raw)/np.array(im.shape) if ratio[0] > 1: im = transform.rescale(im, scale = ratio, mode='reflect', anti_aliasing=True) else: pass Sxx_out.append(im) else: pass return Sxx_out
#%%
[docs] def filter_bank_2d_nodc(frequency, ntheta, bandwidth=1, gamma=2, display=False, **kwargs): r""" Build an ensemble of complex 2D Gabor filters with no DC offset. A Gabor kernel is a Gaussian kernel modulated by a complex harmonic function. Harmonic function consists of an imaginary sine function and a real cosine function. Spatial frequency is inversely proportional to the wavelength of the harmonic and to the standard deviation of a Gaussian kernel. The bandwidth is also inversely proportional to the standard deviation [1]_. These filters have been used to classify textures in images [4]_ and can be applied to characterize spectrograms to classify sounds [5]_. Parameters ---------- frequency: 1d ndarray of scalars Spatial frequencies used to built the Gabor filters. Values should be in the range [0;1] ntheta: int Number of angular steps between 0° to 90° bandwidth: scalar, optional, default is 1 Spatial-frequency bandwidth of the filter. This parameters affect the resolution of the filters. Lower bandwidth result in sharper in filters with more details. gamma: scalar, optional, default is 1 Gaussian window that modulates the continuous sine. For ``gamma = 1``, the result is the same gaussian window in x and y direction (circle). For ``gamma <1``, the result is an increased elongation of the filter size in the y direction (elipsoid). For ``gamma >1``, the result is a reduction of the filter size in the y direction (elipsoid). display: bool Display a visualization of the filter bank. Default is False. Returns ------- params: 2d structured array Parameters used to calculate 2D gabor kernels. Params array has 4 fields (theta, freq, bandwidth, gamma) This can be useful to interpret the result of the filtering process. kernels: 2d ndarray of complex values Complex Gabor kernels References ---------- .. [1] http://en.wikipedia.org/wiki/Gabor_filter .. [2] http://mplab.ucsd.edu/tutorials/gabor.pdf .. [3] http://www.cs.rug.nl/~imaging/simplecell.html .. [4] Sifre, L., & Mallat, S. (2013). Rotation, scaling and deformation invariant scattering for texture discrimination. Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference On, 1233–1240. `DOI: 10.1109/CVPR.2013.163 <https://doi.org/10.1109/CVPR.2013.163>`_ .. [5] Ulloa, J. S., Aubin, T., Llusia, D., Bouveyron, C., & Sueur, J. (2018). Estimating animal acoustic diversity in tropical environments using unsupervised multiresolution analysis. Ecological Indicators, 90, 346–355. `DOI: 10.1016/j.ecolind.2018.03.026 <https://doi.org/10.1016/j.ecolind.2018.03.026>`_ See also -------- filter_multires, opt_shape_presets, shape_features Examples -------- >>> import maad Build a filter bank using predefined parameters with the function maad.features.opt_shape_presets. >>> from maad.features import filter_bank_2d_nodc, opt_shape_presets >>> opt = opt_shape_presets(resolution='med') >>> params, kernels = filter_bank_2d_nodc(opt['frequency'], opt['ntheta'], opt['bandwidth'], opt['gamma'], display=True) Alternatively, custom parameters can be provided to define the filter bank. >>> from maad.features import filter_bank_2d_nodc >>> params, kernels = filter_bank_2d_nodc(frequency=(0.7, 0.5, 0.35, 0.25), ntheta=4, gamma=2, display=True) """ theta = np.arange(ntheta) theta = theta / ntheta * np.pi params=[i for i in it.product(theta,frequency)] kernels = [] for param in params: kernel = _gabor_kernel_nodc(frequency=param[1], theta=param[0], bandwidth=bandwidth, gamma=gamma, offset=0, n_stds=3) kernels.append((kernel, param)) if display: _, fig = _plot_filter_bank(kernels, frequency, ntheta) return params, kernels
#%%
[docs] def opt_shape_presets(resolution, opt_shape=None): """ Set parameters for multiresolution analysis using presets or custom parameters. Parameters ---------- resolution: str Select resolution of analysis using presets. Supported presets are: 'low', 'med', and 'high'. Select 'custom' to select user-defined parameters using a dictionary. opt_shape: dict Key and values for shape settings. Valid keys are: 'ntheta', 'bandwidth', 'frequency', 'gamma', 'npyr' Returns ------- opt_shape: dict A valid dictionary with shape settings See also -------- filter_bank_2d_nodc, filter_multires, shape_features Examples -------- >>> import maad Get parameters to analyse at low, med and high shape resolutions. >>> from maad.features import opt_shape_presets >>> opt_shape_presets('low') # doctest: +NORMALIZE_WHITESPACE {'ntheta': 2, 'bandwidth': 0.8, 'frequency': (0.35, 0.5), 'gamma': 2, 'npyr': 4} >>> opt_shape_presets('med') # doctest: +NORMALIZE_WHITESPACE {'ntheta': 4, 'bandwidth': 0.8, 'frequency': (0.35, 0.5), 'gamma': 2, 'npyr': 6} >>> opt_shape_presets('high') # doctest: +NORMALIZE_WHITESPACE {'ntheta': 8, 'bandwidth': 0.8, 'frequency': (0.35, 0.5), 'gamma': 2, 'npyr': 6} """ # Factory presets opt_shape_low = dict(ntheta=2, bandwidth=0.8, frequency=(0.35, 0.5), gamma=2, npyr = 4) opt_shape_med = dict(ntheta=4, bandwidth=0.8, frequency=(0.35, 0.5), gamma=2, npyr = 6) opt_shape_high = dict(ntheta=8, bandwidth=0.8, frequency=(0.35, 0.5), gamma=2, npyr = 6) if resolution == 'low': opt_shape = opt_shape_low elif resolution == 'med': opt_shape = opt_shape_med elif resolution == 'high': opt_shape = opt_shape_high elif resolution == 'custom': if opt_shape is not None: # check valid values on opt_shape if all (opt in opt_shape for opt in ('ntheta', 'bandwidth', 'frequency', 'gamma', 'npyr')): pass else: print('Warning: opt_shape must have all keys-values pairs:') print('ntheta, bandwidth, frequency, gamma, npyr') print('Setting resolution to low') opt_shape = opt_shape_low else: print('Warning: if resolution is set to custom, a valid opt_shape dictionnary should be provided.') print('Setting resolution to low') opt_shape = opt_shape_low else: raise TypeError("Resolution should be: 'low', 'med' or 'high'") opt_shape = opt_shape_low return opt_shape
#%%
[docs] def shape_features(Sxx, resolution='low', rois=None): """ Computes time-frequency shape coefficients at multiple resolutions using 2D Gabor filters. Parameters ---------- Sxx: 2D array Input image to process resolution: str or dict Specify resolution of shape descriptors. Can be: 'low', 'med', 'high'. Default is 'low'. Alternatively, custom resolution can be provided using a dictionary with options to define the filter bank. Valid keys are: ntheta, bandwidth, frequency, gamma, npyr. rois: pandas DataFrame Regions of interest where descriptors will be computed. Array must have a valid input format with column names: min_t min_f, max_t, max_f. Use format_features(rois,tn,fn) before using shape_features to be sure that the format of the rois DataFrame is correct. Returns ------- shape: pandas DataFrame Shape coefficient for each region of interest params: 2D numpy structured array Corresponding parameters of the 2D fileters used to compute the shape coefficients. Params has 4 fields (theta, freq, pyr_level, scale) See Also -------- maad.features.opt_shape_presets, maad.features.filter_bank_2d_nodc, maad.util.plot_shape References ---------- .. [1] Ulloa, J. S., Aubin, T., Llusia, D., Bouveyron, C., & Sueur, J. (2018). Estimating animal acoustic diversity in tropical environments using unsupervised multiresolution analysis. Ecological Indicators, 90, 346–355. `DOI: 10.1016/j.ecolind.2018.03.026 <https://doi.org/10.1016/j.ecolind.2018.03.026>`_ .. [2] Sifre, L., & Mallat, S. (2013). Rotation, scaling and deformation invariant scattering for texture discrimination. Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference On, 1233–1240. `DOI: 10.1109/CVPR.2013.163 <https://doi.org/10.1109/CVPR.2013.163>`_ .. [3] Mallat, S. (2008). A Wavelet Tour of Signal Processing: The Sparse Way. Academic Press. `DOI: 10.1016/B978-0-12-374370-1.X0001-8 <https://doi.org/10.1016/B978-0-12-374370-1.X0001-8>`_ Examples -------- Get shape features from the whole power spectrogram >>> from maad.sound import load, spectrogram >>> from maad.features import shape_features >>> from maad.util import format_features, power2dB, plot_shape >>> s, fs = load('../data/spinetail.wav') >>> Sxx, tn, fn, ext = spectrogram(s, fs, db_range=100, display=True) >>> Sxx_db = power2dB(Sxx, db_range=100) >>> shape, params = shape_features(Sxx_db, resolution='med') >>> ax = plot_shape(shape.mean(), params) Or get shape features from specific regions of interest >>> from maad.util import read_audacity_annot >>> rois = read_audacity_annot('../data/spinetail.txt') >>> rois = format_features(rois, tn, fn) >>> shape, params = shape_features(Sxx_db, resolution='med', rois=rois) The shape coefficient can be visualized using ``plot_shape``. Note that these coefficient allow to disciminate between the two vocalization types found in the audio file, the spinetail song and the background calls. >>> shape_crer, shape_sp = shape.groupby('label') >>> ax1 = plot_shape(shape_crer[1].filter(regex=r'^shp_').mean(), params) >>> ax1 = ax1.set_title('Spinetail song') >>> ax2 = plot_shape(shape_sp[1].filter(regex=r'^shp_').mean(), params) >>> ax2 = ax2.set_title('Background call') # doctest: +NORMALIZE_WHITESPACE """ # Check input data and unpack settings if type(Sxx) is not np.ndarray and len(Sxx.shape) != 2: raise TypeError('Sxx must be a numpy 2D array') if type(resolution) is str: opt_shape = opt_shape_presets(resolution) elif type(resolution) is dict: opt_shape = opt_shape_presets(resolution='custom', opt_shape=resolution) else: raise TypeError('Resolution must be string or a dictionary. See function documentation.') npyr = opt_shape['npyr'] # build filterbank params, kernels = filter_bank_2d_nodc(ntheta=opt_shape['ntheta'], bandwidth=opt_shape['bandwidth'], frequency=opt_shape['frequency'], gamma=opt_shape['gamma']) # filter spectrogram im_rs = filter_multires(Sxx, kernels, npyr, rescale=True) # If rois are provided get mean intensity for each ROI, # else compute mean intensity for the whole spectrogram if rois is None: shape = [] for im in im_rs: shape.append(np.mean(im)) shape = [shape] # for dataframe formating below else: shape = np.zeros(shape=(len(rois),len(im_rs))) im_rs = np.array(im_rs) index_row = 0 for _, row in rois.iterrows(): im_rs_blobs = im_rs[:, int(row.min_y):int(row.max_y)+1, int(row.min_x):int(row.max_x)+1] roi_mean = np.mean(im_rs_blobs, axis=(1,2)) shape[index_row,:] = roi_mean index_row = index_row + 1 # organise parameters params_multires = _params_to_df(params, npyr) # format shape into dataframe cols=['shp_' + str(idx).zfill(3) for idx in range(1,len(shape[0])+1)] shape = pd.DataFrame(data=np.asarray(shape),columns=cols) if rois is not None: shape.index = rois.index shape = rois.join(shape) return shape, params_multires
#%%
[docs] def shape_features_raw(im, resolution='low', opt_shape=None): """ Computes raw shape of 2D signal (image or spectrogram) at multiple resolutions using 2D Gabor filters. Contrary to ``shape_features``, this function delivers the raw response of the spectrogram to each filter of the filter bank. This function can be used to have access to spectrograms filtered by the filterbank. Parameters ---------- Sxx: 2D array Spectrogram to be analysed resolution: Resolution of analysis, i.e. number of filters used. Three presets are provided, 'low', 'mid' and 'high', which control the number of filters. opt_shape: dictionary (optional) options for the filter bank (kbank_opt) and the number of scales (npyr) Returns ------- shape_raw: list Raw shape response of spectrogram to every filter of the filter bank. Every element of the list is the response to one of the filters. Detail of each filter are available in the ``params`` DataFrame returned. params: pandas DataFrame Corresponding parameters of the 2D fileters used to compute the shape coefficient. Params has 4 fields (theta, freq, pyr_level, scale) See also -------- maad.features.shape_features, maad.features.opt_shape_presets, maad.features.filter_bank_2d_nodc, maad.util.plot_shape Examples -------- >>> import maad Scatter the spectrogram acoustic components using 2D Gabor filters. >>> from maad.sound import load, spectrogram >>> from maad.features import shape_features_raw >>> from maad.util import power2dB, plot2d >>> s, fs = load('../data/spinetail.wav') >>> Sxx, tn, fn, ext = spectrogram(s, fs, db_range=80, display=True) >>> Sxx_db = power2dB(Sxx, db_range=80) >>> shape_raw, params = shape_features_raw(Sxx_db, resolution='low') >>> fig, ax = plot2d(shape_raw[0], **{'extent':ext, 'title': 'High-frequency vertical component'}) >>> fig, ax = plot2d(shape_raw[13], **{'extent':ext, 'title': 'Low-frequency vertical components'}) >>> fig, ax = plot2d(shape_raw[2], **{'extent':ext, 'title': 'High-frequency horizontal components'}) >>> fig, ax = plot2d(shape_raw[15], **{'extent':ext, 'title': 'Low-frequency horizontal components'}) # doctest: +NORMALIZE_WHITESPACE """ # unpack settings opt_shape = opt_shape_presets(resolution, opt_shape) npyr = opt_shape['npyr'] # build filterbank params, kernels = filter_bank_2d_nodc(ntheta=opt_shape['ntheta'], bandwidth=opt_shape['bandwidth'], frequency=opt_shape['frequency'], gamma=opt_shape['gamma']) # filter images shape_raw = filter_multires(im, kernels, npyr, rescale=True) # organise parameters params_multires = _params_to_df(params, npyr) return shape_raw, params_multires
#%%
[docs] def centroid_features(Sxx, rois=None, im_rois=None): """ Computes intensity centroid of a spectrogram. If regions of interest ``rois`` are provided, the centroid is computed for each region. Parameters ---------- Sxx : 2D array Spectrogram rois: pandas DataFrame, default is None Regions of interest where descriptors will be computed. Array must have a valid input format with column names: ``min_t``, ``min_f``, ``max_t``, and ``max_f``. Use the function ``maad.util.format_features`` before using centroid_features to format of the ``rois`` DataFrame correctly. im_rois: 2d ndarray image with labels as values Returns ------- centroid: pandas DataFrame Centroid of each region of interest. See Also -------- maad.features.shape_features, maad.util.overlay_rois, maad.util.format_features Examples -------- >>> import maad Get centroid from the whole power spectrogram >>> from maad.sound import load, spectrogram >>> from maad.features import centroid_features >>> from maad.util import (power2dB, format_features, overlay_rois, plot2d, overlay_centroid) Load audio and compute spectrogram >>> s, fs = load('../data/spinetail.wav') >>> Sxx,tn,fn,ext = spectrogram(s, fs, db_range=80) >>> Sxx = power2dB(Sxx, db_range=80) Load annotations and plot >>> from maad.util import read_audacity_annot >>> rois = read_audacity_annot('../data/spinetail.txt') >>> rois = format_features(rois, tn, fn) >>> ax, fig = plot2d (Sxx, extent=ext) >>> ax, fig = overlay_rois(Sxx,rois, extent=ext, ax=ax, fig=fig) Compute the centroid of each rois, format to get results in the temporal and spectral domain and overlay the centroids. >>> centroid = centroid_features(Sxx, rois) >>> centroid = format_features(centroid, tn, fn) >>> ax, fig = overlay_centroid(Sxx,centroid, extent=ext, ax=ax, fig=fig) """ # Check input data if type(Sxx) is not np.ndarray and len(Sxx.shape) != 2: raise TypeError('Sxx must be an numpy 2D array') # check rois if rois is not None: if not(('min_t' and 'min_f' and 'max_t' and 'max_f') in rois): raise TypeError('Array must be a Pandas DataFrame with column names: min_t, min_f, max_t, max_f. Check example in documentation.') centroid=[] area = [] if rois is None: centroid = ndimage.center_of_mass(Sxx) centroid = pd.DataFrame(np.asarray(centroid)).T centroid.columns = ['centroid_y', 'centroid_x'] centroid['area_xy'] = Sxx.shape[0] * Sxx.shape[1] centroid['duration_x'] = Sxx.shape[1] centroid['bandwidth_y'] = Sxx.shape[0] else: if im_rois is not None : # real centroid and area rprops = measure.regionprops(im_rois, intensity_image=Sxx) centroid = [roi.weighted_centroid for roi in rprops] area = [roi.area for roi in rprops] else: # rectangular area (overestimation) area = (rois.max_y -rois.min_y) * (rois.max_x -rois.min_x) # centroid of rectangular roi for _, row in rois.iterrows() : row = pd.DataFrame(row).T im_blobs = rois_to_imblobs(np.zeros(Sxx.shape), row) rprops = measure.regionprops(im_blobs, intensity_image=Sxx) centroid.append(rprops.pop().weighted_centroid) centroid = pd.DataFrame(centroid, columns=['centroid_y', 'centroid_x'], index=rois.index) ##### duration in number of pixels centroid['duration_x'] = (rois.max_x -rois.min_x) ##### bandwidth in number of pixels centroid['bandwidth_y'] = (rois.max_y -rois.min_y) ##### area centroid['area_xy'] = area # concat rois and centroid dataframes centroid = rois.join(pd.DataFrame(centroid, index=rois.index)) return centroid
#%%
[docs] def all_shape_features(s, fs, rois, resolution='low', display=False, **kwargs): """ Computes shape and central frequency features from signal at specified time-frequency limits defined by regions of interest (ROIs). Parameters ---------- s: ndarray Singal to be analysed fs: int Sampling frequency of the signal rois: pandas DataFrame Time frequency limits for the analysis. Columns should have at least min_t, max_t, min_f, max_f. Can be computed with multiple detection methods, such as find_rois_cwt resolution: str or dict, default is 'med' Specify resolution of shape descriptors. Can be: 'low', 'med', 'high'. Default is 'low'. Alternatively, custom resolution can be provided using a dictionary with options to define the filter bank. Valid keys are: ntheta, bandwidth, frequency, gamma, npyr opt_spec: dictionnary, optional Options for the spectrogram with keys, window length 'nperseg', window overlap in percentage 'overlap' and frequency limits 'flims' which is a list of 2 scalars : minimum and maximum boundary frequency values in Hertz display: boolean, optional, default is False Flag. If display is True, plot results kwargs, optional. This parameter is used by plt.plot and savefig functions - savefilename : str, optional, default :'_spectro_overlaycentroid.png' Postfix of the figure filename - figsize : tuple of integers, optional, default: (4,10) width, height in inches. - title : string, optional, default : 'Spectrogram' title of the figure - xlabel : string, optional, default : 'Time [s]' label of the horizontal axis - ylabel : string, optional, default : 'Amplitude [AU]' label of the vertical axis - cmap : string or Colormap object, optional, default is 'gray' See https://matplotlib.org/examples/color/colormaps_reference.html in order to get all the existing colormaps examples: 'hsv', 'hot', 'bone', 'tab20c', 'jet', 'seismic', 'viridis'... - vmin, vmax : scalar, optional, default: None `vmin` and `vmax` are used in conjunction with norm to normalize luminance data. Note if you pass a `norm` instance, your settings for `vmin` and `vmax` will be ignored. - dpi : integer, optional, default is 96 Dot per inch. For printed version, choose high dpi (i.e. dpi=300) => slow For screen version, choose low dpi (i.e. dpi=96) => fast - format : string, optional, default is 'png' Format to save the figure ... and more, see matplotlib Returns ------- feature_rois: pandas Dataframe A dataframe with each column corresponding to a feature Examples -------- >>> import maad >>> from maad.util import read_audacity_annot >>> from maad.sound import load >>> from maad.features import all_shape_features >>> s, fs = load('../data/spinetail.wav') >>> rois = read_audacity_annot('../data/spinetail.txt') >>> features = all_shape_features(s, fs, rois, resolution='low', display=True) # doctest: +NORMALIZE_WHITESPACE number of rois : 18 >>> print(features.iloc[0]) # doctest: +NORMALIZE_WHITESPACE label SP min_t 0.101385 min_f 6441.064453 max_t 0.36752 max_f 12296.577148 min_y 150 min_x 8 max_y 286 max_x 31 shp_001 0.160656 shp_002 0.155573 shp_003 0.064787 shp_004 0.066215 shp_005 0.158675 shp_006 0.125203 shp_007 0.044987 shp_008 0.043188 shp_009 0.218268 shp_010 0.158725 shp_011 0.043553 shp_012 0.031206 shp_013 0.268222 shp_014 0.21718 shp_015 0.059112 shp_016 0.038394 centroid_y 220.0 centroid_x 19.0 duration_x 135.999274 bandwidth_y 0.0 area_xy 12512 centroid_f 9474.609375 centroid_t 0.2322 duration_t 0.267029 bandwidth_f 5857.03125 area_tf 1564.0 Name: 0, dtype: object >>> features.columns # doctest: +NORMALIZE_WHITESPACE Index(['label', 'min_t', 'min_f', 'max_t', 'max_f', 'min_y', 'min_x', 'max_y', 'max_x', 'shp_001', 'shp_002', 'shp_003', 'shp_004', 'shp_005', 'shp_006', 'shp_007', 'shp_008', 'shp_009', 'shp_010', 'shp_011', 'shp_012', 'shp_013', 'shp_014', 'shp_015', 'shp_016', 'centroid_y', 'centroid_x', 'duration_x', 'bandwidth_y', 'area_xy', 'centroid_f', 'centroid_t', 'duration_t', 'bandwidth_f', 'area_tf'], dtype='object') """ # Spectro computing N=1024 window = kwargs.pop('window','hann') nperseg = kwargs.pop('nperseg',N) overlap = kwargs.pop('overlap',N//2) mode = kwargs.pop('mode', 'psd') fcrop = kwargs.pop('fcrop',None) tcrop = kwargs.pop('tcrop',None) verbose = kwargs.pop('verbose',False) # Spectro display ax = kwargs.pop('ax',None) fig = kwargs.pop('fig',None) vmin = kwargs.pop('vmin',-120) vmax = kwargs.pop('vmax',-20) Sxx, tn, fn, ext = spectrogram(s, fs, window, nperseg, overlap, fcrop, tcrop, mode, verbose, display=False, **kwargs) # dB scale Sxx = 10*np.log10(Sxx) # format rois to bbox rois = format_features(rois, tn, fn) # if rois becomes empty after format_features if rois.empty : print ('rois is empty') features = rois return features else: print('number of rois : %d' % len(rois)) # Compute shape features and centroid features shape, params = shape_features(Sxx, resolution, rois) shape = format_features(shape, tn, fn) if verbose : print('Dataframe with shapes features') print(shape) centroid = centroid_features(Sxx, rois) centroid = format_features(centroid, tn, fn) if verbose : print('Dataframe with centroids features') print(centroid) if display : # view bbox ax, fig = overlay_rois(Sxx, rois, vmin=vmin, vmax=vmax, **kwargs) # view centroids overlay_centroid(Sxx, centroid, savefig=None, fig=fig, ax=ax, **kwargs) # plot shape plot_shape(shape.filter(regex=r'^shp_').mean(), params) # combine into a single df without columns duplication features = pd.concat([rois, shape, centroid], axis=1, sort=False) features = features.loc[:,~features.columns.duplicated()] if verbose : print('list of features') print(list(features)) return features
if __name__ == "__main__": import doctest doctest.testmod()