A Novel Technique for ECG Morphology Interpretation and Arrhythmia Detection Based on Time Series Signal Extracted from Scanned ECG Record

,


Introduction
Cardiovascular disease (CVD) is the one of the biggest health problem in Indian and around the world as well.Electrocardiogram is a traditional method used for the diagnosis of heart diseases for about a century.Maintaining and retrieving patient history during a course of treatment is a essential but a laborious process.More particularly, over a decade ago thermal ECG records were stored physically, off late, due to advancement in technology; it has been stored as scanned ECG images.Storing the scanned ECG trace images requires considerable storage space.This necessitated the development of an automated solution capable of storing the ECG digitally, retrieving it quickly and detecting cardiac arrhythmia automatically.Majority of the ECG's clinical information is said to be found in the intervals and amplitudes defined by its features (characteristic wave peaks and time durations).According to author's knowledge, very few researchers [Lawson et al., 1995, Silva et al., Wang et al., 2009, Chebil et.al., 2008, Kao et al.,2001] have approached the extraction of ECG digital time series signal from scanned ECG trace images.Lawson et al., chose a scanning resolution of 200 dpi and used global thresholds to separate the ECG trace from the background grid lines.The low resolution results in loss of data accuracy and global thresholds results in missing pixels which are replenished by linear interpolation.Fabio Badilini et al., 2005 developed an application for extraction of the ECG trace from the image.But the method requires the user to fix anchor points for missing peaks and thus the accuracy comes down.Shen et al., separated the ECG trace from the background grids using the histogram.The missing pixels are replenished by checking the value of the pixel in the original image.This is a tedious process.Kao et al., employed a color filter to remove the background gridlines in the color image.There was a problem of missing pixels in the process which was replenished by linear interpolation.Jalel Chebil et al., performed a comparative study of the extracted trace accuracy by scanning the image at various resolutions.Global thresholds and median filtering were employed to remove background grids.The threshold to separate the trace from the background should be selected based on the nature of the image to avoid any missing pixels.All conventional techniques use morphological operations such as erosion, dilation, thinning etc [Rafael C. Gonzalez and Richard E.Woods 2008] to extract the ECG trace from the background.However, all the above work addresses the issue of onedimensional time series signal alone.In this work, we propose an improved methodology to extract the digitized ECG time series dignal from scanned ECG records.As a novelty, we are using the Radon transform for deskewing the scanned images.Even though the conventional morphological techniques are adapted, they are applied in an iterative fashion on the binarized de-skewed images.This results in more accurate extraction of the time series trace.Further, a simple and useful way of axis identification is proposed.In addition to ECG digitization, we have extended this work to ECG morphological extraction and report generation.As a novelty, we have applied slope method for morphological extraction that eliminates the pre-processing of noise and baseline wandering technique.This method reduces the retrieval and computational time and improves the accuracy of ECG image.This chapter is divided into three subtopics: Converting thermal ECG trace to Digital ECG signal, Report generation from ECG morphology and automated arrhythmia detection.This chapter explains various techniques, adapted to extract the digital time series signal from scanned thermal ECG records.This process of digitally converting ECG trace reduces the storage space and retrieval time with increased viewing accuracy.The challenges here are as follows.Firstly, the analysis algorithms requires the ECG signal to be digital.Therefore, the conversion of scanned ECG records to digital time series is a pre requisite.The algorithms are also capable of handling data from the digital ECG device that provide a digital signal as an output.Secondly, the digital signal from scanned ECG requires standardization, i.e. based on the ECG records, the ECG digital signal must be re-sampled and voltage levels adjusted automatically.In addition to digital signal conversion, this chapter explains the technique used to interpret ECG morphology and generate reports based on the interpretation.A challenge in report generation is estimation of time and amplitude level from pixel information and ECG morphological interpretation techniques.After ECG morphology analysis, automatic cardiac arrhythmia classification is performed for diagnosis.In this automated cardiac arrhythmia detection we would discuss about various classification technique and there efficiency.

Methodology
The following sections explains in detail the various stages involved in capturing the ECG trace, its storage and retrieval, signal extraction and digital signal generation, report generation and finally abnormality classification.

Data acquisition and image processing
Figure 1 gives an overview of the image processing techniques involved in the signal extraction process.12 lead ECG signals were recorded at a paper speed of 25mm/sec and printed in thermal paper.These stored paper ECG trace is scanned at resolution of 600 dpi (dots per inch) black and white images and stored in jpeg format.Radon transform is applied on these images to detect and correct the skewness, which is incurred during the scanning process.The de-skewed image is adaptively binarized by choosing local thresholds.To limit the area to be binarized, the image is iteratively filtered by morphological filters.Each time the image retained is a cropped version of the original image.An envelope detection operation has been performed on the resultant binary ECG image to yield the upper and lower boundaries.The pixel values are then averaged to obtain the Digital ECG signal.The digital signal is windowed and re-sampled in accordance with the ECG record as shown in Figure .1.

