Extraction of Peak Velocity Profiles from Doppler Echocardiography Using Image Processing

The objective of this study is to extract positive and negative peak velocity profiles from Doppler echocardiographic images. These profiles are currently estimated using tedious manual approaches. Profiles can be used to establish realistic boundary conditions for computational hemodynamic studies and to estimate cardiac time intervals, which are of clinical utility. In the current study, digital image processing algorithms that rely on intensity calculations and two different thresholding methods were proposed and tested. Image intensity histograms were used to guide threshold choices, which were selected such that the resulting velocity profiles appropriately represent Doppler shift envelopes. The resulting peak velocity profiles contained artifacts in the form of sudden velocity changes and possible outliers. To reduce these artifacts, image smoothing using a moving average process was then implemented. Bland–Altman analysis suggested good agreement between the two thresholding methods. Artifacts decreased when image smoothing was performed. Results also suggested that one thresholding method tended to provide the lower limit (i.e., underestimate) of velocities, while the second tended to provide the velocity upper limit (i.e., overestimate). Combining estimates from both methods appeared to provide a smoother peak velocity profile estimate. The proposed automated approach may be useful for objective estimation of peak velocity profiles, which may be helpful for computational hemodynamic studies and may increase the efficiency of current clinical diagnostic tools.


Introduction
Cardiovascular disease is one of the leading causes of death in the world [1]. Doppler echocardiography is a common method used to provide quantitative and qualitative diagnostic information about blood flow velocity and direction. The measured velocity profiles are used in clinical settings to determine cardiac time intervals, which would provide useful information for various cardiovascular conditions [2][3][4][5]. Cardiac time intervals can be used to predict major cardiovascular events (MACE) including cardiac ischemia, heart failure, or cardiac death. These velocity profiles also have utility in research studies of cardiac physiology and hemodynamics. For example, profiles are useful in elucidating the relationship between physiological processes and seismocardiographic signals [6][7][8]. Profiles can also be used to develop realistic boundary conditions in computational fluid dynamics simulations of the heart and great vessels and can potentially lead to more accurate predictions of important local and systemic hemodynamic changes [9,10]. Peak flow velocities in different cardiac phases, including systolic, early, and late diastolic can also be measured from these velocity profiles [11].
Image processing has been used to analyze echocardiographic images including Doppler echocardiography [12][13][14][15][16][17]. In this paper, a study for the automatic extraction of temporal distribution of peak velocity profile from Doppler echocardiography images is presented. The Doppler principle is extremely helpful in evaluating the physiology and pathophysiology of the organ under investigation. Ultrasound machines generally acquire Doppler images of moving targets. However, tedious manual measurements are required to measure peak velocities and various cardiac time intervals. An algorithm, such as the one presented here, can automate these measurements and thereby improve workflow significantly. For example, on average, a complete pediatric echo consists of 20 manual measurements of Doppler velocities and time intervals. It becomes clear how efficient the acquisition would be with an automated algorithm. The proposed algorithm can be especially useful in measurements of aortic ejection time and aortic velocity time integral. These measurements correlate well with cardiac output and are obtained in every diagnostic study. This paper is organized as follows: Section 2 provides image acquisition methods as well as the proposed image processing algorithm; results are reported in Section 3; and conclusions are discussed in Section 4.

