Echocardiography segmentation by Fractional Differential and Improved Canny, analysis by Fourier Descriptor

. The omnidirectional M-mode echocardiogram provides a new method for human heart functional analyses. In this article, to sharpen object edges, we designed image processing kernel based on Fractional differential for image enhancement. After that, the contour of the left ventricle in a short axis is first extracted using both an improved Canny edge detection algorithm and the gray level searching algorithm in the radial direction as auxiliary. The modified Canny edge detection algorithm with the matching method between adjacent frames then is adopted for the subsequent frames to extract the left ventricular contours. The non-functional movements in the B-ultrasonic plane are determined by using the movement extracting method based on Fourier descriptors and the mass center with the inertia axis method, and the movements are removed from a compound motion. The Fourier descriptors are applied to get a series of image contour curves with the principal translation and rotation. Hence the curve of the cardiac motion can accurately show functional movements in any location of the heart. Using our technique, we can reduce multi-lines and excursion, as well as correct the omnidirectional M-mode echocardiography.


Introduction
The heart is a vital organ in a human body and its pumping motion closely relates to its function. It is necessary to analyze the heart motion in the diagnosis of the heart disease [1][2] . The traditional M-mode echocardiography is used to analyze the myocardium function by displaying the movement-curve of the heart structure on the sampling line and by measuring various parameters of the wall motion (such as wall thickness, ventricular diameter, wall short-axis shortening fraction, ejection fraction, and cardiovascular exercise time, etc.) [3][4] . The main weakness of the detection is that the incidence angle of the ultrasound is not consistent with the direction of the movement, so obtained velocity error is great. The Anatomic M-mode (AMM) echocardiography is for post-processing an original image to obtain a motion curve of the interested local myocardium in wall or with the slice captured through the location of a sample line in the image [5][6] . The AMM echocardiography can only be set with three sampling lines in an image, but it cannot get the curves of the myocardial motion of the walls during the same cardiac cycle. Our omnidirectional M-mode echocardiogram system improved the AMM echocardiography system and imaged the myocardial motion of different parts with a number of motion curves synchronously during the same cardiac cycle, and the main work is carried out based on the improved Canny edge detector and the Fourier descriptors. However, with the effect of the non-functional movements in the heart caused by breathing, blood surging and mutual effects among organs, a heart wall does not always move along a sample line. Obviously, the obtained gray level versus time waveform cannot reflect the functional movements of the some parts of the heart accurately [7] .
In this paper, to make keep edges while smooth away noises, we proposed a new Fractional differential template with a circle data structure, the coefficients in a 5x5 template are different in different places. The studied object extraction method, based on the improved Canny edge detection algorithm and the gray scale search algorithm, extracts the left ventricular contours of the short axis, and corrects the non-functional movements of the heart by obtaining the cardiac translation and ventricular deflection in the B-plane by using the movement extracting algorithm based on the Fourier descriptors. The aim is to improve the accuracy of the acquired cardiac waveform graph and then to increase the precision of the obtaining cardiac motion parameters.
While for the B-mode ultrasound images of the cardiac short-axis, the motion tracking based on the ventricular contour is used to correct the non-functional movements of the heart in this study.
There has been extensive literature in the field of this research topic. Mihaela constructed a 3D MR image-based computer model of pathologic hearts, augmented with histology and optical fluorescence imaging to characterize the action potential propagation [8] ; and as the main information sources of magnetic resonance images, the detailed information can be discovered in the recent two review articles [9][10] .
Recently, some researchers were interested in using CT technique to make 3D information analysis [11][12][13][14] . In addition to the magnetic resonance and CT images, Ultrasound images can also be used to do diagnosis in real time, and they are intuitive, no damage, timely accessed and with other characteristics. The ultrasound images can provide the comprehensive information for physicians to understand the anatomical structure and dynamic process of the heart. Besides, the method based on ultrasound images is easy and flexible to operate, and with the low cost in the diagnosis of the heart diseases. So it is still one of the most important ways to check the heart disease [15] . The different model based algorithms or methods have been developed for the heart movement tracing or mapping in real time [16][17] . In addition, it can also record the M-mode ultrasound images of the multiple cardiac cycles in any structure in the heart [18] .
In the B-plane, the non-functional movements include the translation of the cardiac structure and deflection in the direction of the ventricular long axis. The motion tracking, image registration and other functions were commonly used to correct the non-functional movements of the heart [16][17][18] . The image registration was usually used for the medical images from different patients or from the same patient but in different times, which could detect the similarity of images but the result was not satisfactory for the cardiac image sequences with the ventricular systolic and diastolic movements [19] . The image tracking based on features was to track the whole or part movements in images to achieve the registration results. Although a genetic algorithm was adopted to track feature points for correcting the non-functional movements of the heart in the B-mode ultrasound images, it was difficult to find out all the points [20] . Instead, using other image region tracing and detection algorithms could also be valuable, such as the Wavelet based algorithm [21] and the recursive thresholding algorithm [22] .

