Application of Time-Frequency Domain Transform to Three-Dimensional Interpolation of Medical Images

2.1. The method of time-frequency domain transform

Time domain and frequency domain are the basic properties of the signal. They can be used in different ways and aspects to analyze the signal. The method of time-frequency domain transform plays an important role in image processing. In this section, we describe some common methods: Fourier transform, DCT, and wavelet transform.

Fourier transform is a signal analysis method; it can analyze and synthesize the signal components. Fourier transform has been widely used in physics, mathematics, optics, image processing, and other fields. In digital image processing, two-dimensional (2D) Fourier transform is used to transform the 2D gray image f (x, y) to the frequency domain F(u, v). The center of the spectrum is the low-frequency part of the original image, and it reflects the overview of the image. The surrounding part of the spectrum is the low-frequency part of the original image, and it has refection on the sharp edges of the image.

The formula of the 2D Fourier transform is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} F ( u , v ) = \int \nolimits _{ - \infty }^{ + \infty } \int \nolimits _{ - \infty }^{ + \infty } f ( x , y ) {e^{ - j2 \pi ( ux + vy ) }}dxdy. \tag{1} \end{align*} \end{document}

DCT is a special form of Fourier transform. It can transform a set of time-domain data into frequency-domain data. It is considered the best method for the transformation of speech and image signals.

Wavelet analysis is a new field with rapid development in mathematics. This theory is profound and has a wide range of applications. Compared with Fourier transform, it is a kind of flexible time-frequency domain transform. The wavelet transform breaks through the limitation of the traditional Fourier transform; it can make the multi-scale analysis with function and signal through expansion, translation, and other computing functions. In the low-frequency part, signal has high frequency resolution and low time resolution; in the high-frequency part, signal has high time resolution and low frequency resolution.

Wavelet transform is based on Fourier transform. The most important change is: The infinite trigonometric function is replaced by finite decay wavelet bases. The form of wavelet bases is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} WT ( a , \tau ) = \frac { 1 } { { \sqrt { \rm { a } } } } \int \nolimits _ { - \infty } ^ \infty f ( t ) \times \psi ( \frac { { t - \tau } } { a } ) dt. \tag { 2 } \end{align*} \end{document}

In the field of digital images, the main application is the discrete 2D wavelet transform. It is a direct derivation of one-dimensional wavelet transform and its general form is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} W{T_f} ( j , \bar n ) = { \left\vert {{A_0}} \right\vert ^{ - j}} \iint {f ( x , y ) \bullet \psi [ a_{11}^{ ( j ) }} x + a_{12}^{ ( j ) }y - {n_1} , a_{21}^{ ( j ) }x + a_{22}^{ ( j ) }y - {n_2} ] dxdy. \tag{3} \end{align*} \end{document}

When the wavelet transform is used in 2D digital images (Mallat, 1990), it can be regarded as a group of band-pass filters for filtering the original image, and it can get sub-graphs by high- and low-pass filters (Tapia et al., 2004). As shown in Figure 1 (to show more clearly, we do a special treatment for each sub-image), LL is the low-frequency sub-image, and LH, HL, and HH are the high-frequency sub-images.

OPEN IN VIEWER

FIG. 1.

Schematic diagram of a wavelet decomposition.

2.2. Application of time-frequency domain transform in 3D interpolation

The time-frequency domain transform is applied in the 3D interpolation. The main change is that the interpolation of time-domain images turns into the interpolation of frequency-domain images. Figure 2 shows the general methods of 3D interpolation in the time-frequency domain. The left flow in Figure 2 shows the method of 3D interpolation to the time-domain image when there is no time-frequency domain transform; the right flow shows the method of 3D interpolation after frequency domain transform. Before interpolation, the image is transformed from the time domain to the frequency domain, and then, 3D interpolation is performed in the frequency domain. Finally, the interpolation image is obtained by inverse time-frequency domain transform.

OPEN IN VIEWER

FIG. 2.

The general method of 3D interpolation in time-frequency domain. 3D, three-dimensional.

3.1. Calculate the position of candidate corresponding points

