Introduction

Many heart diseases alter form and dynamics of the surface of the endocardium. Reconstructing the endocardium is therefore an important task of medical image processing. Ultrasound imaging is a low-invasive imaging modality widely available due to its low cost. Technical advantages are the high temporal and spatial resolution. A modality for acquisition of volumetric data sets is

Figure 1. Imaging cone in transesophageal ultrasound

fourdimesional transesophageal ultrasound. This modality uses a transducer probe that has to be swallowed by the patient to image the heart from the esophagus. At the tip of the probe a transducer array is positioned that is capable of acquiring two-dimensional sequences. Rotating the tip of the probe the imaging plane is rotated so that a cone-shaped volume can be acquired (Fig.1). For each imaging plane an image sequence of one heart beat is acquired. Resorting the images one obtains for each moment of the heartbeat one volumetric dataset.

The acquisition of the whole four-dimensional dataset takes several minutes. One problem of the imaging modality are motion artifacts that occur due to unintentional movements of the patient during that long examination procedure. They lead to displacements and rotations of the imaging plane relative to their ideal positions.

Methods

Active surfaces

The analysis of endocard motion is based on the reconstruction of a parametrized endocardial surface X(u,v) for each moment of the heartbeat. Cohen et. al. [1] and McInerney et.al. [3] formulated surface reconstruction as an energy-minimizing problem. An energy function E_tot=E_i+E_f  for a surface is set up consisting of an internal energy E_i enforcing smoothness and a feature energy E_f  enforcing similarity between the surface points and corresponding feature points computed from the image. The surface to be reconstructed relates to a local minimum of the total energy.
 

(1)

(2)

X(u_ij,v_ij) are points of the surface, X_fij are the corresponding feature points each numbered by i and j.

In order to solve the problem numerically the surface has to be approximated by a superposition of a finite number of functions. In the case of a finite element approximation these are the piecewise polynomial formfunctions [2]. The expansion coefficients are called knotvariables. So the approximated surface can be written as the product of the rowvector of formfunctions and a columnvector of knotvariables.

(3)

We follow the notation of [2] where columnvectors are symbolized with pointed brackets, rowvectors are symbolized with braces and matrices are symbolized with square brackets. The finite-element approximation casts the function of the internal energy into a quadratic function of the knotvariables:

(4)

At a local minimum of the total energy the partial derivatives after the knotvariables vanish:

(5)

where

(6)

is the vector of so called generalized feature forces. The generalized feature forces depend on the position of the surface given by the knotvariables through eq.3. Eq.5 is a precondition for a local minimum of the total energy. A way of finding the local minimum of the total energy is introducing an iteration time and descending the gradient of the total energy according to:

(7)

This evolution equation converges to a minimum of the total energy. The gradient descent algorithm needs to be given an initial condition to start from. The evolution equation is solved by descretizing eq.7 with finite differences in the time domain and integrating numerically:

(8)

The iteration can be stopped if the change between consecutive time steps falls below a given threshold. The solution of eq.8 corresponds to the evolution of a damped elastic membrane under the influence of a space dependent image force computed from the image according to eq.6.

An integrated pseudomechanical model

The reconstruction of the endocardial surface requires the correction of the above mentioned motion artifacts. Correcting the motion artifacts on the basis of a correlation of neighbouring images is problematic, because speckle noise is uncorrelated between neighbouring images.

Correction of the motion artifacts and surface reconstruction are closely related problems. On the one hand correction of the motion parameters influences the surface to be reconstructed. On the other hand a reconstructed surface can serve as a reference for computing the displacement of the imaging planes [4], [5].

The problem to be solved can be put in the following words:
Compute
- a surface and
- a set of displacement parameters
so that:
- the surface is possibly smooth and
- the surface is possibly similar to the feature points.

These demands shall be put into an energy-minimising framework.
First the displacement parameters of the imaging planes must be specified. Each imaging plane has six degrees of freedom of a solid body. Dr_i and Dz_i are translation coordinates of image plane i in a cylindrical coordinate system. Dt_i is a translation perpendicular to the image plane. Da_i, Db_i and Dj_i are rotations around the r-, z- respectively the t-axis. These parameters are combined to a parameter-vector for each imaging plane. Feature points {r_ij,z_ij} are computed in each imaging plane, i being an index for the imaging plane. The index j indicates the feature points computed in the imaging plane. The location of the feature points in space depends on the location of the feature points in the imaging plane and the displacement parameters of the respective imaging plane:
 

(9)

where  ,   and   are rotation matrices. The requirement of similarity between the reconstructed surface and feature points in space can be formalized in a feature energy:

(10)

that is depending on both, the degrees of freedom {q} of the surface and the degrees of freedom {v} of the displacement parameters of all imaging planes. The task of computing a surface and a set of displacement parameters so that the surface is possibly smooth and possibly near to the feature points can be formulated as an energy minimization problem in the following way: Find a local minimum of a total energy consisting of an internal energy depending on the surface and a feature energy depending on both the surface and the displacement parameters:

(11)

The internal energy is the internal energy of an active surface. At a local minimum of the total energy the derivatives after the degrees of freedom of the surface {q} and the derivatives after the displacement parameters {v}must vanish. The first leads to a condition of balance similar to (eq.5). The letter leads a set of equations:

(12)

If the rotation axis of the ventricle and the rotation axis of the imaging cone are approximately identical, the derivatives after a_i, b_i and t_i are negligible so that one can confine oneself to the computation of displacements Dr_i and Dz_i and rotation Dj_i i.e. motion in the imaging plane. Eq.12 is a condition of balance of the sum of the feature forces of one imaging plane that retroact on the imaging plane. The term containing the derivatives after Dj_i corresponds to the torque of the feature forces. The local minimum of the total energy is again found descending the gradient of the total energy according to:

(13)

(14)

Eq.13 and Eq.14 are coupled because feature forces {f_q} and {f_v} depend on both the position of the imaging plane and an the displacement parameters. This energy minimising model can be interpreted as a pseudomechanical model of an elastic membrane coupled to freely movable imaging planes by feature forces computed in the images (Fig.2). In this picture the gradient descent algorithm is a damped movement of the coupled system to its equilibrium position.
 

Figure 2. Pseudomechanical model

Results

The Algorithm was applied to artificial data sets with given displacement parameters and to real data sets. Applying the algorithm to artificial data sets both the surface and the displacement parameters were determined correctly. In both cases displacement parameters were determined during the first steps of the iteration. Iteration was stopped when the difference of the surface position was below a threshold of 1 pixel. Convergence was reached after approximately 40 iteration steps. Fig.3 is linked with a sequence showing the evolution of the intersecting contour during the iteration. The rigid body movement of the intersecting contour relates to a change of the displacement parameters of the slice plane during the iteration. Fig. 4. is linked with a sequence showing the evolution of the surface during the iteration.
 
 

 

Figure 3. Temporal evolution of the intersecting contour

 

 

Figure 4. Temporal evolution of the surface

 

Conclusions

This paper presents an energy-minimizing algorithm that integrates the tasks of surface detection and correction of motion artifacts in fourdimensional transesophageal ultrasound sequences. The problem arises whenever the raw data for a surface reconstruction scheme is acquired tomographically. The regularization capabilities of the algorithm can be improved if additional temporal smoothness constraints are taken into account. Another way of enhancing stability is starting from the assumption that unintentional movements occur on a slower time scale than heart motion itself. One displacement vector for one heart beat i.e. one displacement parameter per slice plane would then be sufficient. The reduced number of displacement parameters would contribute to additional stability of the algorithm.
 

References