Trajectory optimization for the planning of percutaneous radiofrequency ablation of hepatic tumors

Abstract

Radiofrequency ablation is increasingly used in the treatment of hepatic tumors. Planning the percutaneous intervention is essential and particularly difficult. In this paper, we focus on an automated computation of optimal needle insertion in computer-assisted surgery with 3D visualization. First, we review our method which delineates on the skin of a virtual patient the candidate zones for needle insertion, i.e., those which allow safe access to the tumor. In each case, we look for the trajectory that minimizes the volume of burnt tissue. Secondly, we introduce a quasi-exhaustive method that allies sampling and certified minimization to form a strong argument for the accuracy of our results. We also compare results of applying both methods on 7 representative reconstructed patient cases.

• Radiofrequency ablation
• surgery planning
• optimization

Introduction

The treatment of hepatic cancers is nowadays mainly surgical. However, there are still many contra-indications to liver resection (localization and size of tumor, inoperable patient, etc.). New minimally invasive techniques have been developed in the past few years to enable the treatment of patients for whom open surgery is not feasible. Among these techniques are chemo-embolization, cryotherapy, ethanol injection, and radiofrequency ablation (RFA). The latter appears to be one of the most secure and easily predictable methods, and has seen increasing use.

RFA achieves the destruction of the tumor by hyperthermia, using a high-frequency alternating current delivered via an electrode inserted into the tumor tissues. Planning of such an operation is essential to its success. For percutaneous RFA, the physician must choose an accurate strategy for the insertion of the probe, allowing him to destroy the entire tumor while preserving surrounding vital organs. This planning step is a delicate procedure that is generally performed using only 2D CT slices of the patient, or with 2D ultrasound views.

We developed a specific planning tool called RF-Sim, aimed at guiding the choice of the physician by providing as much information as possible, and thereby reducing the risks of complications and tumor recurrence. This program is based on a 3D reconstruction of the patient's organs, offering an intuitive visualization of his anatomy. RF-Sim is a complete planning tool that integrates both visualization of the lesions induced by the treatment, allowing experimentation with various strategies, and an automatic planning algorithm that proposes optimal trajectories for the needle. A further noteworthy aspect of this tool is the realistic simulation of the surgical gesture by means of a haptic feedback device and virtual reality [1].

In this paper, we will focus exclusively on treatment planning, i.e., computation of an optimal trajectory that avoids vital structures. In the following sections, we present a brief overview of the state of the art and recall the existing functionalities of RF-Sim. We then detail the quasi-exhaustive method introduced to strengthen the arguments for the accuracy of our results. Finally, we discuss the results obtained and possible future extensions of this program.

State of the art

The medical literature concerning RFA is rich [2–5]. These studies detail the principles of the approach, contra-indications [6], characteristics of the different types of devices [7], and possible complications of the operation [8]. They describe common causes of failure [9], and factors that may be modified to improve the effects of treatment [10]. The influence of the human factor has also been evaluated. Poon et al. [11] showed that the chances of success are closely dependent on the experience of the surgeon. Moreover, Antoch et al. [12] underlined the importance of a volumetric visualization of the patient for adequate placement of the probe inside the tumor. The elaboration of a program integrating both 3D visualization and planning assistance should contribute to the training of surgeons and improve the success of such operations.

Numerous works concerning planning and simulation of surgical operations have been published. However, as RFA is a fairly recent technique, there are only a few published studies on the simulation and planning of this specific operation. References [13] and [14] present studies focusing on the simulation of the lesion induced by radiofrequency. By means of simulations based on Finite Element Modeling (FEM), modeling thermal exchanges during the operation, the lesion shape can be evaluated. However, within the framework of our program, such methods would have to be adapted to our constraint requiring real-time visualization of the lesion. Littmann et al. [15] described software for hepatic resection planning, including the possibility of simulating manually an RFA operation in the case of a non-resectable tumor. Simultaneous 2D/3D visualization and realistic simulation of the lesion are possible, but the software does not provide automatic computation of an optimal needle trajectory.

Adhami et al. [16] have presented an optimization method that positions different tools (endoscope and robot arms) in robotically assisted heart surgery. The challenges here are similar to those in RFA, in that the goal is to choose trajectories that avoid organs and minimize a chosen criterion. Their method evaluates a limited set of given trajectories and chooses suitable trajectories using an exhaustive method, while our goal is to take into account all possibilities.

