R/wavelet-T-detection.R

In Teresa00/hfmAnnotation: HFM Internal Package for Annotations

# Teresa updated April 19th 2018

#' @importFrom assertthat assert_that
#' @importFrom magrittr %>%
#'
library("magrittr")
library("purrr")
library("assertthat")


# some convenience constants for describing search directions...
BEFORE <-
  TRUE

AFTER <-
  FALSE

POSITIVE <-
  TRUE

NEGATIVE <-
  FALSE

extrema_around_index_Twave_ <-
  function(w_scale, pivot_index, search_window, search_direction, amplitude_threshold) {

    polarity <- ifelse(search_direction$is_positive, +1, -1)

    slice_start <- ifelse(search_direction$is_before, search_window$begin, pivot_index)

    slice_end <- ifelse(search_direction$is_before, pivot_index, search_window$end)

    wavelet_slice <-
      w_scale[(slice_start + 1):slice_end]

    extrema_set <-
      find_modulus_maxima_PT(wavelet_slice, 0, polarity)

    extrema <-
      if (length(extrema_set) != 0) {
        define_extrema_(w_scale, slice_start + extrema_set, search_direction$is_positive)
      } else {
        if(search_direction$is_before &&
           (pivot_index != search_window$begin) &&
           (sign(w_scale[search_window$begin]) == polarity)){

          list(location = search_window$begin,
               amplitude = w_scale[search_window$begin])

        }else if(! (search_direction$is_before) &&
                 (pivot_index != search_window$end) &&
                 (sign(w_scale[search_window$end]) == polarity)){

          list(location = search_window$end,
               amplitude = w_scale[search_window$end])

        }else{
          list(location = 0, amplitude = 0)
        }
      }

    # Amplitude Protection
    if(extrema$amplitude < amplitude_threshold){
      list(location = 0, amplitude = 0)
    }else{
      extrema
    }
    # ifelse(extrema$amplitude < amplitude_threshold,
    #        list(location = 0, amplitude = 0),
    #        extrema)
}


define_morphology_Twave_ <-
  function(w_scale, w_zerocross, search_window, t_slopes) {

    ts <-
      t_slopes

    PM <-
      ts$positive_maximum

    NM <-
      ts$negative_maximum

    is_minima <-
      ts$abs_max_amp > ts$abs_min_amp

    is_maxima <-
      !is_minima

    is_min_before <-
      PM < NM

    is_min_after <-
      PM >= NM

    is_max_before <-
      NM < PM

    is_max_after <-
      NM >= PM

    amplitude_threshold <-
      0.125 * ifelse(is_minima, ts$abs_max_amp, ts$abs_min_amp)

    find_min_around_pm_ <-
      function(is_before) {
        search_direction <-
          list(is_before = is_before, is_positive = FALSE)

        extrema_around_index_Twave_(
          w_scale,
          PM,
          search_window,
          search_direction,
          amplitude_threshold
        )
      }

    find_max_around_nm_ <-
      function(is_before) {
        search_direction <-
          list(is_before = is_before, is_positive = TRUE)

        extrema_around_index_Twave_(
          w_scale,
          NM,
          search_window,
          search_direction,
          amplitude_threshold
        )
      }

    min_before_pm <-
      if(is_minima){
        find_min_around_pm_(BEFORE)
      }else{
        list(location = 0,amplitude = 0)
      }

    min_after_pm <- if(is_minima){
      find_min_around_pm_(AFTER)
    }else{
      list(location = 0,amplitude = 0)
    }

    max_before_nm <- if(is_maxima){
      find_max_around_nm_(BEFORE)
    }else{
      list(location = 0,amplitude = 0)
    }

    max_after_nm <- if(is_maxima){
      find_max_around_nm_(AFTER)
    }else{
      list(location = 0,amplitude = 0)
    }

    if(is_minima){
      if(min_before_pm$amplitude == 0){
        if(min_after_pm$amplitude == 0){
          # Upwards T wave
          type_of_t <- 5
          if(min_before_pm$location == 0){
            t_wave <- peak_before(w_scale[search_window$begin:PM], PM)

            t_peak_onset <- t_wave
          }else{
            t_wave <- detect_zerocross_PT_(w_zerocross, min_before_pm$location, PM)

            t_peak_onset <- min_before_pm$location
          }

        }else{
          # Normal T wave
          type_of_t <- 1

          t_wave <- detect_zerocross_PT_(w_zerocross, PM, min_after_pm$location)

          t_peak_onset <- PM
        }
      }else{
        if (min_after_pm$amplitude == 0) {
          # Inverted T wave
          type_of_t <- 2

          t_wave <- detect_zerocross_PT_(w_zerocross, min_before_pm$location, PM)

          t_peak_onset <- min_before_pm$location
        } else {
          # Biphasic T wave (-/+)
          type_of_t <- 4

          t_wave <- detect_zerocross_PT_(w_zerocross, min_before_pm$location, PM)

          t_peak_onset <- min_before_pm$location
        }
      }
    }else{
      if (max_before_nm$amplitude == 0) {
        if (max_after_nm$amplitude == 0) {
          # Downwards T wave
          type_of_t <- 6

          if (max_before_nm$location== 0) {
            t_wave <- peak_before(w_scale[search_window$begin : NM], NM)

            t_peak_onset <- t_wave
          } else {
            t_wave <- detect_zerocross_PT_(w_zerocross, max_before_nm$location, NM)

            t_peak_onset <- max_before_nm$location
          }
        } else {
          # Inverted T wave
          type_of_t <- 2

          t_wave <- detect_zerocross_PT_(w_zerocross, NM, max_after_nm$location)

          t_peak_onset <- NM
        }
      } else {
        if (max_after_nm$amplitude == 0) {
          # Normal T wave
          type_of_t <- 1

          t_wave <- detect_zerocross_PT_(w_zerocross, max_before_nm$location, NM)

          t_peak_onset <- max_before_nm$location
        } else {
          # Biphasic T wave (+/-)
          type_of_t <- 3

          t_wave <- detect_zerocross_PT_(w_zerocross, max_before_nm$location, NM)

          t_peak_onset <- max_before_nm$location

        }
      }
    }

    list(
      t_wave = t_wave,
      t_peak_onset = t_peak_onset,
      type = type_of_t
    )
}