Image enhancement by a new algorithm on Fractional differential
As usual, before image segmentation or object extraction, the image preprocessing is needed [23] .
To increase the detailed discontinuity information in our images, we made a new algorithm based on Fractional differential. Normally, the medical images are enhanced by using 1st or 2nd order differential algorithms such as Sobel, Robert, Laplacian and Canny etc. Although they can sharpen the edges and textures in images more or less, but meanwhile they will also increase noise very much, which make further image segmentation hard. As instance, the high pass operators often fail to enhance the vague or weak edges in images [24][25] . The fractional difference (order 0.3-0.6) can detect the vague edges and textures, and keep noise lower; hence, it is suitable for the complicated medical images. The key task is how to design a simple and applicable fractional differential kernel.
In a digital image, there is a large relevance among the pixel grey scales in the neighborhood of the detecting pixel. The pixels are auto-correlated, their fractal geometric information show up the properties of complex textural features, and the fractional differential is one of Fractal Geometry Mathematical Foundations. It is natural for one to think about what effect is gained when applying the fractional differential to detect textural features in an image. In this paper, a new fractional differential algorithm, which can improve detailed texture information remarkably, is studied.
In this way, the v order fractional order of the differential expression in one dimensional signal ( ) s t can be deducted as: These n +1 non-zero coefficient values are in order as: In a traditional Tians 5x5 template of Fractional calculus, as Fig.1 shown, except for the central point, the other points have only two different coefficients " 1 1 a → " and " The coefficient sum in Fig.1 If we set: 1 v = , then, ( ) If we set: It is very different to the Tians template, and the new algorithm has the 5x5 templates as shown in Fig.2. In the new template, with the square data structure, there are 5 different coefficients, and with the circle data structure, there 4 different coefficients. The different coefficient values are expressed in Eq. (7), and the coefficient sums in the square data structure and the circle data structure are 3 S (Eq.(8) and 4 S (Eq.(9)) respectively. Fig.5. New Fractional differential 5x5 templates with square and circle date structures.
The coefficient sum in Fig.2

Cardiac ventricular contour extraction and description
For the object tracking, there is a huge number of algorithms [26][27] , generally, they have been studied completely or partly either based on similarity or based on discontinuity, which is depending on image resolution [28] and the applications [29] . In our application, the algorithm is proposed based on discontinuity, i.e, the improved Canny edge detector + Fractional differential image enhancement template. The details for the improved Canny edge detector is described as the follows.