Mention should also be made of research in the field of cryoablation, which achieves the destruction of tumors by freezing rather than heating. The concerns of both cryoablation and RFA are consequently similar in many respects. Butz et al. [17] have presented a planning tool for cryoablation that is adaptable for RFA, and which integrates the computation of an optimal needle insertion strategy. However, the proposed method is quite time-consuming and requires the prior definition of an insertion window, whereas we wish to acquire a very rapid and fully automatic proposition after a search over the entire abdomen.

The simulator: RF-Sim

The planning software that we are developing relies on a method of 3D reconstruction from 2D CT slices as described in reference [18]. It provides a 3D scene that is very easy to visualize and manipulate. Moreover, it is possible to introduce into the scene one or more RF probes, place them inside the tumor, and visualize the corresponding lesions. The simulation of cellular heating, the impedance variation, and the thermal exchanges within the tissues would allow determination of the exact shape of the lesions produced, but the methods proposed so far do not allow such an estimate to be performed quickly enough to enable the result to be displayed in real time. Mulier et al. [19] have described the properties of the lesions produced by different types of RF needles. The shape of the lesions is predictable and can be approximated in most cases to a spheroid, the size of which depends on the duration of exposure of the tissues to heat. In some cases, the proximity of large vessels cools the surrounding area and deforms the shape of the lesion. Therefore, in our simulator the lesion is approximated by a spheroid, and the “heat-sink effect” caused by the blood flow is represented as a “repulsion” of the vertices of the spheroid [20]. This model takes into account a quick and close approximation of the main characteristics of the RF burn, and allows interactive modification of the configuration of the probe while visualizing the resulting lesion in real time.

The program also integrates a functionality that computes an optimal position for the needle insertion after the user has chosen an initial position. The main criterion for this optimization is the volume of the resulting lesion. Basically, a needle placement will be considered as increasingly interesting when it damages as little healthy tissue as possible while still covering the entire tumor. Therefore, the purpose is to determine the minimum of a function that associates to each candidate needle position the volume of healthy tissue burnt. However, needle placements that cross organs cannot be considered as appropriate trajectories. To eliminate trajectories that would damage vital organs, the function to minimize is modified by adding a penalty volume. Among the classical minimization methods [21], the downhill simplex method appeared to be the most efficient, combining, at best, fast computation with accuracy.

Nevertheless, the result obtained was not entirely satisfactory. The resulting position was always of interest, but it was sometimes restricted to a convex zone containing the initial position from which the process was launched. The process often faced boundaries composed of parts of bones, vessels, or other organs meshes, and was unable to cross them. This was due partly to the optimization method, and partly to the function to minimize. The simplex method is very sensitive to local minima, and the way in which penalties were added to avoid organs produced a few of these and confined the search within the artificial boundary. To find the optimal trajectory efficiently, it was often necessary to try various initial positions. We therefore could not consider the process as being truly fully automatic, though that was one of our goals. In the following sections we describe our proposed solutions to this problem.

Candidate zones for the insertion of a needle

The first step in our approach consists of automatically determining the candidate zones of the skin, i.e., the zones in which it is possible to insert a needle and reach the tumor without encountering any vital organs or bones. We summarize the algorithm below; more details can be found in reference [22].

Simplified algorithm

In a first step, let us consider the tumor as a single point in space. (It will be seen later how our method is adapted to take into account the actual volume of the tumor.) If we place a virtual camera at the tumor location, each point of the skin that is visible to the camera and not hidden by an organ is a possible insertion point. If there is no obstacle between a point and the tumor, it means that the tumor can be reached without causing any damage. Still with a view to efficiency, we considered the skin in terms of a mesh. A triangle of the skin is considered as visible from the tumor, and therefore a candidate, if – and only if – all the points of this triangle are visible. All the points of a candidate triangle are thus possible insertion points.

The visibility of the triangles of the mesh is determined thanks to six renderings of the 3D scene, each corresponding to one face of a virtual cube surrounding the tumor. The vision volumes corresponding to the six views allow coverage of the whole scene (see Figure 1). That way, whatever triangle of the skin is considered, it is included in one of the views, making it possible to compute whether it is partially or totally hidden by an organ.

Figure 1. Covering of the scene with six renderings. [Color version available online.]

This computation is performed very quickly and takes advantage of the possibilities of the graphics card. On a Pentium 4 3.2 GHz with a Geforce 7800 GT and 2 Go RAM, the computation time is approximately 60 ms for scenes containing an average of 2000 triangles for the skin and 150,000 triangles for the other organs.