Scanning and standardization
Table 1, lists out various scanning resolution and formats that that can be handled by the algorithm proposed.Based on the resolution of scanning and paper speed the pixel is defined in terms of time and amplitude units.For example, a resolution of 600 dpi implies 600 pixels in an inch (25.4 mm).Thus the number of pixels per mm can be calculated.Paper speed of 25 mm/sec i.e. 1 sec = 25 mm and calibration mark of 1mV amplitude = 10 mm is used to evaluate the value of each pixel in the time scale and amplitude scale.Pixel value in time scale is found to be 1.693 ms and amplitude scale to be 4.233 mV.These values are used during the digital time series signal generation.

Skew detection and correction using radon transform
Scanning process of paper ECG may results in skewness in scanned images either due to human error or faulty scanners.In order to extract faithfully the ECG signal from images, the skewness has to be eliminated.To remove the skew we have applied Radon transform [Lins R. D. andÁvila B. T., 2004, Prashanth et al., 2010] to find the angle of skewness.The skew angle has been selected, based on the maximum variance.

Adaptive binarization and iterative morphological operations
Our next objective is to binarize the image.Extracting the ECG signal from the image depends on the accuracy with which it is separated from the rest of the attributes present in the image like grid lines, textual characters etc. From elaborate experimentation, it is observed that, using various image processing filters and tweaking the thresholds, could not eliminate the noise completely from the ECG signal.
In this work Otsu's algorithm [Otsu N., 1979] has been performed for image adaptive binarization.Adaptive threshold technique for image binarization yields better results compared to global thresholds.This process of adaptive binarization ensures that the threshold is selected based on an active signal region using morphological operation and not on the entire image.Morphological operations include dilation and erosion as shown in the figure 2 (c).Erosion operation on the binary image results in the loss of ECG signal as shown in Figure 2c.However, during this process of erosion, we record the upper and lower limit in the Cartesian coordinate.The boundary limit values are assigned as threshold and the image clipping operation has been performed on the ECG trace.The clipped image is again fed back to the adaptive binarization algorithm and the whole process is repeated again.This methodology reduces the original image to the requisite binarized image containing the useful information.Further this has been achieved through reduced processing time.

Signal extraction and report generation
Once the ECG waveform in the image is separated from the gridlines, it must be converted to a digital format.Data logged as X-Y coordinates represent the signal.The binarized imaged is subjected to envelope detection to obtain a complete digital signal.

Envelope detection and axis identification
The result as shown in figure 2 (b), contains only the binary ECG trace whose thickness is more than a single pixel.An envelope detector is applied in order to obtain a time series.In an envelope detector, the image is scanned column wise, at each column the uppermost and lowermost non zero values are recorded.Plotting all the upper and lower bound values, we obtain upper and lower envelopes of the ECG signal respectively.Figure 3 (a) shows a original gray scale ECG trace.Figure 3 (b) shows its corresponding envelopes.
The mean of ECG signal is represented as Where, X is the mean ECG signal, Xub and Xlb is upper and lower envelope of ECG signal respectively  Axis identification plays a vital role in further diagnosis and automatic report generation of the ECG records.The test square pulse present in the starting of any ECG trace is used as reference for the axis.However, in most practical scanning procedures and data capture, square pulses are often absent due to the sheer length of the ECG paper.In order to overcome this, a novel and simple technique to identify the axis is proposed in this paper.
The obtained signal X can be represented as X = [x 1 , x 2 , x 3 ....x n ,], where the values of each element in the vector correspond to an ECG signal pixel location.By observation, it was found that the most significant and recurring pixels usually represent the axis of the signal along the horizontal.As the signal obtained can be treated as a vector, the axis is obtained by calculating the mode of the vector.Hence, we can represent it as: In most cases, the axis will be non-zero; therefore there is a need to offset the axis to the horizontal zero in order to standardize the signal.Equation 2 describes the offset process.

Normalising and R-peak detection
The R-peaks which have maximum amplitude in an ECG signal which is extracted by differentiating the ECG signal.Taking the first derivative of the ECG signal and discarding the negative values provides the location of the R-peaks.Subsequently, the number of peaks is used to calculate the heart rate.
Fig. 5. Normalized ECG signal with reference axis and location of the corresponding R-Peaks.Normalized signal with reference axis is used to extract the r-peaks.All peak locations are represented as 1's in the x-axis ranging between -to 3000.