Materials and Methods
In Doppler echocardiography, the speed and direction of blood flow within the heart are determined noninvasively using the Doppler effect (a distinct advantage to invasive cardiac catheterization, expensive cardiac MRI, and without the radiation or dye reaction risks of CT angiography). In Doppler cardiography, the ultrasound beam should be parallel to the blood flow, which is a potential limitation of the method. Ultrasound machines make use of the Doppler shift principle to determine the direction and velocity of a structure (tissue or blood) based on the formula, ∆f = f t − f r , where f t and f r are the transmitted and reflected frequency, respectively. The motion of the structure toward the ultrasound transducer results in compression of the reflected waves, reducing its period (i.e., negative ∆T, where T is the period) and thus increasing their frequency (i.e., negative ∆f ). In the same manner, a positive ∆T would indicate an object moving away from the transducer. By convention, negative ∆T is displayed above the baseline and positive ∆T below the baseline. These principles are used routinely by ultrasound systems.
In this study, an image processing algorithm was developed to extract peak velocity profiles from the echocardiographic images. The algorithm was based on identification and thresholding of gray pixels in the echocardiography images ( Figure 1). Pulsed Doppler echocardiography videos were recorded from a cardiovascular ultrasound machine (Model: EPIQ 5, Philips, Netherlands) without moving or tilting the probe at a frame rate of 30 frames per second. The pulsed Doppler was acquired at sweep speed of 75 mm/s. The echocardiography video frames were then read as eight-bit unsigned integer images using MATLAB (R2015b, The MathWorks, Inc., Natick, MA, USA). Therefore, the color components of these images were integers in the range [0, 255]. Eight-bit images were used since they would provide enough information (e.g., tonal range and resolution) for estimating the peak velocity profiles in this study.
Bioengineering 2019, 6, x FOR PEER REVIEW 2 of 11 velocities in different cardiac phases, including systolic, early, and late diastolic can also be measured from these velocity profiles [11]. Image processing has been used to analyze echocardiographic images including Doppler echocardiography [12][13][14][15][16][17]. In this paper, a study for the automatic extraction of temporal distribution of peak velocity profile from Doppler echocardiography images is presented. The Doppler principle is extremely helpful in evaluating the physiology and pathophysiology of the organ under investigation. Ultrasound machines generally acquire Doppler images of moving targets. However, tedious manual measurements are required to measure peak velocities and various cardiac time intervals. An algorithm, such as the one presented here, can automate these measurements and thereby improve workflow significantly. For example, on average, a complete pediatric echo consists of 20 manual measurements of Doppler velocities and time intervals. It becomes clear how efficient the acquisition would be with an automated algorithm. The proposed algorithm can be especially useful in measurements of aortic ejection time and aortic velocity time integral. These measurements correlate well with cardiac output and are obtained in every diagnostic study. This paper is organized as follows: Section 2 provides image acquisition methods as well as the proposed image processing algorithm; results are reported in Section 3; and conclusions are discussed in Section 4.

Materials and Methods
In Doppler echocardiography, the speed and direction of blood flow within the heart are determined noninvasively using the Doppler effect (a distinct advantage to invasive cardiac catheterization, expensive cardiac MRI, and without the radiation or dye reaction risks of CT angiography). In Doppler cardiography, the ultrasound beam should be parallel to the blood flow, which is a potential limitation of the method. Ultrasound machines make use of the Doppler shift principle to determine the direction and velocity of a structure (tissue or blood) based on the formula, ∆f = ft − fr, where ft and fr are the transmitted and reflected frequency, respectively. The motion of the structure toward the ultrasound transducer results in compression of the reflected waves, reducing its period (i.e., negative ∆T, where T is the period) and thus increasing their frequency (i.e., negative ∆f). In the same manner, a positive ∆T would indicate an object moving away from the transducer. By convention, negative ∆T is displayed above the baseline and positive ∆T below the baseline. These principles are used routinely by ultrasound systems.
In this study, an image processing algorithm was developed to extract peak velocity profiles from the echocardiographic images. The algorithm was based on identification and thresholding of gray pixels in the echocardiography images ( Figure 1). Pulsed Doppler echocardiography videos were recorded from a cardiovascular ultrasound machine (Model: EPIQ 5, Philips, Netherlands) without moving or tilting the probe at a frame rate of 30 frames per second. The pulsed Doppler was acquired at sweep speed of 75 mm/s. The echocardiography video frames were then read as eight-bit unsigned integer images using MATLAB (R2015b, The MathWorks, Inc., Natick, MA, USA). Therefore, the color components of these images were integers in the range [0, 255]. Eight-bit images were used since they would provide enough information (e.g., tonal range and resolution) for estimating the peak velocity profiles in this study.  The images were considered as a scalar function f (X, Y), where X ∈ x = [1, . . . , m] and Y ∈ y = [1, . . . , n], and m and n are the number of pixels in the horizontal and vertical directions, respectively. Figure 2a shows a sample Doppler echocardiography image taken from the center of the mitral valve of a 55-year-old male with no history of heart disease. This image shows the four-chamber view of the heart, where the Doppler sample volume is placed in mitral valve inflow. The Doppler flow above the baseline depicts the inflow (from the left atrium to the left ventricle), and the Doppler flow below the baseline is the outflow (from the left ventricle toward the aorta). In Figure 2a, m and n are 1920 and 1080, respectively. The Doppler baseline (i.e., zero velocity) is shown as a yellow solid line at Y base = 790. Positive and negative shifts are displayed above and below this line and are indicative of blood flow toward and away from the transducer, respectively.
The proposed image processing algorithm for peak velocity profiles detection can be summarized in the following steps: • Step 1: Load the raw echocardiographic image (such as the one shown in Figure 2a).

