Electroporation-Based Treatment Planning for Deep-Seated Tumors Based on Automatic Liver Segmentation of MRI Images

Abstract

Electroporation is the phenomenon that occurs when a cell is exposed to a high electric field, which causes transient cell membrane permeabilization. A paramount electroporation-based application is electrochemotherapy, which is performed by delivering high-voltage electric pulses that enable the chemotherapeutic drug to more effectively destroy the tumor cells. Electrochemotherapy can be used for treating deep-seated metastases (e.g. in the liver, bone, brain, soft tissue) using variable-geometry long-needle electrodes. To treat deep-seated tumors, patient-specific treatment planning of the electroporation-based treatment is required. Treatment planning is based on generating a 3D model of the organ and target tissue subject to electroporation (i.e. tumor nodules). The generation of the 3D model is done by segmentation algorithms. We implemented and evaluated three automatic liver segmentation algorithms: region growing, adaptive threshold, and active contours (snakes). The algorithms were optimized using a seven-case dataset manually segmented by the radiologist as a training set, and finally validated using an additional four-case dataset that was previously not included in the optimization dataset. The presented results demonstrate that patient's medical images that were not included in the training set can be successfully segmented using our three algorithms. Besides electroporation-based treatments, these algorithms can be used in applications where automatic liver segmentation is required.

About the Authors:

Denis Pavliha

Affiliation: University of Ljubljana, Faculty of Electrical Engineering, Ljubljana, Slovenia

Maja M. Mušič

Affiliation: Institute of Oncology Ljubljana, Ljubljana, Slovenia

Gregor Serša

Affiliation: Institute of Oncology Ljubljana, Ljubljana, Slovenia

Damijan Miklavčič

* E-mail: [email protected]

Affiliation: University of Ljubljana, Faculty of Electrical Engineering, Ljubljana, Slovenia

Introduction

Electroporation is the phenomenon that occurs when a biological cell is exposed to an adequately high electric field, which results in the cell membrane becoming transiently permeabilized [1]. Electroporation is considered to be a universal method and platform technology since all types of cells (animal, plant, and microorganisms) can be electroporated [2]. A paramount electroporation-based application is electrochemotherapy [3], [4] which enhances chemotherapy outcome due to transient permeabilization of targeted cell membranes: because of the externally applied electric field, electroporation facilitates the chemotherapeutic drug diffusion through the plasma membrane into the cells, which would be otherwise hampered, because of the impaired or slow transport of the chemotherapeutics that are used in electrochemotherapy [5].

Electrochemotherapy is performed by high-voltage electric pulses delivery using applicators (i.e. electrodes) that are in contact with (or located near the) target tissue. Electrochemotherapy has already been introduced into clinical use for treating skin melanoma using plate or needle electrodes with a fixed geometry; the use of such electrodes imposes following predefined standard operating procedures for a successful treatment [6], [7]. Recently, however, electrochemotherapy has been introduced to clinical trials for treating deep-seated metastases in liver [8], bone [9], brain [10]–[12], and soft tissue [13]. Electrochemotherapy of deep-seated tumors imposes the use of variable-geometry long-needle electrodes introduced either percutaneously or during open-surgery [13]. Hence, only following the standard operating procedures cannot ensure success of the treatment, and patient-specific treatment planning is required for effective electrochemotherapy of deep-seated tumors [14].

Another important electroporation-based application is termed non-thermal irreversible electroporation (N-TIRE) and is used for tissue ablation performed using an externally applied electric field with electric field strengths and higher number of pulses than the values used for electrochemotherapy [15]. Nonetheless, the procedure is technologically very similar to electrochemotherapy and, also, requires patient-specific treatment planning [16].

To prepare a robust treatment plan for electroporation-based treatments, an anatomical model that is built from medical images (Magnetic Resonance Imaging – MRI) needs to be constructed first [17]. Construction of such a model is based on the acquisition of the patient's medical images and relies on processing the images in order to perform relevant-tissue extraction (i.e. image segmentation) [18]. Image segmentation, then, serves as the basis for generating a three-dimensional model consisting of the relevant healthy tissue (e.g. liver) and pathological tissue (i.e. tumors) [19]. Vessels may also be segmented and included into the model [20] since vessel positions have to be taken into account when defining electrodes' entry direction and relative positions. Then, a Finite-Element Model (FEM) is built and using the defined electrode parameters (number, dimensions, position), the distribution of the electrical field is calculated and optimized [21], [22] and finally presented to the attending physician.

In order to establish the concept of electroporation-based treatment planning, we follow radiotherapy treatment planning as the basis [23] using parallelisms and similarities between the planning procedures [14]. Since development of a user-friendly treatment planning would simplify electroporation-based preoperative procedures, we opt towards developing treatment planning software that will not require any prior engineering knowledge from its end-user (e.g. the attending physician). The whole treatment planning software needs to perform as automatically as possible, i.e. with minimum of interaction by the clinician, and the most challenging task is development and implementation of an automatic image segmentation algorithm. Within the clinical study of electrochemotherapy of colorectal metastases in the liver [8], we developed treatment planning procedure that includes liver segmentation. After implementing the segmentation algorithms and concluding the segmentation procedures, the latter were additionally modified using optimization results obtained using a training set of seven cases that were previously manually segmented by a radiologist. Finally, additional four cases were manually segmented by a radiologist and used for the final validation of the segmentation procedures.

In this study, we evaluated three different liver segmentation algorithms that can be used for electroporation-based treatment planning: region growing, adaptive threshold and active contours. Region growing was selected for evaluation because despite its simplicity (i.e. segments are included solely based on their intensities) this algorithm is robust and can provide good results if its basic drawbacks (e.g. oversegmentation due to leakage) are eliminated [24] using a postprocessor. Our implementation of the adaptive threshold algorithm was evaluated because this algorithm is based on a physical property, i.e. continuity of the tissue: segments on two neighboring slices are expected to be minimally different, which although being an intuitive solution which can be used as initialization of other segmentation methods [25], it proved to be accurate enough to be used as a standalone method for liver segmentation. Finally, the active contours algorithm [26] based on the Gradient Vector Flow (GVF) [27] was evaluated because of its insensitivity for intensity-based anomalies (e.g. inhomogeneity, or thin bonds connecting different segments such as the liver and e.g. kidneys) and possibility of influencing the movement of the contour by balancing the coefficient that influence attraction of the contours by the image or by the contour's inner properties. All three segmentation algorithms were optimized on a training set of seven cases, i.e. quantitatively assessed using real case data obtained from a radiologist. Finally, algorithms were validated using additional four real cases obtained from a radiologist, therefore accuracy of how their results are produced is known.