ECG morphological feature extraction
A database of 25 patients paper ECG were recorded from 12 lead ECG machine is created and the digital signals are generated.The obtained ECG signals are processed to extract morphological features.The morphological feature extraction is carried out by two methods: time based method and slope based method.The morphological features extracted are P wave duration, QRS complex duration, T duration, PR interval, QT interval and ST intervals and P, R and T amplitudes.

Time based feature extraction
In this method, the digital ECG was filtered using a bandpass filter designed for a frequency range of 0.05 to 30 Hz. Obtained digital ECG signal is differentiated and the R peaks were identified.Heart Rate is calculated by mean of calculating the distance between two peaks.Heart rate is calculated using the formula 60 Heart Rate = RR Interval (4) In this method, morphological features are extracted by traversing the windowing function on either side of the R peak, based on the ascending and descending nature of the waveform the various peaks and onset and offset of each peak is identified.The window size for Q and S is 120 ms and for P and T is 150 ms.The Q and S peaks are found by traversing on the left and right side of the R peak within the specified window and locating the minimum or negative peak values.From the Q peak, by traversing on its left side, the maximum value is found to be the P peak.Similarly, by traversing to the right side of the S peak, the maximum value is found to be the T peak.The P_on and P_off points are identified by traversing the window function on either side of the P peak until they descend and reach the baseline.Similarly, T_on and T_off points are detected with respect to T peak.Using these data point various morphological features such as the duration of P, QRS and T waves, intervals such as PR, QT and ST and the amplitude of P, R and T waves were identified and marked on the digital signal plot as shown in Figure .6.The accuracy of 94.9% had been obtained for this method for a database of 25 patient's digial ECG records

Slope based feature extraction
In this method, slope of the ECG signal [Damodaran, et al., 2011]within a window size of 'n' number of samples is evaluated to extract the morphological details.The slope of the signal has both positive and negative values due to Increasing and decreasing peaks in an ECG waveform.Slope of the signal is calculated using Equation 5.
where i =1, 2...N-n, S(t) = Extracted ECG Signal with samples 1 to N and n=Window size S slope (t) = Slope signal The window size depends on the number of samples between the Q peak and R peak in the ECG signal.For finding the window size, the R peak is found by differentiating the ECG signal and the Q wave is detected as the negative peak immediately prior to the detected R peak.The window is placed at the 1 st sample and the slope between the 1 st and the (n+1) th sample is found and stored.The window is then placed on the 2 nd sample and the slope between the 2 nd and (n+2) th is found.The window is placed at all samples till the (N-n) th sample and the slope values found is stored as the S slope signal.
A standard range of values is defined for the inclination angle of the P wave, QRS complex and T wave for both normal and abnormal ECG.Thus from the defined range of slope values for the ECG waveform, the slope values between the minimum positive slope value and the maximum negative slope values are removed to eliminate any noise.
For finding the window size, the R peak is found by differentiating the ECG signal and the Q wave is detected as the negative peak immediately prior to the detected R peak.The slope of the signal within this window is found for the entire signal is shown in Figure .7. The first positive peak is the P_on, the first negative peak is P and the following zero crossing is P_off.Similar procedure is followed to identify the Q, R and S peaks and T wave.The features extracted using slope method is marked on the signal plot as shown in Figure .8. Accuracy of 97.09% had been obtained for this method for a database of 25 patient's digital ECG records As an extension, we have also performed a comparitive study between the slope method and the time base method as shown in table 2.

Automated arrhythmia detection
The morphological features extracted were used to detect the arrhythmias.In this study we have considered three different abnormalities, namely, Sinus Bradycardia, Sinus Tachycardia and PVC.Feature parameter chosen were P wave duration, QRS complex durations, T wave duration, PR intervals, QT intervals and ST intervals.Two classifiers are used to detect arrhythmias, namely Dynamic time warping (DTW) and Adaboost and their performance were compared.