First, we select the fixed window width W. To ensure the symmetry of the window, the singular value is set as the window width, for example, 3 × 3, 5 × 5, 7 × 7, and so on. Second, as shown in Figure 4, the distance between adjacent layers is d, and we need to establish the gray value of point w(x′,y′,z′) on the middle layer. We can select a point w(x′,y′,z_k) on the layer Z_k, and connect this point and point w(x′,y′,z′) on the middle layer point with the extension line. This extension line and layer Z_k+1 intersect at a point w(x′,y′,z_k+1). Point w(x′,y′,z_k) and point w(x′,y′,z_k+1) constitute the corresponding point pair. With each point in the layer Z_k, we can find w² pairs of candidate corresponding points.

OPEN IN VIEWER

FIG. 4.

Position of candidate corresponding points of medical image.

3.2. Select the best correspondence point pair

We use a suitable matching criterion to judge candidate corresponding point pairs. The best correspondence point pair is selected by calculating the matching degree. After analysis and improvement, we use the following rules to select the best correspondence point pair:

(1) Gray value f (w_i,j) of the point w_i,j;

Set up the gray value of point w_i,j as f (w_i,j), the Sobel gradient as D_i,j, and the gray gradient direction as θ_i,j. The calculation formula is as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \Delta xf \left( {{w_{i , j}}} \right) = f \left( {{w_{i - }}_{1 , j + 1}} \right) + 2f \left( {{w_{i , j + }}_1} \right) + f \left( {{w_{i + }}_{1 , j + 1}} \right) - f \left( {{w_{i - }}_{1 , j - 1}} \right) - 2f \left( {{w_{i , j - }}_1} \right) - f \left( {{w_{i + }}_{1 , j - 1}} \right) \tag{4} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \Delta xf{ \rm{ ( }}{w_{i , j}}{ \rm{ ) }} = f{ \rm{ ( }}{w_{i - 1 , j + 1}}{ \rm{ ) }} + 2 f{ \rm{ ( }}{w_{i , j + 1}}{ \rm{ ) }} + f{ \rm{ ( }}{w_{i + 1 , j + 1}}{ \rm{ ) }} - f{ \rm{ ( }}{w_{i - 1 , j - 1}}{ \rm{ ) }} - 2f{ \rm{ ( }}{w_{i , j - 1}}{ \rm{ ) }} - f{ \rm{ ( }}{w_{i + 1 , j - 1}}{ \rm{ ) }} \tag{5} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {D_{i , j}} = \sqrt { \Delta xf{{{ \rm{ ( }}{w_{i , j}}{ \rm{ ) }}}^2} + \Delta yf{{{ \rm{ ( }}{w_{i , j}}{ \rm{ ) }}}^2}} \tag{6} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \theta _ { i , j } } = arctg { \rm { ( } } { \frac { \Delta yf { \rm { ( } } { w_ { i , j } } { \rm { ) } } } { \Delta xf { \rm { ( } } { w_ { i , j } } { \rm { ) } } } } { \rm { ) } } \times \frac { { 180 } } { \pi } . \tag { 7 } \end{align*} \end{document}

At the same time, considering that the matching degree between corresponding point pairs is inversely proportional to the distance, the distance exponential distribution function is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} L { \rm { ( } } w { \rm { ) } } = \exp \ { \rm { ( } } - \frac { { { { { \rm { ( } } { x_ { i ^\prime , j ^\prime } } - { x_ { i , j } } { \rm { ) } } } ^2 } + { { { \rm { ( } } { y_ { i ^\prime , j ^\prime } } - { y_ { i , j } } { \rm { ) } } } ^2 } } } { 2 } { \rm { ) } } . \tag { 8 } \end{align*} \end{document}

Comprehensively analyze the characteristics of the correspondence point pair, and add the weight λ₁ of gray value, the weight λ₂ of Sobel gradient, and the weight λ₃ of gradient direction. The matching criterion formula is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm R } { \rm ( i , j , i ^\prime , j ^\prime ) } = { \frac { { \lambda _1 } \left\vert { { { \rm f } _ { i , j } } - { { \rm f } _ { i ^\prime , j ^\prime } } } \right\vert + { \lambda _2 } \left\vert { { { \rm D } _ { i , j } } - { { \rm D } _ { i ^\prime , j^ \prime } } } \right\vert + { \lambda _1 } \left\vert { { \theta _ { i , j } } - { \theta _ { i^ \prime , j ^\prime } } } \right\vert } { { \rm L ( w ) } } } . \tag { 9 } \end{align*} \end{document}

When the minimum value of R(i,j,i′,j′) is obtained, the matching degree of the candidate matching point is the highest, and the best matching point is obtained. The value is 0 when the two medical images in the region are of the same pixel value.