Methodology

2.1. Automatic Liver Segmentation

2.1.1. Importing DICOM Images.

The segmentation procedure begins with importing the patient's images into the treatment planning software. The latter is a MATLAB application, developed in MATLAB R2012a (Mathworks, Nantick, MA, USA) using the Image Processing Toolbox and Parallel Processing Toolbox. The procedure for loading images loads all DICOM (Digital Imaging and Communications in Medicine) [28] files from the user-defined folder, and reads their DICOM headers' SeriesNumber parameter in order to determine the number of different acquisition series present in the folder. Then, an image from every series is presented to the end-user (e.g. the attending physician); an image from every acquisition series is displayed and labeled using the original label that was stored at acquisition time and is read from the DICOM header as the string stored in the SeriesDescription parameter. Finally, the end-user determines which acquisition series will be used for planning of the electroporation-based treatment.

After that, all the images from the selected acquisition series are loaded and sorted according to their spatial location (i.e. according to their Z-index which can be read from the DICOM header as the SliceLocation parameter) using a common bubble-sorting algorithm. If all obtained Z-indexes after bubble-sorting are not evenly distributed, empty slices are inserted where the slices are detected as missing. However, since missing slices may indicate corruption of the patient's images collection, the software does not try to interpolate the missing slices but notifies the end-user instead.

Besides the image data, essential DICOM metadata is loaded; the metadata structure appended to the slices includes these parameters: Width, Height, SliceThickness, PixelSpacing, Modality, AcquisitionDate, BitsAllocated, and Volume of Interest (VOI) parameters WindowCenter (WC) and WindowWidth (WW) which are most important. Namely, WC and WW determine how source image data need to be interpreted when displayed; therefore, an initial sigmoid transformation using parameters WC and WW needs to be performed first for the image data to be displayed correctly. The latter is done within the preprocessing procedure.

2.1.2. Preprocessing.

For segmentation algorithms to perform without problems, the imported slices first need to be preprocessed. Preprocessing is a procedure which is executed on each slice separately; therefore, the procedure is non-recursive and can be run in parallel using multiple processors or processor cores. Since the preprocessing procedure comprises of several steps, the steps are marked for debugging and algorithm evaluation purposes by storing partial results (i.e. partially preprocessed slices) into separate layers, starting with the source slices (i.e. slices stored as DICOM data), which enables the developers of the algorithms to have a clear overview of the whole preprocessing procedure.

Interpretation of imported (i.e. source) slices is defined by the Volume-of-Interest (VOI) parameters (i.e. Window Center – WC, and Window Width – WW) that are stored as DICOM metadata. Since WC and WW differ from slice to slice, each slice first needs to be transformed from the imported data values (i.e. source layer) to the normal values (i.e. original layer) which are defined by the VOI parameters. The transformation can be performed using a sigmoid function as defined by the DICOM standard [29]; the transformation is described using Equation 1.(1)

As seen in Figure 1, such a transformation can be used for multiple purposes. When the WW parameter is small (), the sigmoid function (Figure 1B) changes into an approximation of a step function (Figure 1A) and can be used for thresholding, the WC parameter being the threshold value and the output value being Boolean with possible values (0, output_range). When the WW parameter is large (e.g. WW > output_range), the sigmoid function changes into an approximation of a linear function (Figure 1C) and can be used for linear transformations, the output value residing in the range (0, output_range) depending on the WC and WW parameter values.

[Figure omitted. See PDF.]

Figure 1. Sigmoid transformation function for preprocessing purposes.

Presented is an example sigmoid function (B), with approximations of the step function (A: WC  = 512, WW  = 1) and linear function (C: WC  = 512, WW  = 2048).

https://doi.org/10.1371/journal.pone.0069068.g001

After transforming the source slices into original slices using the transformation from Eq.1 and parameter values (WC, WW, and output_range) from the DICOM metadata, the slices are then de-biased. De-biasing is a procedure that removes intensity inhomogeneity [30] that is caused because the magnetic field in the area where the patient is positioned is not equally intense (i.e. the magnetic field is more homogeneous in the focal part of the device); a publicly available inhomogeneity correction algorithm was implemented for de-biasing [31]. Then, filtering of the slices, which is necessary for noise elimination, is performed by applying an average and a Gaussian blur filter (), both with window sizes of 3×3 pixels. Finally, another sigmoid transformation is applied to the slices using fixed VOI parameters (WC = 20000, WW = 100, output_range = 216) which ensures the intensity distribution of the slices is redistributed in the whole 16-bit range regardless of the source slices' range, and an adequate contrast which is dependent on the WW parameter. The fixed VOI values were selected based on our experience using real case data, and assure that the liver segment will have an appropriate intensity value range for the segmentation to be successful. After these procedures are applied, the resulting slices are stored as preprocessed slices (i.e. on a separate layer) and are ready for segmentation.

2.1.3. Referential Slice.

First, we define the referential slice as the slice in the patient's medical images collection with a high probability to include a large liver segment. We define ZREF as the index (i.e. spatial location) of the referential slice using Equation 2:(2)where ZMAX is the number of all the slices in the patient's medical images collection and the constant 0.65 was found empirically on real case data. Using Equation 2, we have a high probability of obtaining a referential slice with a liver segment that is morphologically similar to the liver segment shown in Figure 2B.

[Figure omitted. See PDF.]

Figure 2. An example three-dimensional liver object.

The presented object consists of 72 slices (ZMAX = 72, A) and includes slice ZREF = 47 as the referential slice (B).

https://doi.org/10.1371/journal.pone.0069068.g002

The identified referential slice serves as the beginning point of segmentation, i.e. the slice where segmentation is initiated, and its identification is independent of the chosen segmentation algorithm. Hence, there is a high probability that the referential slice will include a liver segment as the one shown in Figure 2B. For region growing, the referential slice is used as the slice that is presented to the end-user (e.g. attending physician) in order to place the initial seed on the slice; for adaptive threshold (and also for active contours which use adaptive threshold for contour initialization), the referential slice is used for comparing the dynamically thresholded slice to the presets that include similar liver segments, and marking the threshold value with the highest similarity to any of the presets as the initial liver segment.

2.1.4. Region Growing Segmentation Algorithm.

We first implemented an intensity-based segmentation algorithm known as the region growing algorithm. The latter determines whether voxels are part of the target region or not by comparing their intensities to the intensity of the initial seed. The initial seed is a voxel manually selected by the end-user (e.g. the attending physician) at the beginning of the procedure; in our case, the referential slice is presented to the end-user who is required to click on the liver segment where there are no internal liver structures such as vessels or tumor nodules. The pixel clicked then serves as the initial seed voxel.