Dynamic Time Warping (DTW)
The DTW classifier [Niranjan et.al., 2004, Venkatesh N andSrinivasan J.,2011] is based on the ranking of the prototypes by the distance to the query.Let, F = (f1…..fn) and G = (g1…..gm) be two time series of length n and m, respectively.To align the two sequences using DTW, we construct an n-by-m matrix whose (i,j)th element is the Euclidean distance d(i,j) between two points fi and gj.The (i,j)th matrix element corresponds to the alignment between the points fi and gj.A warping path, R is a contiguous sets of matrix elements that defines a mapping between F and G and is written as R={r1…..rS} where, max(m, n) < S < m + n -1.To limit the warping path, several constraints such as boundary conditions, continuity, monotonicity, and windowing [Bishop, 2006] are used.The DTW algorithm finds the point-to-point correspondence between the curves, which satisfies the above constraints and yields the minimum sum of the costs associated with the matching of the data points.There are exponentially many warping paths that satisfy the above conditions.The path that minimizes the warping cost is, The warping path can be found efficiently using dynamic programming to evaluate a recurrence relation, which defines the cumulative distance   , i j  up to the element (i, j) as the sum of d(i, j), the cost of dissimilarity between the ith and the jth points of the two sequences and the minimum of the cumulative distances up to the adjacent elements: The classification procedure based on DTW yielded the following results.
where, 1 indicates the sample has been correctly classified.In this experiment, stumps are used as a weak classifier.For reassigning the weights to the weak classifier 5000 iterations were performed and this was experimentally found to yield better results.Because it may have potential advantages such as higher classification performance, more rapid recognition process time and extension of recognition features, Adaboost was applied for the detection of cardiac arrhythmia.Each class of ECG type i.e. normal or arrhythmic, a label +1 or -1 is assigned to it.A large number of weak classifiers around 5000 are chosen.Decision stumps are chosen for classification.Decision stumps make prediction based on the value of just a single input feature.The input value if greater than the prediction value then the feature vector belongs to one class else it belongs to another class.Initially a set of training vectors are fed for classification.Labels are assigned for each input.A set of testing vectors are given as inputs for classification.Based on the labels assigned to each of the testing vector, the classification or misclassification is decided.The Adaboost classifier is implemented and the classification results are as shown in Table 4.The sensitivity of the classifier is evaluated and the average sensitivity is found to be 99%.

Report generation
Based on the various morphological features extracted using the proposed method and the arrhythmia detection using classifiers, a report is generated for each patient record.This ECG Report consists of Heart Rate, Morphological features duration and arrhythmias will be listed as a report.

Conclusion
The conversion of the scanned ECG record to a digital time series signal has been performed by an improved method of binarisation accurately.The digital time series data obtained is scaled in terms of amplitude and time.The digital signal is further processed for ECG morphological extraction procedure, by two methods namely, time based and slope based methodology.The accuracy of both the methods is evaluated by comparing the obtained results with manually read data from the paper record.Slope method is more accurate than other methods, in addition, this method eliminate the base line correction issues, noise removal issues with ECG signals.Further, this work has been extended to classification of ECG using DTW and Adaboost classifier for arrhythmia detection.The paper ECG converted will be provided as report with consists of Heart Rate, Morphological features duration and arrhythmias.This could be an aid tool to physician and electronic medical record maintains.This can function as a second option tool for initial screening producer for ECG.

Acknowledgment
The authors wish to thank Mr.Balamuralidhar.P, Head, Innovation Lab, TCS, for his continuous support.

Fig. 3 .
Fig. 3. (a) shows original paper ECG image of size 3000x250 pixels and (b) shows its corresponding digital ECG signal.
Fig. 4. A typical ECG extraction process is as shown in figure 4 (a -c).(d) represents the signal with the reference axis plotted as a dotted horizontal line

Fig. 7 .
Fig. 7. Slope plot for the digital ECG extracted form the paper ECG.Plotting of the slope values result in three peaks, each one for P, QRS and T waves respectively

Table 1 .
Comparison of the time based method and slope method for a single patient ECG patient record In this study, multiclass adaboost has been used to identifying the arrhythmias detection.Adaboost classifier increases the accuracy of weak classifier by reinforcing training on misclassified samples and assigns appropriate weights to each weak classifier.The final classification is given by

Table 5
FA), false acceptance rate (FAR) and false rejection rate (FRR) for different cases.In case-1, normal is made one class with PVC, Sinus Bradycardia and Sinus Tachycardia together is made into another class.Similarly in case-2, Sinus Bradycardia is one class, in case-3 Sinus Tachycardia is one class and in case-4 PVC is one class with the other three types together is the second class respectively.
presents the performance of the classification system for different arrhythmias.The performance of an arrhythmias detection is measured based on the confusion matrix with parameters false rejection (FR), false acceptance (

Table 4 .
Classification of ECG for different arrhythmias, Case -1 is normal, Case-2 is sinus Bradycardia, Case -3 sinus Tachycardia, and Case -4 is PVC.