3.3. Calculate the gray value of the interpolation point

After getting the best corresponding points w_i,j and w_i′,j′, we can obtain the gray value of the target interpolation point through their gray value. It is assumed that the variation between the target interpolation point and the best corresponding point pair is linear. Therefore, the gray value f (w_x,y) of point w_x,y is obtained by the linear interpolation formula. The distance between w_x,y and Z_k is set as d₁, and the distance between w_x,y and Z_k+1 is d₂. The interpolation formula is as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} f { \rm { ( } } { w_ { x , y } } { \rm { ) } } = \frac { 1 } { { { d_1 } + { d_2 } } } { \rm { ( } } f { \rm { ( } } { w_ { i , j } } { \rm { ) } } \times d { } _2 + f { \rm { ( } } { w_ { i^ \prime , j ^\prime } } { \rm { ) } } \times { d_1 } { \rm { ) } } . \tag { 10 } \end{align*} \end{document}

We can get the final image through the interpolation of all points in the target image by using the method described earlier.

4.1. Comparison of interpolation result by different transform

In this article, we choose three different types of time-frequency domain transform: DCT, Fourier transform, and wavelet transform. Due to the lack of general methods for 3D interpolation in the frequency domain, we select the linear interpolation method. We do linear interpolation after the transformation of the graph, and the interpolation image is obtained by the inverse transformation. These images are compared with images performing interpolation directly without transformation. We analyze the values of the four indicators mentioned earlier. Table 1 and Figure 5 show the comparison among various interpolation methods.

OPEN IN VIEWER

FIG. 5.

Comparison chart of interpolation result by different transforms.

As shown by the experimental results mentioned earlier, interpolation without transformation has a higher running speed, but its effect is the worst. If the 3D interpolation method is not optimized and improved, σ(I_d), β(I_d), and λ(I_d) all have corresponding reductions and improved effects through a variety of different time-frequency domain transformations. The effects do not have many differences among various transformations.

This shows that the time-frequency domain transformation improves the quality of 3D interpolation, and in the frequency domain, there are a variety of image processing methods that can further enhance the quality of interpolation.

Wavelet transform decomposition will obtain four sub-images, which is different from other transformations. We can perform different operations on each sub-image. Low-frequency sub-images roughly have the appearances of original images. Many other interpolation methods can be used in 3D interpolation. Therefore, although the wavelet transform takes more time, its advantage is obvious.

4.2. Comparison of interpolation result by different wavelet transforms

In this article, we selected bior3.7, haar, coif1, sym2, dmey, and rbio3.7 wavelet to decompose and reconstruct tomographic images. With the proposed method, we experiment on two tomographic images. The parameters are set as: λ₁ = 1, λ₂ = 2, λ₃ = 1.

We analyze and compare the values of the four indicators mentioned earlier. Table 2 and Figure 6 are given to compare the interpolation effects with different wavelets.

OPEN IN VIEWER

FIG. 6.

Comparison chart of interpolation result by different wavelet transforms.

From the experimental results mentioned earlier, the execution time of the method using wavelet bior3.7 has the best results on mean-square deviation, gray variance values ranging from the total number of points, and gray value difference of the absolute value of the parameter, although it is not the shortest. Thus, we choose bior3.7 as the wavelet to perform wavelet decomposition and reconstruction.

4.3. Comparison of the proposed method and traditional methods

Experiments are executed with linear interpolation, traditional matching interpolation, and the method proposed in this article on five groups of tomographic images. The parameters of the interpolation method proposed in this article are set as λ₁ = 1, λ₂ = 2, and λ₃ = 1, and bior3.7 wavelet is used. Table 3 and Figure 7 compare our proposed method with the others. It can be seen that in the same experimental image and under the same experimental conditions, the average variance of the proposed method, the total number of pixels, and the absolute value of all gray values are significantly lower than those in the traditional methods. But it takes more time than linear interpolation and traditional matching interpolation.

OPEN IN VIEWER

FIG. 7.

Comparison of interpolation result by different methods.

Figure 8 depicts the results of the fifth group experiments; Figure 8a–c indicates the original images; and Figure 8d–f depicts the effects after various methods of interpolation. Compared with the traditional method, the proposed method has a clear outline of the image that is the closest to the original images.

OPEN IN VIEWER

FIG. 8.

The interpolation result by different methods.