Image Segmentation for the Treatment Planning of Magnetic Resonance-Guided High-Intensity Focused Ultrasound (MRgHIFU) Therapy: A Parametric Study

: In the present research work, image segmentation methods were studied to ﬁnd internal parameters that provide an e ﬃ cient identiﬁcation of the regions of interest in Magnetic Resonance (MR) images used for the therapy planning of High-Intensity Focused Ultrasound (HIFU), a minimally invasive therapeutic method used for selective ablation of tissue. The involved image image segmentation for treatment planning of the MRgHIFU therapy. Future work will address the search of an automatic segmentation process, regardless of the experimental setup. and they were applied to 112 transverse and 112 sagittal MR images, obtained from a study for the treatment of abscesses with MRgHIFU in a murine model [3]. In the images, the regions of interest were the air, tissue, gel-pad, transducer and water regions. The evaluation of the segmentation quality was performed using F-measure [4] supported by ground truth images obtained from manual delineation of the regions of interest. A parametric study consisting of tests with different values for the internal parameters of the implemented image segmentation methods was carried out. The F-measure results from the parametric study were analyzed independently by region of interest in order to determine which method is the most appropriate to identify each region.


Introduction
High-Intensity Focused Ultrasound (HIFU) is a minimally invasive therapeutic method in which ultrasound beams are concentrated at a focal region, producing heating and selective ablation within the focal volume without compromising surrounding tissues [1]. HIFU has been proposed for the safe ablation of both malignant and benign tissues and as an agent for drug delivery, and Magnetic Resonance Imaging (MRI) has been proposed for guidance and monitoring of the HIFU therapy [2]. The combination of HIFU and MRI is known as Magnetic Resonance-guided HIFU (MRgHIFU). The HIFU treatment is performed with an ultrasonic transducer, acoustically coupled to the target tissue using water and gel-pads, as shown in Figure 1. In MRgHIFU therapy planning, identifying the regions allows for accurate targeting of tissue, positioning of the transducer and its geometric focus (focal point). If a proper automated identification of the transducer and tissue region can be achieved, with image segmentation, targeting would be facilitated, as well as therapy for planning using numerical methods for the calculation of acoustic energy distribution and thermal response in biological tissue.
In the present research work, segmentation methods were studied to find internal parameters that provide an efficient identification of the regions of interest in MR images used for HIFU treatment planning. The involved segmentation methods were threshold, level set and watershed with markers, and they were applied to 112 transverse and 112 sagittal MR images, obtained from a study for the treatment of abscesses with MRgHIFU in a murine model [3]. In the images, the regions of interest were the air, tissue, gel-pad, transducer and water regions. The evaluation of the segmentation quality was performed using F-measure [4] supported by ground truth images obtained from manual delineation of the regions of interest. A parametric study consisting of tests with different values for the internal parameters of the implemented image segmentation methods was carried out. The F-measure results from the parametric study were analyzed independently by region of interest in order to determine which method is the most appropriate to identify each region.

The Experimental Setup
Experimental data was recovered from a previous study in a murine model [3] where subcutaneous abscesses were induced by methicillin-resistant Staphylococcus aureus and treated with MRgHIFU. Fifty-five female BALB/c mice, aged 7-12 weeks were used. All animal experiments were In MRgHIFU therapy planning, identifying the regions allows for accurate targeting of tissue, positioning of the transducer and its geometric focus (focal point). If a proper automated identification of the transducer and tissue region can be achieved, with image segmentation, targeting would be facilitated, as well as therapy for planning using numerical methods for the calculation of acoustic energy distribution and thermal response in biological tissue.
In the present research work, segmentation methods were studied to find internal parameters that provide an efficient identification of the regions of interest in MR images used for HIFU treatment planning. The involved segmentation methods were threshold, level set and watershed with markers, and they were applied to 112 transverse and 112 sagittal MR images, obtained from a study for the treatment of abscesses with MRgHIFU in a murine model [3]. In the images, the regions of interest were the air, tissue, gel-pad, transducer and water regions. The evaluation of the segmentation quality was performed using F-measure [4] supported by ground truth images obtained from manual delineation of the regions of interest. A parametric study consisting of tests with different values for the internal parameters of the implemented image segmentation methods was carried out. The F-measure results from the parametric study were analyzed independently by region of interest in order to determine which method is the most appropriate to identify each region.