•
Step 2: Define a region of interest (ROI) that encompasses the gray pixels containing the blood velocity information (1 < x < m, n up < y < n lo ), where n up and n lo are the upper and lower borders of the ROI, respectively ( Figure 2a) then remove pixels outside the ROI. For example, all pixels above n up = 575 and below n lo = 1030 in Figure 2a are removed. Figure 2b shows the resulting image. The histogram of the ROI is shown in Figure 2c. Once the ROI is defined, calculate the gray-scale intensity of pixels, µ, at pixel (X ∈ x, Y ∈ y). • Step 3: Find the maximum intensity, µ, at each vertical line (i.e., at each X ∈ x) for n up < y < n lo , which was called µ max,X .

•
Step 4: Detect gray pixel edges. Here, two different thresholding methods are attempted to detect the edge of the gray pixels in the ROI. In the first method (Figure 3a), edge search starts at the upper and lower borders of the ROI and moves toward the baseline. The upper and lower edges of the Doppler shift at each X are chosen as the smallest and largest Y ∈ (n up < y < n lo ) such that µ(X, Y) > θ 1 µ max,X , where θ 1 is the threshold value, µ(X, Y) is the intensity of the pixel located at (X, Y), and µ max,X is the maximum intensity at time X. This method is designed to provide a lower limit of the velocity magnitude.
Bioengineering 2019, 6, x FOR PEER REVIEW 3 of 11 The images were considered as a scalar function f (X, Y), where X ∈ x = [1, …, m] and Y ∈ y = [1, …, n], and m and n are the number of pixels in the horizontal and vertical directions, respectively. Figure 2a shows a sample Doppler echocardiography image taken from the center of the mitral valve of a 55-year-old male with no history of heart disease. This image shows the four-chamber view of the heart, where the Doppler sample volume is placed in mitral valve inflow. The Doppler flow above the baseline depicts the inflow (from the left atrium to the left ventricle), and the Doppler flow below the baseline is the outflow (from the left ventricle toward the aorta). In Figure 2a, m and n are 1920 and 1080, respectively. The Doppler baseline (i.e., zero velocity) is shown as a yellow solid line at Ybase = 790. Positive and negative shifts are displayed above and below this line and are indicative of blood flow toward and away from the transducer, respectively.
The proposed image processing algorithm for peak velocity profiles detection can be summarized in the following steps: • Step 1: Load the raw echocardiographic image (such as the one shown in Figure 2a).