Improved Canny edge detection algorithm
In the B-mode ultrasound image in the cardiac short-axis, the radial variation of the heart chamber's contour reflects the function properties of the systolic and diastolic of the heart, and the overall translation and rotation in the B-plane of the contour reflect the non-functional movements of the heart. At present it is possible to get the cardiac ventricular contours with an improved boundary extraction algorithm. The time interval between two standard DICOM images with 75 frames per second is approximately 13ms, which is much shorter than the cardiac cycle 0.75s~1s of a normal person. So there is a strong spatial correlation between frames. For this reason, one can divide the extraction procedure of the ventricular contours into two stages: (1) the extraction of the ventricular contour in the first frame; and (2) the extraction of the ventricular contours in the subsequent frames. After the comparison among experiments with a variety of boundary extraction algorithms, it proves that a special edge detection based object tracing algorithm [20] with a number of post processing functions are definitely needed [24] . In this paper, an improved Canny edge detection algorithm [30][31] is studied by a modified gray scale search algorithm [30] based on radial, and it is adopted for the first frame to extract the ventricular contour. An improved Canny edge detection algorithm with the matching method between adjacent frames is adopted for the subsequent frames to get the left ventricular contours, which is described as follows.
In the Canny algorithm, the one-dimensinonal Gaussian filter is firstly discussed because the two-dimensional Gaussian filter can be superimposed by two one-dimensional Gaussian filters.
And the mathematical expression is: Where, 2 σ denotes the extent of the smoothed image. If it is too small, the noise cannot be eliminated completely. On the contrary, the false edges are generated. There are two ways to solve the contradiction. One is to use different filters and different smoothing parameters to preprocess images, then to do edge detection, and finally to analyze the results for synthesizing the image edges. The alternative way is to use different smoothing methods for different parts of an image based on different noise and edge types as described in the following.
The Gaussian filter of a signal ( ) can be represented as: where ( ) The expression omits the higher order derivatives as: The error of the pretreatment is defined as: Therefore, The latter method is used to improve: assuming that the available signal ( ) Generally, the human visual system is more sensitive to edges and noise of flat portions. When the second derivative of edges is fairly large, 2 σ must be as small as possible. When the second derivative of the flat portions is fairly small, 2 σ must be as large as possible. The equation (16) is used instead of equation (10), Here, ( ) Where, k is a scale factor. When x is in the image edges, ( ) σ is changed with the change of ( y x, ).
According to the features, a different 2 σ in a different position is applied to realize the purpose. Figure 3 shows The threshold determination is often made according to the principle of the suitable optimization in a thresholding algorithm. An image is divided into two classes: objects (o) and background (b), and each of the classes are supposed to have a normal distribution, in which the parameters can be obtained from the histogram of an original image: Where, t is a threshold, g is a gray scale value (0-255). The variances in the two categories are estimated respectively as follows: Where, the priori probability of the object class is The difference between the objects and the background in an image is measured by using the cross-entropy, combined with the Bayesian judgment [32] ; and the difference might be represented The difference between classes is: Replacing pixel gray scale s with the gray value g for simplifying the calculation is for replacing a probability distribution with a gray scale histogram. The equation (25) can be re-written as: Where, L is a grey value for upper bound, and T is a grey value threshold.
To obtain the optimal threshold T* based on the maximum cross-entropy between classes, it can be achieved through a searching operation: ( : ; *) max ( : ; ) The edge detection and post processing results are displayed in Fig.1, and the three sequences of processed images are illustrated in Figs.4-6.