#' Search for T wave in ECG signals based on wavelet transform
#' Return a vector of index (based on processed ECG signal)
#' and corresponding type of morphology
#' 1 - normal; 2 - inverse; 3 - biphasic +/-; 4 - biphasic -/+;
#' 5 - upwards; 6 - downwards
#'
#' TODO: Add details, perhaps better description
#'
#' @param R_peaks An index vector of R peaks
#' @param wavelets The n * 5 wavelet transform matrix of ECG signal
#' @param alignments An index vector for aligning indexes
#' @param sampling_frequency signal alignmentsling rate, in alignmentsles per second
#'
#' @export
#'

wavelet_T_detection <-
  function(R_peaks, wavelets, alignments, sampling_frequency = 250) {
    scale4_wavelets <- wavelets[, 4]
    scale5_wavelets <- wavelets[, 5]

    R_peaks <-
      R_peaks - alignments[1] + 1

    T_wave <-
      c()

    T_onset <-
      c()

    type_of_T <-
      c()

    for (i in 1:length(R_peaks)) {
      search_window <-
        list(
          begin =
            R_peaks[i] + round(0.1 * sampling_frequency),

          end =
            R_peaks[i] + round(0.4 * sampling_frequency)
        )

      fourth_scale_slopes <-
        PT_slopes_searching_(scale4_wavelets, R_peaks, i, search_window, 4)

      t_slopes <-
        if (fourth_scale_slopes$does_wave_exist) {
        fourth_scale_slopes
        } else {
        PT_slopes_searching_(scale5_wavelets, R_peaks, i, search_window, 5)
        }

      if (t_slopes$does_wave_exist) {
        zerocross_scale <-
          t_slopes$search_scale - 1

        w_scale <-
          wavelets[, t_slopes$search_scale]

        w_zerocross <-
          wavelets[, zerocross_scale]

        ##################### Define T wave morphologies
        t_list <-
          define_morphology_Twave_(w_scale, w_zerocross, search_window, t_slopes)

        # search onset of T wave
        onset_search_window <-
          w_scale[max(search_window$begin, t_list$t_peak_onset - round(0.1 * sampling_frequency)):t_list$t_peak_onset]

        t_wave_onset <-
          wavelet_onset_search(t_list$t_peak_onset, onset_search_window, 0.25)




        T_wave[i] <- t_list$t_wave
        T_onset[i] <- t_wave_onset

        # c(T_wave[i], T_onset[i], type_of_T[i]) %<-%
        #   if (!is.nan(t_list$t_wave)) {
        #     list(T_detect, t_wave_onset, t_list$type)
        #   } else {
        #     define_morphology_Twave_(w_scale, w_zerocross, search_window, t_slopes)
        #   }

        type_of_T[i] <-
          if (is.nan(T_wave[i])) {
            NaN
          } else {
            # type_of_T[i]
            t_list$type
          }
      } else {
        T_wave[i] <-
          NaN

        T_onset[i] <-
          NaN

        type_of_T[i] <-
          NaN
      }
    }

    list(
      T_wave = T_wave + alignments[1] - 1,
      T_onset = T_onset + alignments[1] - 1,
      type_of_T = type_of_T
    )
  }