•
Step 2: Define a region of interest (ROI) that encompasses the gray pixels containing the blood velocity information (1 < x < m, nup < y < nlo), where nup and nlo are the upper and lower borders of the ROI, respectively ( Figure 2a) then remove pixels outside the ROI. For example, all pixels above nup = 575 and below nlo = 1030 in Figure 2a are removed. Figure 2b shows the resulting image. The histogram of the ROI is shown in Figure 2c. Once the ROI is defined, calculate the gray-scale intensity of pixels, µ, at pixel (X ∈ x, Y ∈ y). • Step 3: Find the maximum intensity, µ, at each vertical line (i.e., at each X ∈ x) for nup < y < nlo, which was called µmax,X. • Step 4: Detect gray pixel edges. Here, two different thresholding methods are attempted to detect the edge of the gray pixels in the ROI. In the first method (Figure 3a), edge search starts at the upper and lower borders of the ROI and moves toward the baseline. The upper and lower edges of the Doppler shift at each X are chosen as the smallest and largest Y ∈ (nup < y < nlo) such that µ(X, Y) > θ1 µmax,X, where θ1 is the threshold value, µ(X, Y) is the intensity of the pixel located at (X, Y), and µmax,X is the maximum intensity at time X. This method is designed to provide a lower limit of the velocity magnitude.   In the second method (Figure 3b), the edge search starts at the baseline and moves toward the ROI borders. The upper edge at each X is found as the smallest Y ∈ (nup < y < Ybase) such that µ(X, Y) < θ2 µmax,X. The lower edge is calculated as the largest Y ∈ (Ybase < y < nlo) such that µ(X, Y) < θ2 µmax,X. This method is designed to provide an upper limit of the velocity magnitude.
The intensity histograms of the ROI (Figure 2c) consisted of two main peaks around gray intensities of 0 and 175. These peaks are marked with red arrows in Figure 2c. The first peak corresponds to the black background. The threshold values for both methods, θ1 and θ2, were chosen to be between the two peaks such that they result in velocity profiles that agree with the expert opinion of two study specialists with more than 30 years of experience.

•
Step 5: Construct the positive peak velocity profile by connecting the upper edges detected in Step 4 above. Then construct the negative peak velocity profile by connecting the lower edges. The pseudo-code for the proposed algorithm is provided in Algorithm 1.  In the second method (Figure 3b), the edge search starts at the baseline and moves toward the ROI borders. The upper edge at each X is found as the smallest Y ∈ (n up < y < Y base ) such that µ(X, Y) < θ 2 µ max,X . The lower edge is calculated as the largest Y ∈ (Y base < y < n lo ) such that µ(X, Y) < θ 2 µ max,X . This method is designed to provide an upper limit of the velocity magnitude.
The intensity histograms of the ROI (Figure 2c) consisted of two main peaks around gray intensities of 0 and 175. These peaks are marked with red arrows in Figure 2c. The first peak corresponds to the black background. The threshold values for both methods, θ 1 and θ 2 , were chosen to be between the two peaks such that they result in velocity profiles that agree with the expert opinion of two study specialists with more than 30 years of experience.