Algorithm taking tumor volume into account

We will now detail the modifications to the previous algorithm required to take into account the actual volume of the tumor. We could simply use the algorithm successively with each voxel constituting the tumor as a point of view. By doing this, we would obtain for each triangle an accessibility rank that would correspond to the ratio of voxels of the tumor that can be seen from the triangle. However, in order to be less time-consuming, we only consider voxels of the tumor's boundary, and among those we only render views corresponding to faces that are not adjacent to an internal voxel. A visibility ray coming from an internal voxel would inevitably pass through an external face, so it is not useful to consider them. Figures 2a and b show the faces considered for a voxel, and Figure 2c shows all considered faces for this case.

Figure 2. Computation of the 100% zone. (a, b) Visibility zone for 2 different voxels and corresponding views for voxel faces in bold (hatched areas are not considered for these specific voxels). (c) Resulting 100% zone once all boundary voxels are studied. (d) 100% zones on a slice of a virtual patient, reconstructed from a real case. [Color version available online.]

With this optimization, only an estimation of the accessibility rank is obtained; however, this does not have any undesirable consequences, as for the rest of our investigation it was decided to only consider zones from which it was possible to reach the whole tumor; we will call these 100% zones (Figure 3). Even if it were technically possible to reach parts of the tumor from some zones outside the 100% zones (i.e., from partially shaded zones), in practice this would induce trajectories that would be very risky, as they would pass extremely close to vital organs. In Figure 2c we emphasized the 100% zones that give full access to the tumor. In Table I, the surfaces of the 100% zones corresponding to each case are quite variable in size, but they are always sufficiently large to enable them to receive the needle insertion safely.

Figure 3. Example of a 100% zone. [Color version available online.]

Table I.  Results of computation of 100% zones.

Case #	No. of tumor external voxels	Execution time (s)	Size of 100% zone (cm^2)
1	220	17	234
2	325	21	52
3	401	28	118
4	417	31	165
5	484	32	176
6	1198	82	44
7	1220	81	254

It can be seen from Table I that the time taken to compute the 100% zones correlates strongly with the number of voxels of the tumor's boundary. In fact, this number determines the number of views that have to be computed. The candidate zones on the skin are computed in less than 2 minutes for skin meshes having an average of 2000 triangles, and for tumors having an average of approximately 600 voxels.

Targeted search within the connected components

Once the connected components of the 100% zones are computed, the minimization process has to be launched in each one of them in order to compare the minima and choose the optimal trajectory corresponding to the best result. To restrict the search within one of the connected components, we add a penalty to the volume returned by the function to minimize when the trajectory reaches the boundary, as in our earlier method, turning its drawback into an advantage.

In an initial approach, we launched one minimization per connected component of the 100% zone with a random initial position for each. As suspected, the results were satisfying for the convex zones, but in more complex zones we often fell into local minima. To avoid this problem, an extra initialization phase was added to the minimization process. The goal of this initialization phase is to perform a quick evaluation of the best strategy for needle insertion in the component, in order to initialize the position as close as possible to the minimum to be found. Thus, as the minimization process starts from the appropriate valley, it converges quickly to the accurate minimum, and is prevented from moving away to a wrong local minimum.

In choosing an initial position, the needle tip must first be positioned in the center of the bounding box of the tumor; then the axis is successively rotated to cross the barycenter of each triangle of the connected component, and the volume of the corresponding lesion is evaluated. We choose as an initial axis the one providing the smallest volume. The minimization process is then launched from this position. Figure 4 represents the different trajectories obtained in each of the connected components of the 100% zone.

Figure 4. Optimal trajectories in each connected component of the 100% zone. [Color version available online.]

Validation of the quality of the results

In order to estimate the reliability and accuracy of our results and to validate our fast method, we studied a way to compute an accurate approximation of the theoretical minimum of the volume function in a reasonable time. The computation of the theoretical minimum must undoubtedly avoid all local minima phenomena.

To this end, in a first step, we will show that, in our case, the problem of local minima of the minimization concerns only the orientation of the needle axis. This provides a certified minimization if the needle axis is fixed. Then, in a second step, we will explain how this property was used to compute a theoretical minimum with a quasi-exhaustive method allying a sampling of the axis with the certified minimization, and will then use this to evaluate the accuracy of our results.

No local minima concerning needle tip position