The Experimental Setup
Experimental data was recovered from a previous study in a murine model [3] where subcutaneous abscesses were induced by methicillin-resistant Staphylococcus aureus and treated with MRgHIFU. Fifty-five female BALB/c mice, aged 7-12 weeks were used. All animal experiments were run following an approved protocol according to institutional and Canadian Council of Animal Care guidelines (Lakehead University protocol AUP #14 2011/ROMEO #1461984). A MR compatible single-element transducer (model 10-09-11TBHC, FUS Instruments, Toronto, Canada) with a focal length of 50 mm, diameter of 32 mm and operating at 3 MHz was used in the study. The focal point was positioned 2 mm under the skin at the abscess center. The experimental setup for the study is shown in Figure 2.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 24 run following an approved protocol according to institutional and Canadian Council of Animal Care guidelines (Lakehead University protocol AUP #14 2011/ROMEO #1461984). A MR compatible single-element transducer (model 10-09-11TBHC, FUS Instruments, Toronto, Canada) with a focal length of 50 mm, diameter of 32 mm and operating at 3 MHz was used in the study. The focal point was positioned 2 mm under the skin at the abscess center. The experimental setup for the study is shown in Figure 2.

Magnetic Resonance Images
Transverse and sagittal T1-weighted MR images of the experimental setup were obtained with a 3T MRI scanner (Achieva, Philips Healthcare) using a Flex-S coil (Philips Healthcare). The parameters for the generation of the MR images were a gradient echo sequence (GRE), a field-of-view (FOV) of 120 × 120 × 48 mm, pixel size of 0.5 mm, slice thickness of 2 mm, echo time/repetition time (TE/TR) of 2.5/4.9 msec, flip angle of 35°, acquisition matrix of 120 × 100, reconstruction matrix of 240 × 240 and number of excitations (NEX) of 2 [3]. The regions of interest to be identified are shown in transverse and sagittal images in Figures 3 and 4, respectively.

Magnetic Resonance Images
Transverse and sagittal T1-weighted MR images of the experimental setup were obtained with a 3T MRI scanner (Achieva, Philips Healthcare) using a Flex-S coil (Philips Healthcare). The parameters for the generation of the MR images were a gradient echo sequence (GRE), a field-of-view (FOV) of 120 × 120 × 48 mm, pixel size of 0.5 mm, slice thickness of 2 mm, echo time/repetition time (TE/TR) of 2.5/4.9 msec, flip angle of 35 • , acquisition matrix of 120 × 100, reconstruction matrix of 240 × 240 and number of excitations (NEX) of 2 [3]. The regions of interest to be identified are shown in transverse and sagittal images in Figures 3 and 4, respectively.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 24 run following an approved protocol according to institutional and Canadian Council of Animal Care guidelines (Lakehead University protocol AUP #14 2011/ROMEO #1461984). A MR compatible single-element transducer (model 10-09-11TBHC, FUS Instruments, Toronto, Canada) with a focal length of 50 mm, diameter of 32 mm and operating at 3 MHz was used in the study. The focal point was positioned 2 mm under the skin at the abscess center. The experimental setup for the study is shown in Figure 2.

Magnetic Resonance Images
Transverse and sagittal T1-weighted MR images of the experimental setup were obtained with a 3T MRI scanner (Achieva, Philips Healthcare) using a Flex-S coil (Philips Healthcare). The parameters for the generation of the MR images were a gradient echo sequence (GRE), a field-of-view (FOV) of 120 × 120 × 48 mm, pixel size of 0.5 mm, slice thickness of 2 mm, echo time/repetition time (TE/TR) of 2.5/4.9 msec, flip angle of 35°, acquisition matrix of 120 × 100, reconstruction matrix of 240 × 240 and number of excitations (NEX) of 2 [3]. The regions of interest to be identified are shown in transverse and sagittal images in Figures 3 and 4, respectively.

Ground Truth Images
The regions of interest of 112 transverse and 112 sagittal MR images were manually segmented (delineated) to obtain ground truth images as shown in Figure 5. This information was used as a reference to evaluate the results, by region of interest, from the image segmentation methods described in Section 2.4.

Image Segmentation Methods for the Identification of Regions of Interest
Two image segmentation approaches for the identification of the regions of interest in the MR image were implemented: image segmentation with threshold and level set methods (TLSM) and image segmentation with watershed segmentation algorithm with markers (WSAM).

Image Segmentation with Threshold and Level Set Methods (TLSM)
In this approach, each region of interest was segmented using different image segmentation methods.

Division of the MR Image (Preprocessing)
The original MR image was divided into two parts: the upper part where the air, gel-pad and tissue regions are located and the lower part where the water and the transducer regions are located as shown in Figure 6a,b, respectively. The user is required to perform this division in order to facilitate the application of segmentation methods.

Ground Truth Images
The regions of interest of 112 transverse and 112 sagittal MR images were manually segmented (delineated) to obtain ground truth images as shown in Figure 5. This information was used as a reference to evaluate the results, by region of interest, from the image segmentation methods described in Section 2.4.

Ground Truth Images
The regions of interest of 112 transverse and 112 sagittal MR images were manually segmented (delineated) to obtain ground truth images as shown in Figure 5. This information was used as a reference to evaluate the results, by region of interest, from the image segmentation methods described in Section 2.4.

Image Segmentation Methods for the Identification of Regions of Interest
Two image segmentation approaches for the identification of the regions of interest in the MR image were implemented: image segmentation with threshold and level set methods (TLSM) and image segmentation with watershed segmentation algorithm with markers (WSAM).

Image Segmentation with Threshold and Level Set Methods (TLSM)
In this approach, each region of interest was segmented using different image segmentation methods.

Division of the MR Image (Preprocessing)
The original MR image was divided into two parts: the upper part where the air, gel-pad and tissue regions are located and the lower part where the water and the transducer regions are located as shown in Figure 6a,b, respectively. The user is required to perform this division in order to facilitate the application of segmentation methods.

Image Segmentation Methods for the Identification of Regions of Interest
Two image segmentation approaches for the identification of the regions of interest in the MR image were implemented: image segmentation with threshold and level set methods (TLSM) and image segmentation with watershed segmentation algorithm with markers (WSAM).

Image Segmentation with Threshold and Level Set Methods (TLSM)
In this approach, each region of interest was segmented using different image segmentation methods.

Division of the MR Image (Preprocessing)
The original MR image was divided into two parts: the upper part where the air, gel-pad and tissue regions are located and the lower part where the water and the transducer regions are located as shown in Figure 6a,b, respectively. The user is required to perform this division in order to facilitate the application of segmentation methods. Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 24

Segmentation of Air Region
The air region in the upper part of the MR image, shown in Figure 6a, appeared as an almost homogeneous dark region, so it was easily identifiable using the threshold method [5] as follows: where: is the threshold that separates the image into two dominant modes. The value of this parameter can be varied, so it was considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called T_air.
From the upper part of the MR image, shown in Figure 7a, a binary image was generated with the application of the threshold method as shown in Figure 7b, with an area corresponding to a waterbased heating pad appearing as multiple separate areas around the tissue region. These areas were removed by using the region growing algorithm (RGA). The RGA works on the original input image ( , ), using a similarity criterion and a set of seed points [5]. The RGA used in the current study was a MATLAB ® (2017 MathWorks, Inc., Natick, MA, USA) implementation that worked on the grayscale input MR image ( , ), a similarity criterion and a single seed point [6]. The RGA yielded a binary image ( , ) with the segmented region [5,6]. In the current study the seed point was set inside the largest separate area in the binary image, shown in Figure 7b, and the binary image could be obtained without the areas corresponding to the heating pad as shown in Figure 7c. The RGA had a similarity criterion defined as the difference between the intensity value of the pixel and the mean intensity of the region. The pixels with a difference between their intensity values and the mean intensity below a fixed value of 0.1 on a normalized image were considered to meet this criterion.

Segmentation of Air Region
The air region in the upper part of the MR image, shown in Figure 6a, appeared as an almost homogeneous dark region, so it was easily identifiable using the threshold method [5] as follows: where: T is the threshold that separates the image into two dominant modes. The value of this parameter can be varied, so it was considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called T_air.
From the upper part of the MR image, shown in Figure 7a, a binary image was generated with the application of the threshold method as shown in Figure 7b, with an area corresponding to a water-based heating pad appearing as multiple separate areas around the tissue region. These areas were removed by using the region growing algorithm (RGA). The RGA works on the original input image I(x, y), using a similarity criterion Q and a set of seed points [5]. The RGA used in the current study was a MATLAB ® (2017 MathWorks, Inc., Natick, MA, USA) implementation that worked on the grayscale input MR image I(x, y), a similarity criterion and a single seed point [6]. The RGA yielded a binary image J(x, y) with the segmented region [5,6]. In the current study the seed point was set inside the largest separate area in the binary image, shown in Figure 7b, and the binary image could be obtained without the areas corresponding to the heating pad as shown in Figure 7c. The RGA had a similarity criterion defined as the difference between the intensity value of the pixel and the mean intensity of the region. The pixels with a difference between their intensity values and the mean intensity below a fixed value of 0.1 on a normalized image were considered to meet this criterion.

Segmentation of Air Region
The air region in the upper part of the MR image, shown in Figure 6a, appeared as an almost homogeneous dark region, so it was easily identifiable using the threshold method [5] as follows: where: is the threshold that separates the image into two dominant modes. The value of this parameter can be varied, so it was considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called T_air.
From the upper part of the MR image, shown in Figure 7a, a binary image was generated with the application of the threshold method as shown in Figure 7b, with an area corresponding to a waterbased heating pad appearing as multiple separate areas around the tissue region. These areas were removed by using the region growing algorithm (RGA). The RGA works on the original input image ( , ), using a similarity criterion and a set of seed points [5]. The RGA used in the current study was a MATLAB ® (2017 MathWorks, Inc., Natick, MA, USA) implementation that worked on the grayscale input MR image ( , ), a similarity criterion and a single seed point [6]. The RGA yielded a binary image ( , ) with the segmented region [5,6]. In the current study the seed point was set inside the largest separate area in the binary image, shown in Figure 7b, and the binary image could be obtained without the areas corresponding to the heating pad as shown in Figure 7c. The RGA had a similarity criterion defined as the difference between the intensity value of the pixel and the mean intensity of the region. The pixels with a difference between their intensity values and the mean intensity below a fixed value of 0.1 on a normalized image were considered to meet this criterion.

Segmentation of Gel-Pad Region
The gel-pad region in the upper part of the MR image, shown in Figure 6a, was segmented using two level set methods independently: the Geodesic Active Contours (GAC) method [7] and the Distance Regularized Level Set Evolution (DRLSE) method [8]. These methods are useful for the detection of boundaries. The GAC method [8] is defined as follows: where: u is the level set function. The initial form of this function is the initial contour, denoted as u 0 . Its shape can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called u_0. c is the constant velocity. The value for this parameter can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called c_gelpad. I is the image. g(I) is the stopping function.
The stopping function is given by the following Equation: where: I is the smoothed version of the image I using a Gaussian filter. The filter size and its standard deviation σ are fixed values and they were set to 5 × 5 and 1 respectively.
p is an exponent with a fixed value of 2.
The implementation for the GAC method used in the current study, and the values of its fixed parameters, are based on the implementation reported in [9].
The DRLSE method [8] is defined as follows where: φ is the level set function. The initial form of this contour is the initial contour, denoted as φ 0 . Its shape can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called phi_0. µ is the coefficient of the distance regularization term. For this parameter, a fixed value of 0.2 was considered. λ is the coefficient of the weighted length term. For this parameter, a fixed value of 5 was considered.
g is the edge indicator function. α is the weighted area coefficient. The value for this parameter can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called a_gelpad. δ ε is the Dirac delta approximation function.
d p is a function in terms of potential function.
The edge indicator function g is defined as where: Appl. Sci. 2019, 9, 5296 7 of 23 G σ is a Gaussian kernel with standard deviation σ. The filter size and its standard deviation σ are fixed values and they were set to 5 × 5 and 1, respectively. The values of these parameters were chosen to be the same as those of the stopping function in the implementation of the GAC method.
I is the image.
The Dirac delta approximation function δ ε is defined as where: ε is a parameter with a fixed value of 1.5.
d p is a function in terms of potential function defined as where: s is the magnitude of the gradient of the level set function φ defined as The implementation for the DRLSE method used in the current study and the values of its fixed parameters are based on the implementation reported in [10]. For the gel-pad detection, both the GAC and DRLSE methods require an initial contour that is based on the calculation of the maximum rectangle inside an arbitrary shape defined by a binary image [11]. The algorithm used in the current study is a MATLAB ® (2017 MathWorks, Inc., Natick, MA, USA) implementation that works with a binary input image I(x, y), shown in Figure 7c, the specification of an optimization criteria and the minimum dimensions to restrict the calculation. Finally, the algorithm yields a binary image M(x, y) with the mask of the maximum rectangle inside the arbitrary shape.
The obtained maximum rectangle is shown in Figure 8a. The initial contour, u_0 for the GAC method and phi_0 for the DRLSE method, was a reduced version of the maximum rectangle and its dimensions can be varied, so it will be considered in the parametric study described in Section 2.6. An example of one of the reduced versions is shown in Figure 8b. In both the GAC and DRLSE methods, the initial contour changes its shape iteratively, and the algorithm stops at the detected borders. The final contour has the shape of the found object as shown in Figure 8c. The segmentation of the gel-pad region, together with the previous segmentation of the air region, gives an enclosed region corresponding to the tissue region. The resulting segmented image is shown in Figure 8d. The reduced version of the maximum rectangle consisted of a fixed reduction of its width by five pixels from the right and left sides, and a variable reduction of its height from the lower side by 25 to 32 pixels. Hereinafter, the set of the different reduced versions of the maximum rectangle will be called W − i, where i ∈ Z : i ∈ [25,32] is the number of pixels used for the reduction. Another set of reduced versions was also considered in the parametric study described in Section 2.6. This set has the same characteristics of the previous set and, in addition, features a displacement of five pixels from the upper side. Hereinafter, this set will be called Another approach for the generation of the initial contour was explored. An initial contour based on the edge indicator function of Equation (5) was generated. The process is shown in Figure 9 and was carried out as follows: The edge indicator function was applied to the upper part of the original MR image shown in Figure 9a. The result is shown in Figure 9b. A threshold of 0.2 in a grayscale of 256 levels between 0 and 1 was applied to obtain the image shown in Figure 9c. A mask from the image shown in Figure 9c was obtained using the RGA with a seed point in (3,3). The complement of the resulting image was calculated, and the result is shown in Figure 9d. The image shown in Figure 9c was multiplied by this image to limit the number of areas. The result is shown in Figure 9e where the larger area corresponds to the gel-pad region. The centroid of this large area is shown as a red spot. The gel pad area was then extracted using the RGA with the centroid as the seed point. The extracted region is shown in Figure 9f, and it was used as a mask for the extraction of the corresponding region in the original MR image shown in Figure 9a. The result is shown in Figure 9g. The initial contour was obtained after the application of a threshold of 0 in a grayscale of 256 levels between 0 and 255, in the image shown in Figure 9g. In Figure 9h the initial contour is shown superimposed in the original MR image. This approach for the initial contour (hereinafter this approach will be called the automatic initial contour) was used by both the GAC and DRLSE methods to obtain a final contour with the shape of the found gel-pad region as shown in Figure 9i. The segmentation of the gel-pad region, together with the previous segmentation of the air region, gives an enclosed region corresponding to the tissue region. The resulting segmented image is shown in Figure 9j.  Another approach for the generation of the initial contour was explored. An initial contour based on the edge indicator function of Equation (5) was generated. The process is shown in Figure 9 and was carried out as follows: The edge indicator function was applied to the upper part of the original MR image shown in Figure 9a. The result is shown in Figure 9b. A threshold of 0.2 in a grayscale of 256 levels between 0 and 1 was applied to obtain the image shown in Figure 9c. A mask from the image shown in Figure 9c was obtained using the RGA with a seed point in (3,3). The complement of the resulting image was calculated, and the result is shown in Figure 9d. The image shown in Figure 9c was multiplied by this image to limit the number of areas. The result is shown in Figure 9e where the larger area corresponds to the gel-pad region. The centroid of this large area is shown as a red spot. The gel pad area was then extracted using the RGA with the centroid as the seed point. The extracted region is shown in Figure 9f, and it was used as a mask for the extraction of the corresponding region in the original MR image shown in Figure 9a. The result is shown in Figure 9g. The initial contour was obtained after the application of a threshold of 0 in a grayscale of 256 levels between 0 and 255, in the image shown in Figure 9g. In Figure 9h the initial contour is shown superimposed in the original MR image. This approach for the initial contour (hereinafter this approach will be called the automatic initial contour) was used by both the GAC and DRLSE methods to obtain a final contour with the shape of the found gel-pad region as shown in Figure 9i. The segmentation of the gel-pad region, together with the previous segmentation of the air region, gives an enclosed region corresponding to the tissue region. The resulting segmented image is shown in Figure 9j. Another approach for the generation of the initial contour was explored. An initial contour based on the edge indicator function of Equation (5) was generated. The process is shown in Figure 9 and was carried out as follows: The edge indicator function was applied to the upper part of the original MR image shown in Figure 9a. The result is shown in Figure 9b. A threshold of 0.2 in a grayscale of 256 levels between 0 and 1 was applied to obtain the image shown in Figure 9c. A mask from the image shown in Figure 9c was obtained using the RGA with a seed point in (3,3). The complement of the resulting image was calculated, and the result is shown in Figure 9d. The image shown in Figure 9c was multiplied by this image to limit the number of areas. The result is shown in Figure 9e where the larger area corresponds to the gel-pad region. The centroid of this large area is shown as a red spot. The gel pad area was then extracted using the RGA with the centroid as the seed point. The extracted region is shown in Figure 9f, and it was used as a mask for the extraction of the corresponding region in the original MR image shown in Figure 9a. The result is shown in Figure 9g. The initial contour was obtained after the application of a threshold of 0 in a grayscale of 256 levels between 0 and 255, in the image shown in Figure 9g. In Figure 9h the initial contour is shown superimposed in the original MR image. This approach for the initial contour (hereinafter this approach will be called the automatic initial contour) was used by both the GAC and DRLSE methods to obtain a final contour with the shape of the found gel-pad region as shown in Figure 9i. The segmentation of the gel-pad region, together with the previous segmentation of the air region, gives an enclosed region corresponding to the tissue region. The resulting segmented image is shown in Figure 9j.

Segmentation of Tissue Region
After the segmentation of the air and gel-pad regions, a resulting region is enclosed between them. This region is the segmented tissue region, as shown in Figures 8d and 9j. Therefore, an appropriate segmentation of the tissue region depends on the performance of the segmentation methods for the air and gel-pad regions.

Segmentation of Transducer and Water Regions
The transducer region in the lower part of the MR image, shown in Figure 6b, was segmented in the same way as the gel-pad region, using two level set methods independently: the GAC method [7] and the DRLSE method [8]. In the GAC method, the constant velocity c can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called c_transducer. In the DRLSE method, the weighted area coefficient α can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called a_transducer.
The initial contour φ 0 , for both the GAC and DRLSE methods, was the binary step function suggested in [8] and is defined as where: c 0 is a constant and c 0 > 0. For this constant, a fixed value of 2 was considered. This value was used in the implementation for the DRLSE method reported in [10]. R 0 is a region in the domain of the image.
The obtained initial contour with the binary step function is shown in Figure 10a. The final contour has the shape of the found object as shown in Figure 10b, and the segmented region is shown in Figure 10c. The background of the segmented region is the water region.

Segmentation of Tissue Region
After the segmentation of the air and gel-pad regions, a resulting region is enclosed between them. This region is the segmented tissue region, as shown in Figures 8d and 9j. Therefore, an appropriate segmentation of the tissue region depends on the performance of the segmentation methods for the air and gel-pad regions.

Segmentation of Transducer and Water Regions
The transducer region in the lower part of the MR image, shown in Figure 6b, was segmented in the same way as the gel-pad region, using two level set methods independently: the GAC method [7] and the DRLSE method [8]. In the GAC method, the constant velocity can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called c_transducer. In the DRLSE method, the weighted area coefficient can be varied, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called a_transducer.
The initial contour , for both the GAC and DRLSE methods, was the binary step function suggested in [8] and is defined as where: is a constant and > 0. For this constant, a fixed value of 2 was considered. This value was used in the implementation for the DRLSE method reported in [10]. is a region in the domain of the image. The obtained initial contour with the binary step function is shown in Figure 10a. The final contour has the shape of the found object as shown in Figure 10b, and the segmented region is shown in Figure 10c. The background of the segmented region is the water region. The threshold method was also used for the segmentation of the transducer region in the lower part of the MR image. From Figure 11a, a total of nine thresholds were equidistantly defined between the two maximal peaks in the histogram of the image, with the first and last defined thresholds as the two maximal peaks. A threshold, between the nine defined, can be chosen to perform segmentation as shown in Figure 11b, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called Tn_transducer. Finally, the segmented transducer is shown in Figure 11c. The threshold method was also used for the segmentation of the transducer region in the lower part of the MR image. From Figure 11a, a total of nine thresholds were equidistantly defined between the two maximal peaks in the histogram of the image, with the first and last defined thresholds as the two maximal peaks. A threshold, between the nine defined, can be chosen to perform segmentation as shown in Figure 11b, so it will be considered in the parametric study described in Section 2.6. Hereinafter this parameter will be called Tn_transducer. Finally, the segmented transducer is shown in Figure 11c.

Image Segmentation with Watershed Segmentation Algorithm with Markers (WSAM)
An approach based on the WSAM was developed. The image has sections of regional minima, and the sets of coordinates of the points within them are , ,…, . There is a catchment basin associated with each regional minimum ; then, ( ) is the set of coordinates of the points in the catchment basin. In the algorithm, the topography of the image is flooded from a minimum level = ( , ) + 1 to a maximum level = ( , ) + 1 to find regions and finally establish their borders with watershed lines. Markers are used in the algorithm for an efficient identification of the desired regions, preventing oversegmentation [5]. The WSAM method used in the current study was implemented in OpenCV [12,13]. The original MR image was again divided into an upper and a lower part of the original MR image. In the following paragraphs, the process to generate markers for the WSAM for each region is described. The process is also illustrated in the diagram shown in Figure A1 (Appendix).

Definition of Markers for Gel-Pad and Tissue Regions
The marker for the gel-pad region was based on the maximum rectangle inside an arbitrary shape, defined by a binary image [14]. The binary image was obtained from the upper part of the original MR image, shown in Figure 6a, after the application of the threshold method, shown in Equation (1) (threshold of 13 in a grayscale of 256 levels). The marker of the gel-pad region consisted The final segmented image is built from the individual segmented regions. This image is shown in Figure 12.

Image Segmentation with Watershed Segmentation Algorithm with Markers (WSAM)
An approach based on the WSAM was developed. The image has sections of regional minima, and the sets of coordinates of the points within them are , ,…, . There is a catchment basin associated with each regional minimum ; then, ( ) is the set of coordinates of the points in the catchment basin. In the algorithm, the topography of the image is flooded from a minimum level = ( , ) + 1 to a maximum level = ( , ) + 1 to find regions and finally establish their borders with watershed lines. Markers are used in the algorithm for an efficient identification of the desired regions, preventing oversegmentation [5]. The WSAM method used in the current study was implemented in OpenCV [12,13]. The original MR image was again divided into an upper and a lower part of the original MR image. In the following paragraphs, the process to generate markers for the WSAM for each region is described. The process is also illustrated in the diagram shown in Figure A1 (Appendix).

Definition of Markers for Gel-Pad and Tissue Regions
The marker for the gel-pad region was based on the maximum rectangle inside an arbitrary shape, defined by a binary image [14]. The binary image was obtained from the upper part of the original MR image, shown in Figure 6a, after the application of the threshold method, shown in Equation (1) (threshold of 13 in a grayscale of 256 levels). The marker of the gel-pad region consisted

Image Segmentation with Watershed Segmentation Algorithm with Markers (WSAM)
An approach based on the WSAM was developed. The image has sections of regional minima, and the sets of coordinates of the points within them are M 1 , M 2 , . . . , M R . There is a catchment basin associated with each regional minimum M i ; then, C(M i ) is the set of coordinates of the points in the catchment basin. In the algorithm, the topography of the image is flooded from a minimum level n = min( f (x, y)) + 1 to a maximum level n = max( f (x, y)) + 1 to find regions and finally establish their borders with watershed lines. Markers are used in the algorithm for an efficient identification of the desired regions, preventing oversegmentation [5]. The WSAM method used in the current study was implemented in OpenCV [12,13].
The original MR image was again divided into an upper and a lower part of the original MR image. In the following paragraphs, the process to generate markers for the WSAM for each region is described. The process is also illustrated in the diagram shown in Figure A1 (Appendix A).

Definition of Markers for Gel-Pad and Tissue Regions
The marker for the gel-pad region was based on the maximum rectangle inside an arbitrary shape, defined by a binary image [14]. The binary image was obtained from the upper part of the original MR image, shown in Figure 6a, after the application of the threshold method, shown in Equation (1) (threshold of 13 in a grayscale of 256 levels). The marker of the gel-pad region consisted of a reduced version of the maximum rectangle. Its width was reduced by 10 pixels from the right and left sides and its height by 10 pixels from the upper and lower sides.
The marker for the tissue region was obtained from the upper part of the original MR image, shown in Figure 6a. This image was divided taking the upper side of the maximum rectangle as frontier. A binary image was obtained after the use of the threshold method considering (threshold of 13 in a grayscale of 256). In this image, a median blur filter [15] was used to smooth the tissue contours. The resulting image is used to define the marker of the tissue region using the border following algorithm [16,17] to detect contours and the Ramer-Douglas-Peucker algorithm [18][19][20] to obtain a curve with fewer points from the detected contours. Finally, the obtained curve is drawn as using the 'drawContours' function in OpenCV [21].

Definition of Marker for Air Region
The binary image that was used for the calculation of the marker for the tissue region was revisited to obtain the marker for the air region, using the border following an algorithm [16,17] to detect contours and the Ramer-Douglas-Peucker algorithm [18][19][20] to obtain a curve with fewer points from the detected contours. Finally, the obtained curve was drawn as using the 'drawContours' function in OpenCV [21]. The resulting image was filtered using a convolution filter (filter2D) [22]. The filtering was repeated until there were no black pixels within the tissue. After this process, the threshold method in inverted threshold configuration type [23] was applied. From the resulting image, a rectangular region was extracted for the definition of an image with the air region marker. Finally, this image was eroded [24] two times to get a separation between the markers of the air and tissue regions.

Definition of Markers for Transducer and Water Regions
To obtain the marker of the transducer region, a binary image was obtained from the lower part of the original MR image, shown in Figure 6b, after the application of the threshold method in inverted threshold configuration type. A median blur filter was used to remove pixels surrounding the transducer region in the binary image. In the same way as in the definition of the marker for the tissue region, the border following algorithm [16,17], the Ramer-Douglas-Peucker algorithm [18][19][20] and the 'drawContours' function in OpenCV [21] were used for the detection of contours, for the calculation of a curve with fewer points from the detected contours and for the drawing of the obtained curve, respectively. A rectangular region was extracted for the definition of an image with the marker of the transducer region.
To obtain the marker of the water region, these algorithms were also used followed by the application of the threshold method in inverted threshold configuration type [23]. Finally, the extraction of a rectangular region and an erosion [24] were performed for the definition of an image with the marker for the water region.
To obtain the marker of the hole inside the transducer region, a binary image was obtained from the lower part of the original MR image, shown in Figure 6b, after the application of the threshold method in inverted threshold configuration type. Holes in the image were filled, and a maximum rectangle inside the arbitrary shape in the resulting binary image was calculated [14]. A rectangular region with the hole was extracted from the original MR image using the calculated maximum rectangle, and a binary image was obtained after the application of the threshold method in inverted threshold configuration type. The pixels of the hole region are identified, and the marker was defined after the application of a median blur filter (kernel of 3).

Elements of the Definition of Markers in OpenCV
The border following algorithm was used in the definition of markers for the air, tissue, transducer and water regions. The retrieval mode was set to CV_RETR_CCOMP to retrieve and organize the contours in a hierarchy of two levels, and the contour approximation method was set to CV_CHAIN_APPROX_SIMPLE to obtain the end points of encoded linear segments.
The Ramer-Douglas-Peucker algorithm was used in the definition of markers for the air, tissue, transducer and water regions The parameter of accuracy was 0, and a closed curve was considered for the approximation.
The 'drawContours' function was used in the definition of markers for the air, tissue, transducer and water regions. The thickness was set to a value of 1 in the definition of the markers of the tissue and transducer regions, to a value of −1 for the marker of the air region and to a value of CV_FILLED for the marker of the water region. The maximal level of drawn contours was set to 8.
The threshold method in inverted threshold configuration type was used in the definition of markers for the air, hole, transducer and water regions. The threshold was set to a value of 0 in the definition of the markers of the air and water regions, to a value of 13 for the marker of the hole region and to a value of 10 for the marker of the transducer region. The maximum values of the grayscale were 1 in the definition of the air region marker, 255 for the markers of the hole and transducer regions and 2 for the marker of the water region. An erosion filter was used in the definition of the markers for the air and the water regions. A fixed ellipse was used as structuring element with size 8 × 8 for the air region and 9 × 9 for the water region.

Quantitative Evaluation of Image Segmentation and Ground Truth Images
F-measure was considered for the evaluation of the image segmentation quality [4], and is defined as where: α is the optimal detector threshold.
R is the recall. P is the precision.
The recall [25] is defined as where: T r are the pixels of the region in the ground truth image. S r are the pixels of the region in the segmented image. |T r | is the total number of pixels of the specific region. |T r ∩ S r | is the set of detected pixels that were assigned to the specific region correctly.
The precision is defined as where: |T r ∩ S r | are the detected pixels that were assigned to the specific region correctly. |T r \S r | are the pixels that were detected incorrectly for the region. These are the pixels that are present in S r but not in T r .
A MATLAB ® (2017 MathWorks, Inc., Natick, MA, USA) script developed inhouse was used to generate masks from the delineated regions of interest in each of the ground truth images. Masks for each region were generated from the segmentation and the F-measure was obtained per region.

Parametric Study of Image Segmentation Methods
To evaluate the impact of parameters on the resulting TLSM and WSAM segmentation the impact of different values on the F-measure were obtained (Table 1). A statistical analysis was carried out in order to find the image segmentation methods and their internal parameters that provide the most efficient segmentation of each region of the MR images. Table 1. Segmentation of the regions of interest with threshold and level set methods and parameters to involve in the parametric study along with their values. 1 The initial contour, for both the GAC and DRLSE methods, was given by the binary step function defined in Equation (10), so this parameter was not varied. Therefore, it was not considered in the parametric study. 2 The threshold used for the segmentation of the transducer region can be chosen between nine defined options between the two maximal peaks in the histogram of the Figure 6b. Each defined option is a threshold number as shown in Figure 10a. The selected threshold number to perform segmentation is Tn_transducer. The background of the image with the segmented transducer is the segmented water region.

Statistics
Data analysis was performed using R [26]. Data were tested for normality using the Shapiro-Wilk test [27]. The Bartlett's test [28] was used to test homoscedasticity. Data were analyzed using Welch's ANOVA [29] to confirm the existence of significant differences among the means of the measurement variables corresponding to the nominal variables. Welch's ANOVA was followed by post hoc Games-Howell test [30] to determine the groups of nominal variables that are significantly equal or different from each other. A p-value of less than 0.05 was considered to be significant.

Segmentation of the Air Region with Threshold Method (TLSM Approach)
In the following tables, the highest values for the F-measure median are reported. Table 2 contains information of transverse images and Table 3 contains information of sagittal images. The maximum is shown in bold at the top of the second column. Only the groups with different values for T_air and no significant difference (p > 0.05) are reported in the subsequent lines. The values for the F-measure medians are shown along with the values for quartile 1 (Q1), quartile 3 (Q3) and Interquartile Range (IQR). Table 2. Highest values for the F-measure median obtained in the segmentation of the air region with the threshold segmentation method in transverse images. Maximum is shown in bold.

Parameter Value
Highest

Segmentation of the Gel-Pad Region with GAC and DRLSE Methods (TLSM Approach)
In the following tables, the highest values for the F-measure median are reported. Table 4 contains information of transverse images; Table 5 contains information of sagittal images segmented with the GAC method and Table 6 contains information of sagittal images segmented with the DRLSE method. The maximum is shown in bold at the top of the second column. Only the groups with different values for u_0 and c_gelpad (GAC method) and phi_0 and a_gelpad (DRLSE method) and no significant difference (p > 0.05) are reported in the subsequent lines. The values for the F-measure medians are shown along with the values for Q1, Q3 and IQR. In transverse images the maximum was 0.9589 (0.9446-0.9698, IQR 0.0252). Given this maximum value, the corresponding maximum value of the F-measure median in the segmentation of the tissue region was 0.9133 (0.8427-0.9408, IQR 0.0981). In sagittal images the maximum was 0.9524 (0.9365-0.9610, IQR 0.0245). Given this maximum value, the corresponding maximum value of the F-measure median in the segmentation of the tissue region was 0.9264 (0.8889-0.9470, IQR 0.0581). Table 4. Highest values for the F-measure median obtained in the segmentation of the gel-pad region with the GAC segmentation method in transverse images. Maximum is shown in bold.

Parameters Values
Highest  Table 6. Highest values for the F-measure median obtained in the segmentation of the gel-pad region with the DRLSE segmentation method in sagittal images. Maximum is shown in bold and it was obtained with the DRLSE segmentation method.

Segmentation of the Transducer Region (TLSM Approach)
In the following tables, the highest values for the F-measure median are reported. Table 7 contains information of transverse images and Table 8 contains information of sagittal images. The maximum is shown in bold at the top of the second column. Only the groups with different values for a_gelpad (DRLSE method), c_gelpad (GAC method) and Tn_transducer (Threshold method) and no significant difference (p > 0.05) are reported in the subsequent lines. The values for the F-measure medians are shown along with the values for Q1, Q3 and IQR. In transverse images the maximum was 0.9305 (0.9213-0.9365, IQR 0.0153). Given this maximum value, the corresponding maximum value of the F-measure median in the segmentation of the water region was 0.9609 (0.8931-0.9768, IQR 0.0837). In sagittal images the maximum was 0.9323 (0.9221-0.9402, IQR 0.0181). Given this maximum value, the corresponding maximum value of the F-measure median in the segmentation of the water region was 0.9681 (0.9627-0.9715, IQR 0.088). Table 7. Highest values for the F-measure median obtained in the segmentation of the transducer region in transverse images. Maximum is shown in bold.

Segmentation Method and Parameter Value
Highest In the following tables, the values for the F-measure median are reported by region of interest. Table 9 contains information of transverse images and Table 10 contains information of sagittal images. The values for the F-measure medians are shown along with the values for Q1, Q3 and IQR. Table 9. Values for the F-measure median obtained in the segmentation of the regions of interest with the WSAM segmentation method in transverse images.

Region of Interest
Values for the F-Measure Median After the evaluation of the segmentation quality with F-measure using the WSAM approach for each region, a comparison with the maximum results of the TLSM approach (reported in Sections 3.1-3.3) was carried out.
In Figure 13, for each region, the obtained results with the WSAM and the TLSM approaches and the statistical significance between them are reported. For the TLSM approach parameters that yielded the maximum results are reported. Suitable methods and parameters for each region and plane are chosen according to the obtained value for the median and minimum of the F-measure.

Discussion
In the present research work, the problem of segmentation in planning MR images for the HIFU therapy was addressed. A set of methods and parameters to perform identification of objects in the planning MR images was analyzed in a parametric study supported by the evaluation of segmentation quality with F-measure and the analysis of means of the measurement variables with Welch's ANOVA.
For the segmentation of the air region, the results given by the threshold method yielded a greater median value than the WSAM in both the transverse and sagittal images. Moreover, there was a significant difference between the groups of results given by both the threshold method and the WSAM. This fact leads to consider the threshold method as the more convenient choice over the WSAM to perform the segmentation of this region.
For the segmentation of the gel-pad region, the results given by the WSAM did not present significant differences between the groups of results obtained with the GAC method in transverse images and with the DRLSE method in the sagittal images. Despite this fact, the results given by the WSAM yielded F-measure values that were higher than those obtained with the other methods. WSAM was then considered as a suitable choice to perform the segmentation of this region.
For the segmentation of the tissue region, the WSAM is the choice to carry out the segmentation of this region, given the fact that this segmentation strongly depends on the gel-pad region where WSAM was deemed the most appropriate method.
For the segmentation of the transducer region, the results given by the WSAM did not present significant difference with the group of results obtained with the DRLSE method in transverse images. However, the WSAM yielded a higher F-measure median value than the obtained with the DRLSE method, so the WSAM was considered the most appropriate choice to perform the segmentation. For sagittal images, the GAC method was the most suitable choice as it yielded a greater F-measure median value than the WSAM with a significant difference between the groups of results given by both methods.
The planning MR images for the treatment of abscesses with MRgHIFU in a murine model, used in the present research work, have specific characteristics such as a FOV of 120 × 120 × 48 mm, a pixel size of 0.5 mm and the presence of intensity inhomogeneity and five regions of interest. However, the planning MR images for other experimental setups may present differences. The dimensions of the FOV and the pixel size may be different from those used in the current experimental setup. Moreover, the tissue region may be larger if another animal model is used in the experimental setup, and the gel-pad region may be absent or if the target region is located a different anatomical area of the animal model. In addition, the problem of intensity inhomogeneity may be absent. These conditions may facilitate the application of the implemented segmentation methods and even better results could be achieved.
Level set methods (GAC and DRLSE) and WSAM proved to have enough capacity to solve the problem of segmentation in the current experimental setup, with functional results. However, there are limitations in the implemented segmentation methods such as the need to establish an initial contour for both the GAC and the DRLSE methods. Currently, this procedure is performed after the division of the MR image, an action that is performed by the user. Another limitation is the relative complexity in the procedure for the definition of markers of the WSAM. There is no limitation in the shape for the required markers for the WSAM, so simpler markers could be defined in a less complex procedure. Once this action has been achieved, an image segmentation solution based on a single method, the WSAM, would be available. However, the definition of markers is a procedure that previously requires the division of the MR image. The division of the MR image is a minor intervention that can be integrated as a simple step in a user interface in future developments.
The use of methods such as Scale Invariant Feature Transform (SIFT) [31] or Speeded-Up Robust Features (SURF) [32] could perform the detection of the regions of interest, based on the extraction of invariant features in the image. Both SIFT and SURF are invariant to rotation and scale. SIFT also has a remarkable performance if affine distortion and illumination changes are present. A specific training image is required to be defined for each region of interest in a given setup to perform detection. This approach would not require the division of the MR image, so it could be considered as an improvement of the segmentation procedure from the point of view of simplification of the overall process.

Conclusions
The present work was carried out to study the problem of segmentation in planning MR images for the HIFU therapy, using a murine model, to achieve an efficient identification of objects. A set of methods, with their parameters, to identify the regions of interest in the MR image were subjected to a parametric study. The findings of the study indicate that the segmentation of the regions of interest can be achieved with the proposed segmentation methods, given the characteristics of the MR images. However, the WSAM proved to be a remarkable method to carry out the identification of the regions of interest, although it is not the only one with which this task can be carried out. The threshold, the GAC and DRLSE methods can yield similar results to those obtained with the WSAM. However, some F-measure minimum values obtained with the GAC and DRLSE methods are close to zero.
This research work integrates preliminary results to generate more efficient procedures of image segmentation for treatment planning of the MRgHIFU therapy. Moreover, it also represents a contribution to the efforts that have been performed for the enhancement of the process of thermal ablation in the MRgHIFU therapy as an alternative for the non-invasive treatment of cancer.
There is no unique method that is able to perform the required segmentation task. The achievement of an efficient segmentation process is feasible, but it requires previous knowledge acquired from a preliminary segmentation step, application of filters or edge detection. However, the search for an automatic segmentation process remains. Funding: This work was partially supported through NSERC Discovery funds.