•
Step 5: Construct the positive peak velocity profile by connecting the upper edges detected in Step 4 above. Then construct the negative peak velocity profile by connecting the lower edges. The pseudo-code for the proposed algorithm is provided in Algorithm 1. The code for the proposed algorithm in this paper is available online at: https://github.com/mirtatae/ DopplerVelocityProfileExtraction.
The images were also smoothed by convolving with a square step kernel and the peak velocity profiles were estimated using the proposed algorithm. The pixels' intensity was smoothed by a moving average process. Here, the smoothed image intensity, µ smooth , at (X, Y) was calculated as the mean intensity of a (2p + 1) × (2q + 1) rectangle centered at (X, Y): where µ was the image intensity before smoothing. The peak velocity profiles detected from the original and smoothed images were compared in order to investigate the effect of image smoothing. The E waves and A waves were then measured from the extracted profiles using the proposed algorithm and compared with the values manually calculated by the study specialists. Edges of the Doppler shift (i.e., peak velocity profiles) were also detected using a standard edge detection method, Canny approximation. In this method, the gradient of the image intensity function is calculated using the derivative of a Gaussian filter [18]. The edges of the Doppler shift were then found by looking for local maxima of the gradient in the vertical direction. Canny approximation was selected since it is one of the edge detection methods that is less sensitive to noise [18]. The results were qualitatively compared with the results of our proposed algorithm.  = [1, . . . , m] and Y ∈ y = [1, . . . , n] 6: 7: # first criteria 8: for ∀ (X,Y) | X ∈ x and n up < Y < n lo 9: µ max,X ← max( µ(X,:) ) 10: if µ(X,Y) < θ 1 × µ max,X 11: µ(X,Y) ← 0 12: end if 13: k ← find all Y that µ(X,Y) > 0 14: µ max,X,up ← max( µ(X,Y base − a:Y base ) ) 22: i up ← Y | µ(X,Y) = µ max,X,up 23: k up ← find all Y that µ(X, n up : i up + Y base − a) < θ 2 × µ max,X,up 24: P 2,up ← max(k up ) 25: 26: µ max,X,lo ← max( µ(X,Y base :Y base + a) ) 27: i lo ← Y | µ(X,Y) = µ max,X,lo 28: k lo ← find all Y that µ(X, i lo + Y base : n lo ) < θ 2 × µ max,X,lo 29: P 2,lo ← min(k lo ) 30: end for