Assuming that the radiofrequency lesion has a spheroid shape and that the direction of the needle axis is fixed, the function that associates each needle tip position to the volume of the minimal lesion containing the tumor has only one local minimum, which is, of course, the absolute minimum. To prove this fact, we use a convexity argument.

The radiofrequency lesion is assimilated to a spheroid centered on the needle tip, its main axis corresponding to the needle axis. The ratio k between the short and long radius is fixed. Consider a Cartesian coordinate system (O, x, y, z) such that the x-axis corresponds to the fixed needle direction; the equation of the spheroid family that represents the lesion is then:where the center of the spheroid p(x_p, y_p, z_p) corresponds to the needle tip. Suppose that point p moves on a line l (not necessarily parallel to the needle axis) that is parametrically defined byIf we consider the function f_i that associates to each point p on l the radius of the minimal lesion centered in p including a point, P_i (x_i, \, y_i, \, z_i), we have:Then, the curve Γ_i that graphically represents the function , is a concave up parabola (see Figure 5). We conclude that there is a single value t for which g_i(t) is a minimum, i.e., there is a single position of the needle tip on l that corresponds to a minimal spheroidal lesion including P_i.

Figure 5. Graphic representation of the function g and its interpretation in space. Spheroids A, B and C in the diagram at left correspond to points A, B and C on the graph. [Color version available online.]

Spheroids that include all points P_i are also represented by a convex part of the plane because it is the intersection of all convex parts of the plane over parabolas Γ_i. We can then conclude that, for a fixed needle axis, there is a unique minimal spheroid placed on l and including all points P_i.

This can be extended by considering that the needle tip is moving in the entire space (the direction of the needle is already fixed). If we suppose that two different needle tip placements correspond to a local minimal volume, then these placements are two points that define a line. This is in contradiction with the above argument. Thus, we can conclude that, for a fixed needle axis, there is a unique needle tip position in the whole space for which the lesion including the entire tumor is a minimum.

Computation of a theoretical minimum

We have shown that, in our case, local minima affect only the choice of the needle axis. This information is useful for the computation of a theoretical minimum of our function.

One secure way to find this would be to use an exhaustive method by discretizing the whole parameter space. We tested this method on one patient using 3° increments for the needle axis and considering each tumor voxel as a possible needle tip position. Computation of this minimum took 12 hours and the precision was still low due to the chosen discretization steps.

In practice, a method based on discretization of the whole parameter space would have to sample five parameters (three for translation and two for rotation). Such a global sampling leads to a complexity containing an n⁵ factor, where n is the chosen number of samples in one dimension. This complexity induces unreasonably long computation times. For instance, if we want to double the precision of our sampling, we have to multiply the computation time by 2⁵ = 32, giving us a computation time of 384 hours (16 days) instead of 12 hours.

That is why we developed a hybrid method, to raise the precision of our theoretical minimum in a reasonable computation time. Our nearly exhaustive method mixes sampling on the two parameters representing the angles with minimization on the three parameters representing the position of the needle tip. All angles must be tried for the axis while a minimization method is satisfactory for the needle tip position. By comparing the volumes corresponding to each angle, we can obtain the best one.

Now that we are convinced that when angles are fixed minimization is not subjected to local minima anymore, we can say that the error obtained is mainly due to the precision of the sampling on the angles; the precision of the search will depend mainly on the discretization step for the needle axis. With this method, doubling the precision makes the process only four times slower.

Results

Table II summarizes the results obtained for 7 test scenes reconstructed from real patients from a radiology service at the Civil Hospital in Strasbourg, each of whom presented with a single hepatic tumor. The table shows the lesion volume resulting from our optimization method (“minimal volume”), i.e., the search for a minimum in each connected component as detailed above in the section entitled Targeted search within the connected components. The total computation time corresponding to this method is shown in the third column. Note that this computation took less than 1 minute for all test patients, with an average computation time of 4 seconds for each connected component.

Table II.  Results of minimization in each connected component of the 100% zone.

		Our method		Semi-exhaustive method
Case #	No. of connected components of 100% zone	Execution time (s)	Minimal volume (ml)	100% zone min. volume (ml)	Global min. volume (ml)
1	2	6.4	3.07	3.08	2.73
2	5	21.6	2.49	2.48	2.48
3	8	40.6	6.84	6.83	6.83
4	7	23.0	3.83	3.85	3.06
5	5	19.0	3.17	3.16	3.16
6	4	13.0	10.62	10.62	8.86
7	3	11.8	9.32	9.30	9.30