Non-functional movement parameters of ventricular contours
The heart rate of a normal adult is about 75 times per minute and the average cardiac cycle is about 0.8 second. In a cardiac cycle, the systole is 0.1 second and the diastole is 0.7 second for the atrial. And for the ventricle, the systole is 0.3 second and the distole is 0.5 second. In each period of the cardiac systole and diastole, both the direction and the force for the ventricle are various, so that the overall displacement of each frame and the deflection in the B-mode irradiation plane are different. As it can be seen from the contour map, the obtained contour is of irregular curved, which is called elliptic. It contains the higher harmonic components, which is caused by the movement of the ventricle and the existence of the papillary muscle with the noise in the B-mode ultrasound images. With the systole and diastole of the heart, the size of the elliptic varies constantly. The changes in the direction of the shaft diameter reflect the functional movements of the heart -the systole and diastole. And the overall displacement and deflection reflect the non-functional movements of the heart. As the Fourier descriptors and the centroid principal axes algorithm [33][34] can both describe the overall movement of the object and ignore the details (high frequency components), the Fourier descriptors and the classical centroid principal axes algorithm [35][36][37] are adopted in analyzing and calculating the contour of the left ventricle to obtain the overall movement data.
Suppose the abscissa of an image is in the real axis, the ordinates of an image are the imaginary axes and the ventricular contour is composed of N points. Any point in the ventricular contour is around a circle curve, which can be used to get a one-dimensional complex sequence: The Discrete Fourier Transform (DFT) of ) k ( s is: Specifically, u is the frequency, As the high-frequency components of the Fourier transform correspond to the details of the contour and the low-frequency components of the Fourier transform related to the shape of the contour, we can use the first M Fourier coefficients based on the low-frequency components to rebuild the boundary to get an approximate contour with a closed curve: In light of the reconstruction experiment of the contour, (33) If the translation of the contour in the plane is (Δx, Δy), then When the size of the ventricular contour varies with the cardiac systole and diastole, the contour curve can be expressed as is that the Fourier coefficients are magnified C times. Thus, the impact of the ventricular size change caused by the heart systole and diastole can be eliminated by the normalized Fourier coefficients.
The coordinate origin is positioned at the centroid of the ventricular area. The contour curve turned around by angle θ in a space domain can be expressed as . Its Fourier transform is: The Fourier coefficients after rotation are equal to that of the origin Fourier coefficients times exp(jθ). To be precise, the phases of different frequency components in the Fourier coefficients are increased by an angle θ. So the rotation of the boundary can be detected by the phase changes of the Fourier coefficients.
Eliminating the DC component by the equation (26), the ratio of the first M Fourier coefficients from contours of two continuous images is: In Equation (37), C is the amplitude ratio of the Fourier coefficients, which respects the scale of the boundary. And ) u ( ϕ ∆ represent the phase changes of the Fourier coefficients with different frequencies. The phase changes of the Fourier coefficients for the boundaries of two continuous images can be got by the Equations (25), (26) and (27).
The lines in the parameter z , θ space corresponding to the collinear points in  (29) can be obtained, the rotation angle of the boundary curve can be determined in a plane [36][37] .
The translation of a contour in a plane can be detected by the DC component changes in the

Experiments and discussion
In the experimental, all the images are preprocessed by a new kernel of Fractional differential, then they are segmented by the improved Canny edge detector, finally they are analyzed by using Fourier descriptor. The details are as follows.

Video processing and analysis
The comparative data are corresponding to the short axis B-images for the left ventricular, the numbers of which are in 11-20, 22, 24, 26, 28, 30, 32, 38, 40, 42, 44, 46,      In Fig.7(b), when the ventricle contracts to its smallest size, the y-axis coordinate in the contour center is decreased significantly, which is because the lower parts of the ventricular contour are contracted very narrowly (see Figs. [4][5][6] so that the y-axis coordinate used to calculate the contour centroid is decreased (the centroid position is shift up). Using the Fourier descriptor to reconstruct the boundary curve, the center coordinates of the contour are determined as shown in Fig.7(a), which reduces the effects of the heart structure like the papillary muscle and brings the ventricular contour close to the actual case.  Fig.7, the centroid of the ventricular contour shifts slightly in the process of the ventricular systole and diastole, that is, the ventricular wall does the same shift. Though the translation of each frame is not more than two pixels in distance, where one pixel is equal to 0.44 mm, the cumulative displacement of multiple frames still has a significant impact on the sampling waveform. So it is necessary to correct the translation motion of the heart.
In order to compare the deflection angles of the principal axes, the principal direction changing in the ventricular contours are calculated in Table 1. Where, y a is the principal angle calculated by the Fourier descriptors, b y is the principal angle calculated by the centroid principal axes algorithm. As it can be seen from Table 1, the principal direction of the contour is changed a little during the ventricular systole and diastole. When the principal angle increases in the systole, the angle in the diastole will decrease. When the ventricle contracts to its smallest size, the principal axis is almost vertical. The variation rule and the numerical value of the principal axis calculated by two ways are almost the same.
From the above results, both methods can be adopted to isolate the functional movements and the non-functional movements of the ventricle. The non-functional movement parameters of the ventricular contours calculated by the Fourier descriptors are close to that of the actual movement of the heart.

Correction of omnidirectional M-mode echocardiography
The (41)   The method is applied to construct the omnidirectional M-mode echocardiography. Figure 7 shows the omnidirectional M-mode echocardiography sampled by the fixed direction line (the direction line of the non-functional movements is not removed) and the tracking sampling line (the direction line of the non-functional movements is removed). Figure 9 also shows the comparison results in four sampling positions. In each group, the left picture is the location of the sampling line, and the right one is the sampling result. From Fig.9, the sampling line basically is for tracking the specific structure after the correction of the non-functional movements. In Fig.9(a), the movement in the overall cardiac cycle can be obtained. In Fig.9(a, b, d), the obtained omnidirectional M-mode echocardiography cycle is more apparent, and the waveform has the strong repeat-ability. In Fig.9(b), the obtained contour is clear and there is no specific point lost or wandered basically. In Fig.9(b, c), the Echocardiography is complete, and no sampling structure is lost or other organizational structures appear. The crest of the Echocardiography is sharper, which reflects the movement when the heart is impulse and the movement is not smoothed by the overall movement.

Conclusions
From our experiments, we proved that it is feasible to extract the left ventricular contours of the short axis with the Fractional differential and improved Canny edge detection algorithms. This technique is supported by the gray scale search algorithm based on radial, and we demonstrate the accuracy of our final results. For correcting the non-functional movements, we proposed to use the studied motion tracking method based on the ventricular contours. By correcting the non-functional movements, the sampling direction lines can basically be used to track the movement of a certain part in the heart of the same structure at different times, and the functional movements in the specific parts of the cardiac structures can be reflected accurately. The non-functional movements in the B ultrasonic plane can be represented by the overall displacement and the principal axis deflection of the ventricular contours. The Fourier descriptors are adopted to obtain the principal translation and rotation of a series of object contour curves accurately. However, this technique tends to exhaust huge computational expenses.