Results and Discussion
Two different thresholding methods were attempted to estimate the lower and upper limit of the peak velocity profile. Figure 4a shows the peak velocity profiles extracted from the Doppler image using both proposed thresholding methods. It can be seen that the lower limit of the peak velocity profile (red line in Figure 4a) was incorrectly estimated at the upper left edge of the ROI. In this portion of the ROI (1 < x < 200, 500 < y < 700 of the original image), there was some white text. This error resulted because the lower-limit thresholding method scanned the pixels at each vertical line starting from the borders of the ROI. Outside this ROI portion, the peak velocity profiles obtained from the two thresholding criteria showed general agreement, with the upper-limit velocity magnitude estimate being higher most of the time.
Although there was general agreement between the two estimates, both peak velocity profiles appeared to contain potential outliers and rapid velocity changes (as seen in Figure 4a). These results suggested adding a smoothing step to the algorithm as described in Section 2. Figure 4b shows the velocity profiles resulting from the smoothed Doppler image for both proposed thresholding methods. Here, the profiles appear to have less sudden changes compared to those of Figure 4a.
To quantify agreement between the two thresholding methods (i.e., lower and upper velocity limits, respectively), Bland-Altman analysis was carried out. Results are shown in Figure 5a,b for the negative and positive peak velocity profiles extracted from the original image, respectively. Information is also shown for the profiles extracted from smoothed image in Figure 5c,d, respectively. Table 1 lists the bias and standard deviation of the difference between the two estimates. It can be seen that the standard deviation decreased with smoothing possibly due to random noise reduction. The bias increasing with smoothing is an indication of more separation between the upper and lower limits, which would be consistent with reduced noise in the estimates.
The agreement was higher for the negative velocities, which may be due to the reduced outliers that are represented by the points outside the confidence interval in Figure 5. These outliers are more prominent for the positive velocities. Outliers may be found from these plots and eliminated from the data to reduce errors in future studies.
Bioengineering 2019, 6, x FOR PEER REVIEW 6 of 11 The agreement was higher for the negative velocities, which may be due to the reduced outliers that are represented by the points outside the confidence interval in Figure 5. These outliers are more prominent for the positive velocities. Outliers may be found from these plots and eliminated from the data to reduce errors in future studies.      . Fast Fourier transform of (a), (b) negative and positive peak velocity profiles using thresholding method 1, respectively; and (c), (d) negative and positive peak velocity profiles using thresholding method 2, respectively. Black and red lines show the spectra before and after smoothing, respectively. Each pixel in the Doppler image corresponded to 1.9 ms, which resulted in a sampling frequency of 532 Hz.
The spectra of the velocity profile of the non-smoothed and smoothed profiles are shown in Figure 6 and show significant reduction in high frequencies with smoothing. This further demonstrates the role of smoothing in reducing sudden velocity changes as expected.
It is worth mentioning that the two proposed thresholding methods estimate the peak velocity differently and necessitate different threshold values. The first method starts searching for the peak velocity at the ROI borders and uses a relatively larger threshold to avoid artifacts at the borders. This method estimates a lower limit (underestimate) for the peak velocity magnitude. In Figure 5. Bland-Altman analysis to assess the agreement between the two thresholding methods for (a), (b) negative and positive peak velocity profiles before smoothing, respectively; and (c), (d) negative and positive peak velocity profiles after smoothing, respectively. The solid line represents the mean value (i.e., bias) of the differences between two estimated profiles. The dashed lines show the 95% confidence interval (i.e., bias ± 1.96 SD). The unit of the horizontal and vertical axes is pixel, where each pixel corresponds to 0.34 cm/s. The spectra of the velocity profile of the non-smoothed and smoothed profiles are shown in Figure 6 and show significant reduction in high frequencies with smoothing. This further demonstrates the role of smoothing in reducing sudden velocity changes as expected.
Bioengineering 2019, 6, x FOR PEER REVIEW 7 of 11 Figure 5. Bland-Altman analysis to assess the agreement between the two thresholding methods for (a), (b) negative and positive peak velocity profiles before smoothing, respectively; and (c), (d) negative and positive peak velocity profiles after smoothing, respectively. The solid line represents the mean value (i.e., bias) of the differences between two estimated profiles. The dashed lines show the 95% confidence interval (i.e., bias ± 1.96 SD). The unit of the horizontal and vertical axes is pixel, where each pixel corresponds to 0.34 cm/s. Figure 6. Fast Fourier transform of (a), (b) negative and positive peak velocity profiles using thresholding method 1, respectively; and (c), (d) negative and positive peak velocity profiles using thresholding method 2, respectively. Black and red lines show the spectra before and after smoothing, respectively. Each pixel in the Doppler image corresponded to 1.9 ms, which resulted in a sampling frequency of 532 Hz.
The spectra of the velocity profile of the non-smoothed and smoothed profiles are shown in Figure 6 and show significant reduction in high frequencies with smoothing. This further demonstrates the role of smoothing in reducing sudden velocity changes as expected.
It is worth mentioning that the two proposed thresholding methods estimate the peak velocity differently and necessitate different threshold values. The first method starts searching for the peak velocity at the ROI borders and uses a relatively larger threshold to avoid artifacts at the borders. This method estimates a lower limit (underestimate) for the peak velocity magnitude. In Figure 6. Fast Fourier transform of (a), (b) negative and positive peak velocity profiles using thresholding method 1, respectively; and (c), (d) negative and positive peak velocity profiles using thresholding method 2, respectively. Black and red lines show the spectra before and after smoothing, respectively. Each pixel in the Doppler image corresponded to 1.9 ms, which resulted in a sampling frequency of 532 Hz.
It is worth mentioning that the two proposed thresholding methods estimate the peak velocity differently and necessitate different threshold values. The first method starts searching for the peak velocity at the ROI borders and uses a relatively larger threshold to avoid artifacts at the borders. This method estimates a lower limit (underestimate) for the peak velocity magnitude. In contradistinction, the second method starts searching for the peak velocity from the baseline (i.e., zero velocity) and uses a smaller threshold, which would estimate an upper limit (overestimate) of the peak velocity profile. A combined method was also tested in the current study where the estimated velocity profiles from both methods are averaged (Figure 7). The averaged profiles, P avg , were calculated as: P avg,up = P 1,up + P 2,up /2, P avg,lo = P 1,lo + P 2,lo /2, where P 1,up and P 1,lo are the upper and lower limits of the peak velocity profile from thresholding method 1, respectively. P 2,up and P 2,lo are the upper and lower limits of the peak velocity profile from thresholding method 2, respectively. The values of A and E waves are listed in Table 2. The last two columns show the R 2 and Bland-Altman bias ± 1.96 SD between the study specialists and other methods. Results suggested that the bias between the specialist velocity measurements and the proposed method was 1.98 ± 11.20 cm/s for velocities in the range of about 65-75 cm/s.   Figure 8b shows the detected edges using a higher threshold value of 0.20. Although the number of artifacts decreased, the algorithm still resulted in significant number of artifacts especially in the positive Doppler shift. An increased threshold (i.e., 0.30) was then used to reduce these extra unwanted edges (Figure 8c). As a result, the number of artifacts drastically decreased. However, there were still some artifacts around the baseline and Doppler shift borders. Figure 8d shows the results when a threshold of 0.40 was used. In this case, most artifacts were removed. However, due to a high threshold value, the algorithm missed many parts of the positive and negative peak velocity profiles resulting in disconnected profiles. Therefore, in this method, there is a compromise between avoiding artifacts and losing the peak velocity profile information. In contrast to the proposed method in this paper, the Canny method may find less or more than two edge points at each time (i.e., each vertical line). Finding more than two edge points at the same time is indicative of artifacts (such as artifacts seen in Figure 8a-c). Finding less than two points, indicates that the method did not properly estimate either one or both positive and negative Doppler shift envelopes. The proposed algorithm in this study, on the other hand, always estimates a lower and an upper edge in each vertical line. The proposed algorithm also resulted in continuous velocity profiles and less artifacts compared to the Canny method. In addition, it was simpler and less computationally expensive since it was only based on two simple thresholding operations. The proposed algorithm was also tested on 13 more images from the same subject. The results were consistent with the results shown in Figures 4 and 7.