These results are compared with our computation of a theoretical minimum. We launched the computation first in a scene containing only the tumor, to obtain an idea of the volume of the minimal spheroid covering the whole tumor if there were no surrounding organ ("global min. volume"), and then only within the 100% zones ("100% zone min. volume"). Execution time was approximately 2 hours for the global minimum and between 30 min and 1 hour for the minimum restricted to the 100% zones, which is much faster than the exhaustive method while demonstrating better precision.

Comparing theoretical volumes obtained in 100% zones with global theoretical volumes, we observe that the search in 100% zones provides volumes having an average increase of +8% with respect to the global theoretical minimum volumes. This phenomenon is mostly inevitable, as it is due to the presence of organs obstructing the optimal trajectories. It can be mentioned that the restriction to the 100% zones adds only a small penalty to the obtained volume.

The comparison between results of our minimization and of the computation of a theoretical minimum shows the efficiency of our method, as the results are significantly close (average difference of +0.003 ml). In some cases, our method even exceeds the precision obtained by the discretization process used to compute the theoretical minimum. We conclude that our method quickly provides an estimation of the minimal volume of the lesion and of the corresponding needle placement in the 100% zone with good precision.

Discussion and future work

As mentioned earlier, in RF-Sim the computation of the optimal position of the needle is mainly based on a volume criterion. However, according to radiologists and surgeons, this is not the only criterion that must be taken into account for the planning of an RFA operation.

A trajectory considered to be optimal by the planning tool could be eliminated by the surgeon if he judges it to be unrealizable, too risky, difficult to reproduce, or simply not practical. Conversely, the physician, because of other, more interesting properties, could choose a trajectory that would lead to a less interesting volume. Various other criteria need to be integrated in the planning tool in order to provide more realistic assistance, and our tool is only a first step towards this goal.

For instance, we could integrate a degree of risk to needle trajectories that would depend on the distance between the trajectory and the vital organs, allowing a surgeon to choose to give priority to the minimization of risk instead of volume. Alternatively, a surgeon may wish to choose the most direct path for the needle, if minimizing the distance between the insertion point and the tumor is a major concern. However, when the tumor is close to Glisson's capsule (the liver boundary), the needle has to pass first through a minimum amount of healthy liver to avoid hemorrhage. Moreover, the surgeon could choose a totally different solution that has a quite analogous quality but is more compatible with his habits.

Consequently, as many as possible of these different criteria involved in the choice of a strategy have to be exhaustively determined in order to be taken into account during the process and thereby enhance our tool with an increased amount of information that is useful in the elaboration of the planning. These various forms of information come directly from the expertise of the surgeon, which is why our next step will be to collect from physicians all possible relevant information concerning the other decisive factors involved in the planning process, and to formalize them so that they will be exploitable in our tool. We also plan to develop methods to materialize intuitively the influence of these factors, either with visual and/or haptic representations.

Furthermore, we would like to develop a method that provides a realistic simulation of the lesion by simulating heat production and exchange within the organs. Our current estimation of the RF lesion is advantageous for the computation of the optimal trajectory because of its speed, but it would be valuable if a more precise simulation were provided, at least once a trajectory has been chosen.

Finally, it would be useful to reduce the total time necessary to plan an operation. Thanks to the method we presented here, RF-Sim automatically proposes a strategy in less than 5 minutes, but this is added to the time required by the 3D reconstruction from 2D CT slices. Moreover, in the future, we want our tool to be capable of handling many other criteria and planning multiple needle insertions, although we would like the entire process to take less than 15 min (our estimate of a reasonable duration). It will be necessary to keep this constraint in mind in our future developments.

Conclusion

We have presented a method for automatically and quickly computing an optimal trajectory for the insertion of a needle within the framework of percutaneous RFA. In a first step, our algorithm determines the zones over the skin that are candidates for insertion of the needle, then in a second phase it computes an optimal trajectory within each of these candidate zones, minimizing the volume of the induced lesion, and proposes the best one among them. The suggested trajectory is appropriate with respect to the criteria that we imposed–whole tumor burnt, minimum damage to healthy tissue, vital organs preserved–but this method only constitutes a first approach. In our future work we will adapt it to take into account various other factors that usually influence the decision of the surgeon during the elaboration of a planning strategy, allowing a more complete approach for each specific operation, and thereby increasing success rates.

Acknowledgment

We would like to thank Region Alsace for its financial support.