The region growing segmentation algorithm works in three dimensions and evaluates the voxels that are scheduled into the queue. At the beginning, a single voxel is added to the queue, namely the initial seed. The algorithm examines the current voxel in the queue by comparing the intensity of every current voxel's neighbor to the intensity of the current voxel, as shown in Figure 3 where an array of 3×3×3 voxels is displayed and the middle voxel represents the current voxel (i.e. 26-connected neighbors).

[Figure omitted. See PDF.]

Figure 3. Representation of the initial seed voxel and its 26 neighbors.

https://doi.org/10.1371/journal.pone.0069068.g003

Since intensity of the target region varies in all three dimensions due to inhomogeneity, it is imperative to allow some intensity deviation when evaluating if the neighboring voxels belong to the target region. The allowed intensity deviation is defined using a threshold deviation value (e.g. setting the threshold deviation value to 0.20, which is the value we used, determines the intensities that are acceptable for inclusion into the target region; the determined intensities reside in the range , where I denotes the voxel's intensity and a bit is its unit). Therefore, any of the evaluated neighbor voxels that have the intensity in the defined range are marked as part of the target region by being added to the queue. After all the neighbors of the current voxel are evaluated, the algorithm evaluates the next voxel in the queue; the next voxel becomes the current voxel and its neighbors are evaluated. The procedure is repeated until there are no voxels left in the queue. Finally, all the voxels that are stored in the queue represent the target region which in our case is the liver. Figure 4 displays progress of the segmentation based on region growing after 40.000 evaluated voxels (A), after 400.000 evaluated voxels (B), and after all the voxels in the queue have been evaluated (C).

[Figure omitted. See PDF.]

Figure 4. Progress of the region growing algorithm performing liver segmentation.

Presented is the functioning of the algorithm after 40.000 evaluated voxels (A), after 400.000 evaluated voxels (B), and after all the voxels have been evaluated (C). The initial seed is located at X = 192, Y = 209, Z = 47.

https://doi.org/10.1371/journal.pone.0069068.g004

Due to leakage, region growing may include unwanted segments (e.g. the lower part of the heart, as seen in the upper right part of Figure 4C) which are later eliminated by the postprocessing procedure.

2.1.5. Adaptive Threshold Segmentation Algorithm.

The second liver segmentation algorithm that we evaluated is adaptive threshold algorithm. We developed this algorithm as an upgraded threshold-based algorithm that executes filtering of the current slice using a threshold function while sweeping the intensity threshold value and at the same time comparing the currently filtered slice to the previous already-segmented slice (i.e. the maximum similarity criterion). The intensity threshold value is swept over the whole intensity range, and the similarity comparison is done using normalized cross-correlation which performs segment area comparison (i.e. surface overlap). Similarity to the previous slice is chosen as the criterion for segment determination of adaptive threshold algorithm because the difference in liver shape and size between two neighboring slices is expected to be minimal; therefore, choosing similarity with the previous properly segmented slice as the criterion ensures that the current slice will also be segmented properly. Hence, the error of such a procedure is cumulative and shall the segmentation fail on one slice, all the following slices will be improperly segmented as well. The procedure is shown in Figure 5, where the current slice with three different threshold values (Figures 5A, 5B, and 5C) is compared to the previous slice (Figure 5D).

[Figure omitted. See PDF.]

Figure 5. Adaptive threshold algorithm functioning.

Demonstrated is the functioning of the algorithm comparing a slice while sweeping the intensity threshold value; presented are three examples where the intensity threshold value is set to 15.000 bits and similarity is 0.12 (A), 25.000 bits and similarity is 0.81 (B) and 35.000 bits and similarity is 0.95 (C); the comparison is done to the previous already-segmented slice (D).

https://doi.org/10.1371/journal.pone.0069068.g005

The comparison results (i.e. similarity) which are obtained using normalized cross-correlation are stored during sweep for each intensity threshold value. After the intensity threshold value sweep is done, the intensity threshold value with the highest similarity (i.e. the maximum normalized cross-correlation) is selected, and the current slice is, finally, transformed using the selected intensity threshold value. The procedure is started from the referential slice and repeated on the following slices until the end of the slices collection; then, the procedure is restarted from the referential slice to the beginning of the slices collection. If the highest similarity is lower than a certain similarity threshold, the algorithms determines the segment has ended and empties the current and all the following (or previous, depending on the segmentation Z-direction) slices; in our case, the similarity threshold was set to 0.70 based on our experience on real case data. Moreover, since the referential slice has no prior slices it could be compared to, a set of six presets that include various possible liver segments is used instead, and the maximum similarity to any of the presets indicates the final intensity threshold value for the referential slice. Three out of six example presets that are used for thresholding the referential slice are shown in Figure 6.

[Figure omitted. See PDF.]

Figure 6. Three example liver presets used by the adaptive threshold algorithm.

The preset are used for determining the final intensity threshold value of the referential slice.

https://doi.org/10.1371/journal.pone.0069068.g006

2.1.6. Active Contours Segmentation Algorithm.

The third algorithm for performing liver segmentation that we evaluated for electroporation-based treatment planning is the active contours segmentation, sometimes referred to as the snakes segmentation algorithm. This algorithm is based on placing a deformable (i.e. active) contour, which is a closed curve made of points, on the same location as certain voxels (i.e. the initial contour position). Then, for each point of the active contour (located at the current voxel), the energy of the current and all its neighboring voxels is calculated based on four energy contributions: elasticity of the contour's point in the current voxel, curvature of the contour's point in the current voxel, magnitude of intensity-based energy in the neighboring voxels, and direction of the intensity-based energy in the current voxel. Each point of the active contour is, then, moved to the voxel with the lowest energy. The procedure is repeated until the active contour reaches the desired location (e.g. after a defined number of iterations).

The active contours segmentation algorithm is initiated by placing an initial active contour on top of the referential slice. The active contour is attracted by the edges in the image [26]; therefore, it is imperative to initialize the active contour by placing it near the edge of the desired segment (i.e. the liver). Hence, the adaptive threshold segmentation algorithm is used to generate the initial active contour by segmenting the referential slice and transforming the edge of the segment identified on the referential slice into a closed curve with points sorted according to their location on the segment's circumference. Moreover, the number of all the points in the active contour is reduced by decimation; in our case, the number of circumference pixels between two active contour points is limited to 8. Next, the image energy is calculated as the Gradient Vector Flow (GVF) of the image; a publicly-available GVF calculation algorithm [27] has been implemented using parameter µ = 0.2 and run in 1.000 iterations. The calculation of the GVF is based on the edge map deriving from intensities in a slice; individual steps of this procedure are shown in Figure 7.

[Figure omitted. See PDF.]

Figure 7. Active contours algorithm and the Gradient Vector Flow map.

Presented is an example original liver slice (A) with its edge map (B) and overlaid with a calculated Gradient Vector Flow (GVF) map (C).

https://doi.org/10.1371/journal.pone.0069068.g007

The energy contributions to the total voxel energy are, based on our experience with optimizing algorithms using real case data, balanced using the coefficients 1, 3, 9, 3 for curve elasticity, curve curvature, GVF magnitude, GVF direction, respectively. All the energy contributions are normalized to reside within the range (0,1) in order for the energy coefficients to be properly balanced. In our segmentation algorithm, the curve elasticity energy contribution is calculated based on Equation 3, while the curve curvature energy contribution is calculated based on Equation 4. The implemented GVF algorithm already ensures normalized energy results; the GVF magnitude energy contribution is calculated as the magnitude of the GVF vector in the evaluated voxel, while the GVF direction energy contribution is defined as representing low energy (0) in the neighboring voxel that is located in the direction the GVF vector of the evaluated voxel is pointing to, while all the other voxels have a high energy (1). In Equations 3 and 4, E denotes the energy contribution, while PCURR, PPREV and PNEXT denote the current, the previous and the next points of the active contour, respectively.(3)(4)

After the total energy on and around each active contour point is calculated, and balanced using energy contribution coefficients, the active contour points iteratively move toward the voxel with the lowest total energy. Since the energy depends on the active contour points' locations, the curve elasticity and curve curvature contributions are recalculated after every iteration, and the total energy is recalculated as well. The active contour movement is stopped after a predefined number of iterations (e.g. in our case, the active contour movement is limited to 100 iterations).

When the active contour movement is stopped on one slice, the procedure is repeated on another slice. The procedure is started from the referential slice and is repeated on the next slices until the end of the slices collection; then, the procedure is restarted from the referential slice to the beginning of the slices collection; therefore, the active contour segmentation algorithm may be split into two processing threads.

In order to perform segmentation using the active contours algorithm in three dimensions, the initial contour on the current slice is the same as the final active contour on the previous slice. Because the difference in liver shape and size between two neighboring slices is expected to be minimal on the slices where the organ that is subject to segmentation is present, the previous slice's final active contour is a good initial contour for the current slice. Moreover, due to the expected minimal active contour movement, less iterations are required when calculating the image's GVF since edges will attract active contours in their vicinity even when the GVF is calculated in less iterations; therefore, the processing time of the active contour segmentation algorithm is reduced (e.g. in our case, we are calculating the GVF in 1000 iterations).

2.1.7. Postprocessing.

After segmentation, a postprocessing algorithm needs to be executed to eliminate possible anomalies that may occur during segmentation (e.g. segment leakage). Postprocessing eliminates redundant segments that cannot be part of the final results; elimination is based on comparing neighboring slices in the direction of the third dimension (i.e. component Z) using normalized cross-correlation.

The postprocessing algorithm is initiated on the referential slice; namely, the probability that a slice includes only one identified segment is the highest on the referential slice, since segmentation was initiated on this slice regardless of the chosen segmentation algorithm: for region growing, the seed was placed on this slice; for adaptive threshold, the comparison with presets was made on this slice and also, the active contour was initiated using the adaptive threshold segmentation algorithm on the referential slice as well. Shall the morphological operations during segmentation split the segment on the referential slice into multiple segments, the first step of postprocessing eliminates such redundant segments by only keeping the largest segment on the referential slice.

Then, the postprocessed referential slice, i.e. ZCURR = ZREF, is used as the basis for performing normalized cross-correlation with the next, i.e. ZCURR+1 or the previous, i.e. ZCURR−1 slice, respectively. The template for the normalized cross-correlation is generated by intersecting each segment on the current slice, i.e. ZCURR, with the finally postprocessed previous, i.e. ZCURR-1 or next, i.e. ZCURR+1 slice, respectively (depending on the postprocessing Z-direction). Then, each segment on the current slice is compared to its corresponding template generated from its neighboring slice using normalized cross-correlation; if the result of the comparison exceeds a predefined threshold, the segment is kept on the slice, else it is discarded. In our case, we set the comparison threshold to 0.65 which was found empirically on real case data. The postprocessing procedure may be split into two processing threads, since it is symmetrically executed from Z = ZREF+1 to Z = ZMAX, and from Z = ZREF-1 to Z = 1 (i.e. the postprocessing Z-direction).

Moreover, if we compared the neighboring slices only by using the whole next or previous slice as the template for comparison to the current slice (i.e. without comparing separate segments on a slice), we would be unable to extract these separate segments and determine whether they derive from the target tissue (e.g. liver) or not; therefore, such comparison enables us to keep multiple segments on a slice with the possibility of eliminating segments that are not part of the target tissue. Besides, in order to allow separation of single segments that are in fact multiple segments connected by a thin bond (possibly due to leakage), all the slices are eroded before and dilated after the postprocessing procedure using a disk structuring element of 3×3 pixel size.

2.2. Validation

2.2.1. Optimization using Radiologist Data Set as a Training Set.

Since segmentation algorithms are required to produce not only meaningful but also accurate results, validation of the algorithms is an imperative. In our case, validation was performed as a two-step procedure. In the first step, the algorithms were optimized using radiologist dataset as a training set, and then in the second step, the algorithms were validated after being optimized using another radiologist dataset.

In order to perform segmentation algorithms' optimization, seven sets of patient's liver manually segmented by a radiologist were used as a training set. The patient images that were used for development and optimization of segmentation algorithms belong to patients from the clinical study “Treatment of Liver Metastases with Electrochemotherapy (ECTJ)” (EudraCt no. 2008-008290-54, registered at Clinicaltrials.gov no. NCT01264952). The study was prospective, phase I/II, conducted at the Institute of Oncology Ljubljana, Ljubljana, Slovenia. Regulatory approvals from the Institutional board, as well as from the National Medical Ethics Committee were obtained. Written consents of the patients were obtained. Additional 4 patients were included for final validation of segmentation algorithms: these patients were under standard treatment, according to Declaration of Helsinki, and all patient data were anonymized before processing. For each segmentation algorithm, changeable parameters that significantly influence the functioning of the algorithms were defined, and their possible value ranges were defined based on our previous experience using real-case data. These changeable parameters were then subject to variation within optimization iterations; every of the seven cases that were already manually segmented by the radiologist was re-segmented using our three segmentation algorithms (i.e. region growing, adaptive threshold, and active contours) in the optimization process. During this optimization process, variation of the defined changeable parameters was performed in order to each time automatically obtain a liver object that is most similar to the one segmented by the radiologist. Figure 8 demonstrates how manual segmentation was performed by the radiologist: an example slice during manual segmentation can be seen (Figure 8A), and a final three-dimensional liver object as a result of the manual segmentation by the radiologist (Figure 8B).

[Figure omitted. See PDF.]

Figure 8. Radiologist manual segmentation procedure.

Presented is a defined liver segment on a slice (A) and a final three-dimensional liver object (B).

https://doi.org/10.1371/journal.pone.0069068.g008

After acquiring the data obtained from the radiologist and arranging them into the form that was applicable to optimization and validation (i.e. changing the data syntax by converting them to a raw format, so that inclusion into optimization algorithms was seamless), value ranges of the changeable parameters were defined. For region growing, the parameters subject to variation during optimization were the size of the noise-elimination filter mask during preprocessing, and the threshold deviation value which determines the range of the intensities that are acceptable for inclusion into the target region during segmentation (for optimization purposes, the initial seed of the region growing algorithm was chosen manually and fixed for each segmented case separately). For adaptive threshold, the parameters subject to variation during optimization were the size of the noise-elimination filter mask during preprocessing, and the initial coefficient that determines the targeted size of the referential segment (i.e. the initial liver segment on the referential slice) in the beginning of segmentation. For active contours, all four energy-contribution coefficients were varied during optimization.

Optimization was run on a workstation with Intel Core-i7 965 Extreme Edition processor (maximum frequency 3.46 GHz, 4 cores, 8 threads) with 12 GB of DDR3 memory (frequency 1.333 MHz), two 300 GB hard disk drives (velocity 10.000 rpm) running in stripe mode, and operating system Windows 7 Enterprise (64-bit). The segmented three-dimensional liver object that was produced in each iteration of optimization was compared to the corresponding model that was manually segmented by the radiologist; comparison was done on a slice-by-slice basis using normalized cross-correlation. Similarity (i.e. the results of the normalized cross-correlation) was stored together with the current values of the changeable parameters, and finally the iteration with the maximum similarity (i.e. the mean value of similarity of all slices), was chosen as the optimum (i.e. case-specific similarity, SC in Table 1 and 2); the changeable parameters used in that iteration were marked as optimal parameters for the current case and currently evaluated segmentation algorithm. Moreover, optimal parameters of every evaluated segmentation algorithm were also determined for the whole evaluated series (i.e. all the seven cases) by comparing similarity of different cases with the same changeable parameters to the corresponding model that was manually segmented by the radiologist. These optimized parameters can be defined as globally-optimized parameters since their similarity (SG in Tables 1 and 2) to the training set was evaluated globally (i.e. for all the cases, and not using separate parameters for each case) and is, therefore, estimated that these globally-optimized parameters are optimal for every possible liver that needs to be segmented using our algorithms.

[Figure omitted. See PDF.]

Table 1. Optimization results of seven cases compared to radiologist data by mean values.

https://doi.org/10.1371/journal.pone.0069068.t001

[Figure omitted. See PDF.]

Table 2. Optimization results of seven cases compared to radiologist data by median values.

https://doi.org/10.1371/journal.pone.0069068.t002

2.2.2. Final validation.

After globally optimizing changeable parameters of all three evaluated segmentation algorithms, final validation was performed on additional four cases that were manually segmented by the radiologist. As opposed to the optimization procedure, the changeable parameters were fixed during validation using values obtained from optimization and previously defined as globally-optimized parameters. Then, segmentation was performed and results were compared (i.e. validation similarity, SMN and SMD in Table 3) to the reference model that was manually segmented by the radiologist in the same manner as when performing optimization.

[Figure omitted. See PDF.]

Table 3. Validation results of four cases compared to radiologist data by mean and median values.

https://doi.org/10.1371/journal.pone.0069068.t003

Results

The results of optimization, displayed as similarity to the training set, are shown in Tables 1 and 2. Both the similarity using case-specific optimal parameters (SC in Tables 1 and 2) and using globally-optimized parameters (SG in Tables 1 and 2) are presented. The similarity was evaluated using two possible criteria: as the mean (Table 1) or the median similarity (Table 2) of all the slices within a case of the training set.

As presented in Table 1, after being globally optimized (i.e. for all cases) our implementation of region growing algorithm provides mean slice similarities (SG in Table 1) from 61.1% to 94.2% with the mean value of 79.8% (standard deviation 12.7%) which classifies the region growing as the most accurate algorithm evaluated based on the mean and also the median values of all the slices' similarities. Based on the data from Table 2, median slice similarity values for the region growing algorithm vary from 84.8% to 94.4% with the median value of 89.2%, which is the highest among all three evaluated algorithms.

Our implementation of adaptive threshold algorithm provides globally optimized (i.e. for all cases) mean slice similarities (SG in Table 1) from 40.1% to 80.2% with the mean value of 67.6% (standard deviation 12.8%) which classifies the adaptive threshold algorithm as the least accurate algorithm evaluated based on the mean values of all the slices' similarities. Based on the data from Table 2, median slice similarity values for the adaptive threshold algorithm vary from 55.4% to 92.0% with the median value of 87.6%.

Our implementation of active contours algorithm provides globally optimized (i.e. for all cases) mean slice similarities (SG in Table 1) from 54.8% to 80.3% with the mean value of 68.8% (standard deviation 8.6%). Based on the data from Table 2, median slice similarity values for the active contours algorithm vary from 75.6% to 91.5% with the median value of 86.4%.

The results of validation of the four additional cases manually segmented by the radiologist are shown in Table 3. The similarities are non-optimized: the segmentation results (SV in Table 3) are validated using globally-optimized parameters (SG in Tables 1 and 2) without further modifications of these parameters in order to show functioning of our segmentation algorithms on models which previously were not used as a training set during the optimization procedure.

As presented in Table 3, the region growing algorithm achieves highest similarities (mean 71.9% with standard deviation 16.2%, or median 86.3%) of the validated models. The standard deviation between different models segmented using the same algorithm is the smallest (i.e. 1.4%) using the adaptive threshold algorithm. Individual validation similarities (SV in Table 3) on a slice basis are also presented in Figure 9.

[Figure omitted. See PDF.]

Figure 9. Final validation similarities of four cases (A, B, C, and D).

The similarities presented are on a slice-by-slice basis, i.e. the abscissa represents the location of a slice, while the value of the graphs is the similarity of the slice with the manually validated slice.

https://doi.org/10.1371/journal.pone.0069068.g009

Figure 9 shows that although the four validation cases (A, B, C, and D) were not included in the training set, they could be segmented using our implementations of the region growing, adaptive threshold and active contours algorithms.

Discussion

The aim of this paper was to present three possible automatic liver segmentation algorithms and perform their evaluation by optimizing their functioning for an accurate liver segmentation; the optimization was based on seven-patients dataset segmented by the radiologist and used as a training set, while validation was performed using additional radiologist dataset consisting of four additional patients. As presented in the Results section, developing an automatic segmentation algorithm that generates accurate three-dimensional liver models is demanding, especially because of the variability of input data (i.e. liver in the patient's medical images). Nevertheless, because the region growing requires the end-user to place an initial seed (i.e. to click on the liver on one image), this segmentation algorithm cannot be classified as a fully-automatic algorithm but rather semi-automatic; however, an algorithm for automatically placing the initial seed could upgrade the region growing algorithm to a fully-automatic algorithm. The initial seed-placing algorithm can be developed based on the adaptive threshold segmentation algorithm (described in 2.1.5.) on the referential slice, with an additional task of selecting a pixel with a median intensity of the segment on the referential slice which should ensure avoidance to selecting vessels or nodules as the target region.

The liver includes or may include various inner structures (e.g. blood vessels, metastases, hemangioma, etc.) that do not have a predictable intensity range or texture; these structures may directly interfere with segmentation algorithms by possibly influencing the generation of segments. Also, certain organs near the liver (i.e. the spleen, the heart, the kidneys) have similar intensity ranges to the liver [32], which results in segmentation leakage for threshold-based segmentation methods (i.e. region growing and adaptive threshold). Moreover, active contours are attracted to edges, including the edges of the inner structures; therefore, straight-forward three-dimensional segmentation of the liver using active contours would only be possible if the three-dimensional contour (i.e. the deformable surface) was initialized entirely near the edge of the liver. Besides, the soft transitions of the organs from one slice to another often prevent detection of the organ edges, which leads to an insufficient Gradient Vector Flow (GVF) map and consequently corrupting the generated segments of the active contours algorithm. All the identified difficulties demonstrate that accurate automatic liver segmentation is indeed a challenging task. Fortunately, highest accuracy of liver segmentation is not required for electroporation-based treatment planning since tissue of interest for electric field distribution actually represents the tumor tissue, and possibly large blood vessels which must be taken into account because of electric field calculations and electrode placement. The liver tissue, however, is only a medium that surrounds the tissue of interest, i.e. the tumor that resides in the liver, since electroporation-based treatment planning is required for tumor nodules that are inside or on the edge of the liver. Therefore, the main function of the generated three-dimensional liver model besides being the medium surrounding the tissue of interest is to provide instructions to the attending physician on how the electroporation-based treatment will be performed (i.e. where the tumor is located, since it cannot be seen as it is deep-seated), and serves as electrode-insertion approximate offline navigation for the clinician.

Since similarity comparisons using normalized cross-correlation were performed on each slice of a case separately, final results of each case are presented using two representations: the mean similarity of all the slices, and the median similarity of all the slices in a case. Two representations of the data were presented since there are many individual slices with a similarity of 0%, which means whether the radiologist marked a segment on that slice and the segmentation algorithm missegmented it, or vice-versa (i.e. the segmentation algorithm detected a segment that the radiologist did not identify). Hence, such slices significantly contribute to the final results of the comparison regardless of the size of the segment that caused the 0% similarity on the slice (i.e. even a segment of only few pixels detected by the segmentation algorithms and not identified by the radiologist would produce a 0% similarity of its slice). Since such slices significantly impact the quality of the final results, the mean similarity of all the slices in a case does not reflect success of the segmentation algorithm enough. Therefore, the median similarity of the slices in a case was presented as well.

Optimization results presented in Tables 1 and 2 show that using case-specific optimization parameters (i.e. parameters optimized for every case separately, SC in Tables 1 and 2), very good similarities to the training set can be obtained (median similarities of up to 95.9% using region growing, up to 80.2% using adaptive threshold, and up to 80.3% using active contours). However, results applying to case-specific parameters can only be achieved if a reference model which is part of the training set is available. Therefore, globally-optimized parameters are the parameters that give meaningful information on how accurate a segmentation algorithm is. Our implementation of region growing algorithm could be optimized to achieve an 89.2% median similarity (79.8% mean similarity with 12.7% standard deviation) to the training set, which classifies this algorithm as the most optimization-prone algorithm evaluated. The adaptive threshold algorithm could be optimized to achieve an 87.6% median similarity (67.6% mean similarity with 12.8% standard deviation), and our implementation of the active contours algorithm an 86.4% median similarity (68.8% mean similarity with 8.6% standard deviation), which shows that despite being more sophisticated, these two algorithms achieved lower results than region growing. The main reason is that region growing algorithm can be effectively optimized since the optimizable threshold-deviation parameter majorly influences segmentation (i.e. it significantly influences which voxels will be part of the final tissue segment). The adaptive threshold only has one most relevant parameter that could be optimized (i.e. the initial coefficient that determines the targeted size of the initial segment on the referential slice), which does not influence segmentation on further slices. Although the active contours algorithm includes four optimizable parameters (i.e. the four coefficients that balance energy contributions), it is almost impossible to influence the movement of the active contour (i.e. the snake). Namely, even if the four energy contributions are ideally balanced, the active contour movement needs to mostly rely on the edges in the image (i.e. the two coefficients representing GVF magnitude and directions). Therefore, if the edges of the target tissue were improperly detected after preprocessing (e.g. the target tissue does not have clearly detectable edges), the movement of the active contour will produce unwanted results.

Results of the final validation that are presented in Table 3 show that new patient images that were not part of the optimization training set can be segmented using our implementations of segmentation algorithms (region growing, adaptive threshold, and active contours). The best results are provided by the region growing algorithm (median similarity 86.3%, mean similarity 71.9% with standard deviation 16.2%), while adaptive threshold (median similarity 78.5%, mean similarity 66.3% with standard deviation 1.4%) and active contours (median similarity 80.7%, mean similarity 59.3% with standard deviation 9.4%) provide lower similarity values. As it can be seen in Figure 9, there were many slices missegmented in all four cases. Most missegmented slices were produced using the active contours algorithm, while the least standard deviation is provided using the adaptive threshold algorithm. Although such missegmented slices would negatively influence not only generation of the liver model, but also detection of the tumors, this drawback is avoided using manual validation by the attending physician, which is discussed in the following paragraph. Finally, although the region growing was expected to be the least sophisticated and accurate algorithm, it proved to be the most robust of all the evaluated algorithms, providing highest accuracy and least missegmented slices.

Tumor detection from the patient's medical images is currently still in development and will be implemented as a semi-automatic procedure: identification of structures within the segmented liver (including its edge area) will be done automatically, but finally the end-user (i.e. the attending physician) will be required to manually determine which identified structures are tumors and are, therefore, subject to electroporation-based treatment. Therefore, the whole liver needs to be properly segmented, since tumors are detected as structures within or on the edge of the liver, and could be missed if some slices at the top or at the bottom of the liver are not segmented. Figure 10 demonstrates possible situations that are related to missegmenting liver segments at the beginning or at the end of the liver: Figure 10A shows an ideal case where all the slices are segmented with a high similarity to the radiologist data; in Figure 10B, there are three slices that do not include liver segments but should include them (similarity of these slices is zero); Figure 10C shows a bad case where many slices that should include liver segments do not include them. The main reason for missing segments on such slices is intensity inhomogeneity in the Z direction (i.e. not on a single slice but through slices); although the latter is actually corrected when preprocessing the images applying the sigmoid transformation using Volume of Interest (VOI) parameters from the DICOM header, the anomaly is not completely removed, especially in the beginning and in the end of the series. Since this phenomenon cannot therefore be fully avoided, and because there are rare cases that can be automatically segmented as the case from Fig 10A, manual validation of the segmented slices is required at the end of electroporation-based treatment planning procedure. The validation can be combined with the attending physician's manual identification of the structures whether they are tumors or healthy tissue, which reduces the time needed to execute the whole treatment planning procedure. Hence, the attending physician manually validates the segmented images to ensure proper liver and tumor model generation, and manually corrects them if required.

[Figure omitted. See PDF.]

Figure 10. Possible miss situations when segmenting the liver.

Presented are the situations where all the slices are segmented (A), three slices have been missegmented (B) or many slices have been missegmented (C).

https://doi.org/10.1371/journal.pone.0069068.g010

Currently, only Magnetic Resonance Imaging (MRI) images are segmented for electroporation-based treatment planning, as MRI is preferred in clinical practice for the liver. Namely, MRI is a non-invasive procedure (i.e. there is no ionizing radiation present during MRI) and using a contrast medium it can provide satisfactory images of the liver for electroporation-based treatment planning. Hence, our algorithms are currently written and optimized only for MRI images; the preprocessing procedures are prepared for expected intensity ranges and Volume-of-Interest (VOI) values that derive from MRI sources, and could therefore not perform segmentation of e.g. Computed Tomography (CT) images. However, since the modality of obtained images can be easily determined from the header of the DICOM files, and because only minor modifications of the preprocessing algorithms would be required, a modification allowing non-MRI image segmentation would not be demanding to implement.

Further steps of our research in the field of electroporation-based treatment planning include improved algorithms for tumor detection and, also, automatic vessel segmentation. Namely, segmentation of tumors and vessels needs to be also developed since algorithms for liver segmentation (i.e. region growing, adaptive threshold, and active contours) cannot be used for this purpose without thorough modifications of the algorithms. We expect to develop vessel segmentation algorithms as a combination of intensity- and morphology-based methods with the aim of extracting line-like structures from the liver.

We plan to develop electroporation-based treatment planning as a web application that will be remotely accessible from a web browser and will allow generating treatment plans by allowing the attending physician to upload the DICOM images of the patient and, after calculations and manual validation of the results, obtaining a directly applicable treatment plan. Hence, besides the automatic liver segmentation, the treatment planning software will also need to include electrode insertion and calculation of the electric field distribution with electrode position optimization in order for the treatment plan to comprise all the required information. The treatment planning procedure will be simplified in order to minimize the input of the clinician: because the segmentation of the tissue is automatic, the clinician will only need to validate the segmentation results (and if required, correct the segmentation by dragging the produced contours towards desired positions on each slice) and, finally, determine the entry direction of the electrodes. Namely, the electrodes will be automatically inserted towards the gravitational center of the tumor (or tumors), and the required number of electrodes will be proposed by the software based on the shape and size of the tumor (or tumors). Therefore, the clinician will only need to rotate the electrode array towards the excepted intraoperative entry direction, which will simplify the procedure.

Finally, since automatic liver segmentation can be implemented for other applications [32]–[34] beside electroporation-based treatments, we opt towards extending the functioning of our web-based treatment planning software for electroporation-based treatments onto related fields of surgical liver intervention planning.

Acknowledgments

The first author would like to thank prof. dr. Bostjan Likar for extensive discussion on image segmentation and dr. Barbara Mali for the brief discussion on statistical methods.

Author Contributions

Conceived and designed the experiments: DP GS DM. Performed the experiments: DP MMM. Analyzed the data: DP MMM. Contributed reagents/materials/analysis tools: DP MMM GS. Wrote the paper: DP.

Citation: Pavliha D, Mušič MM, Serša G, Miklavčič D (2013) Electroporation-Based Treatment Planning for Deep-Seated Tumors Based on Automatic Liver Segmentation of MRI Images. PLoS ONE 8(8): e69068. https://doi.org/10.1371/journal.pone.0069068

References

1. Kotnik T, Pucihar G, Miklavcic D (2010) Induced Transmembrane Voltage and Its Correlation with Electroporation-Mediated Molecular Transport. Journal of Membrane Biology 236: 3–13.

2. Miklavcic D (2012) Network for Development of Electroporation-Based Technologies and Treatments: COST TD1104. Journal of Membrane Biology 245: 591–598 .

3. Sersa G, Miklavcic D (2008) Electrochemotherapy of Tumors. Journal of Visualized Experiments 22.

4. Mali B, Jarm T, Snoj M, Sersa G, Miklavcic D (2013) Antitumor effectiveness of electrochemotherapy: A systematic review and meta-analysis. European Journal of Surgical Oncology (EJSO) 39: 4 – 16. doi:10.1016/j.ejso.2012.08.016.

5. Sersa G, Miklavcic D, Cemazar M, Rudolf Z, Pucihar G, et al. (2008) Electrochemotherapy in Treatment of Tumours. EJSO 34: 232–240.

6. Marty M, Sersa G, Garbay J, Gehl J, Collins C, et al. (2006) Electrochemotherapy – an Easy, Highly Effective and Safe Treatment of Cutaneous and Subcutaneous Metastases: Results of ESOPE (European Standard Operating Procedures of Electrochemotherapy) Study. EJC Supplements 4: 3–13.

7. Mir L, Gehl J, Sersa G, Collins C, Garbay J, et al. (2006) Standard operating procedures of the electrochemotherapy: Instructions for the use of bleomycin or cisplatin administered either systemically or locally and electric pulses delivered by the Cliniporator (TM) by means of invasive or non-invasive electrodes. EJC Supplements 4: 14–25 .

8. Edhemovic I, Gadzijev EM, Brecelj E, Miklavcic D, Kos B, et al. (2011) Electrochemotherapy: A New Technological Approach in Treatment of Metastases in the Liver. Technol Cancer Res Treat 10: 475–485.

9. Fini M, Tschon M, Alberghini M, Bianchi G, Mercuri M, et al.. (2011) Cell Electroporation in Bone Tissue. In: Kee ST, Gehl J, Lee EW, editors. Clinical Aspects of Electroporation. New York, NY: Springer New York.

10. Agerholm-Larsen B, Iversen HK, Ibsen P, Moller JM, Mahmood F, et al. (2011) Preclinical Validation of Electrochemotherapy as an Effective Treatment for Brain Tumors. Cancer Research 71: 3753–3762 .

11. Mahmood F, Gehl J (2011) Optimizing clinical performance and geometrical robustness of a new electrode device for intracranial tumor electroporation. Bioelectrochemistry 81: 10–16 .

12. Linnert M, Iversen H, Gehl J (2012) Multiple brain metastases – current management and perspectives for treatment with electrochemotherapy. Radiology and Oncology 2012: 1–8 .

13. Neal II RE, Rossmeisl Jr JH, Garcia PA, Lanz OI, Henao-Guerrero N, et al. (2011) Successful Treatment of a Large Soft Tissue Sarcoma With Irreversible Electroporation. Journal of Clinical Oncology 29: E372–E377 .

14. Miklavcic D, Sersa G, Brecelj E, Gehl J, Soden D, et al. (2012) Electrochemotherapy: technological advancements for efficient electroporation-based treatment of internal tumors. Medical & Biological Engineering & Computing 50: 1213–1225 .

15. Pavliha D, Kos B, Zupanic A, Marcan M, Sersa G, et al. (2012) Patient-specific treatment planning of electrochemotherapy: Procedure design and possible pitfalls. Bioelectrochemistry 87: 265–273.

16. Davalos R, Mir LM, Rubinsky B (2005) Tissue ablation with irreversible electroporation. Annals of Biomedical Engineering 33: 223–231.

17. Zupanic A, Miklavcic D (2009) Optimization and Numerical Modeling in Irreversible Electroporation Treatment Planning. Irreversible Electroporation. Berlin: Springer. 203–222.

18. Miklavcic D, Snoj M, Zupanic A, Kos B, Cemazar M, et al.. (2010) Towards Treatment Planning and Treatment of Deep-Seated Solid Tumors by Electrochemotherapy. Biomedical Engineering Online 9: -.

19. Linguraru MG, Pura JA, Pamulapati V, Summers RM (2012) Statistical 4D graphs for multi-organ abdominal segmentation from multiphase CT. Medical Image Analysis 16: 904–914 .

20. Pahr D, Zysset P (2009) From high-resolution CT data to finite element models: development of an integrated modular framework. Comput Methods Biomech Biomed Engin 12: 45–57.

21. Chi Y, Liu J, Venkatesh SK, Huang S, Zhou J, et al.. (2011) Segmentation of Liver Vasculature From Contrast Enhanced CT Images Using Context-Based Voting. IEEE Transactions on Biomedical Engineering 58: 2144 –2153. doi:10.1109/TBME.2010.2093523.

22. Zupanic A, Corovic S, Miklavcic D (2008) Optimization of Electrode Position and Electric Pulse Amplitude in Electrochemotherapy. Radiology and Oncology 42: 93–101.

23. Zupanic A, Kos B, Miklavcic D (2012) Treatment planning of electroporation-based medical interventions: electrochemotherapy, gene electrotransfer and irreversible electroporation. Physics in Medicine and Biology 57: 5425–5440 .

24. Lecchi M, Fossati P, Elisei F, Orecchia R, Lucignani G (2008) Current Concepts on Imaging in Radiotherapy. European Journal of Nuclear Medicine and Molecular Imaging 35: 821–837.

25. Heimann T, Thorn M, Kunert T, Meinzer HP (2004) New methods for leak detection and contour correction in seeded region growing segmentation. International Archives of the Photogrammetry, Remote Sensing, Vol. 35. Istanbul, Turkey, Vol. 2004: 317–322.

26. Casciaro S, Franchini R, Massoptier L, Casciaro E, Conversano F, et al. (2012) Fully Automatic Segmentations of Liver and Hepatic Tumors From 3-D Computed Tomography Abdominal Images: Comparative Evaluation of Two Automatic Methods. IEEE Sensors Journal 12: 464–473 .

27. Kass M, Witkin A, Terzopoulos D (1987) Snakes – active contour models. Int J of Comp Vis 1: 321–331.

28. Xu C, Prince J (1998) Snakes, shapes, and gradient vector flow. IEEE Trans Imag Proc 7: 359–369.

29. National Electrical Manufacturers Association (2009) Digital Imaging and Communications in Medicine (DICOM).

30. National Electrical Manufacturers Association (2009) CP-467 addition to DICOM. Available: ftp://medical.nema.org/medical/Dicom/Final/cp467_ft.pdf.

31. Vovk U, Pernus F, Likar B (2007) A review of methods for correction of intensity inhomogeneity in MRI. IEEE Transactions on Image Processing 26: 405–421 .

32. Zheng Y, Grossman M, Awate SP, Gee JC (2009) Automatic Correction of Intensity Nonuniformity from Sparseness of Gradient Distribution in Medical Images. Proceedings of the 12th International Conference on Medical Image Computing and Computer-Assisted Intervention: Part II. MICCAI '09. Berlin, Heidelberg: Springer-Verlag. 852–859. Available: http://dx.doi.org/10.1007/978-3-642-04271-3_103.

33. Massoptier L, Casciaro S (2008) A new fully automatic and robust algorithm for fast segmentation of liver tissue and tumors from CT scans. European Radiology 18: 1658–1665 .

34. Crocetti L, Lencioni R, DeBeni S, See T, Della Pina C, et al.. (2008) Targeting Liver Lesions for Radipfrequency Ablation – an Experimental Feasibility Study Using a CT-US Fusion Imaging System. Investigative Radiology: 33–39.