Limitations and Future Work
The proposed algorithm in this study was employed to detect peak velocity profiles from the images recorded by only one echo machine. In future studies, this algorithm should be tested on more subjects and different echocardiography machines. In addition, the performance of this algorithm was compared with the manual calculations of the study clinicians. The results were also qualitatively compared with a standard edge detection method. However, there are other algorithms for Doppler trace detection in literature [19][20][21] that outperform the Canny edge detector. The proposed algorithm should be compared against these algorithms.

Limitations and Future Work
The proposed algorithm in this study was employed to detect peak velocity profiles from the images recorded by only one echo machine. In future studies, this algorithm should be tested on more subjects and different echocardiography machines. In addition, the performance of this algorithm was compared with the manual calculations of the study clinicians. The results were also qualitatively compared with a standard edge detection method. However, there are other algorithms for Doppler trace detection in literature [19][20][21] that outperform the Canny edge detector. The proposed algorithm should be compared against these algorithms.

Limitations and Future Work
The proposed algorithm in this study was employed to detect peak velocity profiles from the images recorded by only one echo machine. In future studies, this algorithm should be tested on more subjects and different echocardiography machines. In addition, the performance of this algorithm was compared with the manual calculations of the study clinicians. The results were also qualitatively compared with a standard edge detection method. However, there are other algorithms for Doppler trace detection in literature [19][20][21] that outperform the Canny edge detector. The proposed algorithm should be compared against these algorithms.

Conclusions
An image processing algorithm based on thresholding was proposed to extract peak velocity profiles from Doppler echocardiography images. Two different thresholding methods were used to estimate the upper and lower limits of the peak velocity profiles. Doppler images were smoothed using the moving average process to improve the quality of the extracted velocity profiles. The upper and lower limits of the velocity profiles were also averaged, and a combining method was proposed. Bland−Altman analysis suggested that the bias between the specialist velocity measurements and the proposed method was 1.98 ± 11.20 cm/s. The velocity profiles extracted by the proposed method can be used in clinical and research studies including cardiac time interval estimations. It can also be useful in providing boundary conditions for computational cardiovascular hemodynamics studies. Future work would compare the estimates obtained in the current study with additional standard edge detection methods along with validating these profiles using more accurate velocity profiles from invasive cardiac catheterization.