CN104657984B - The extraction method of three-D ultrasonic mammary gland total volume interesting image regions - Google Patents

Abstract

The invention belongs to the field of image processing, in particular to an automatic extraction method for a region of interest in three-dimensional ultrasonic breast volume imaging (ABVS). The present invention uses a method based on maximum directional phase information to process continuous cross-sectional two-dimensional images in three-dimensional ABVS images to obtain candidate regions of interest on each cross-sectional image; Remove irrelevant regions with prior knowledge such as gender and location characteristics; obtain shape and texture features for the remaining suspected tumor regions, input them into a binary logistic regression classifier to obtain the probability that each region may be a tumor, and select the region with the highest probability as a tumor Region: According to the selected region, the smallest ellipsoid containing the region of interest is obtained, which is the region of interest. The invention can realize the automatic extraction of the tumor region of interest in the three-dimensional ABVS image, obtain the accurate position of the tumor, reduce the workload of manual operation, and provide an important reference for further tumor detection.

Description

Automatic Extraction of Regions of Interest in 3D Breast Ultrasound Full Volume Images

technical field

The invention belongs to the technical field of image processing, in particular to an automatic extraction method for a region of interest in a three-dimensional ultrasonic breast full-volume image.

Background technique

Ultrasound imaging has important clinical applications due to its advantages of non-invasive, real-time, strong repeatability, and low cost. Compared with the traditional handheld two-dimensional ultrasound imaging, ABVS has a new imaging mode, which can standardize and automatically scan the breast, digitally process the image, and avoid individual differences among users; ABVS can scan the whole breast, compared with conventional ultrasound, The reconstructed coronal section is added, which can provide more information than two-dimensional images, and has good repeatability.

Since the volume of the tumor is relatively small compared to the entire ABVS image, the accuracy of direct tumor segmentation is low. Therefore, it usually requires the user to manually mark the central position of the tumor or the region of interest in hundreds of cross-sectional images for further analysis. This manual calibration method is very time-consuming and depends on user experience.

Aiming at this problem, the present invention proposes a method for fully automatic extraction of the ABVS region of interest. This method does not require the user to mark the tumor in advance, and can automatically find the region of interest, and finally obtain the smallest ellipsoid containing the region of interest. Applying the method of the present invention to an automatic analysis system of ABVS images can improve the accuracy of the entire system.

Contents of the invention

The purpose of the present invention is to propose a method for automatically extracting a region of interest in a three-dimensional ultrasonic breast full-volume image.

The present invention proposes a method for automatically extracting a region of interest in a three-dimensional ultrasonic breast full-volume image, the specific steps of which are as follows:

1. First, the ABVS image in DICOM format is reconstructed according to the distance between pixels in the three-dimensional direction, so that it corresponds to the size of the actual breast image; after reconstruction, 820 images of the transverse section and 750 images of the sagittal plane are obtained ; Depending on the scanning depth, 98~294 coronal image slices are obtained; on each slice, the distance between two adjacent images is 0.2 mm, and the distance between adjacent pixels on each slice image Also 0.2 mm;

2. According to the feature that the mammary gland is generally elliptical on the coronal plane image, calculate the minimum value mapping image of the first ten coronal plane images of the reconstructed ABVS image, and use the Hough ellipse transform to obtain the breast on the coronal plane template and apply it to all coronal image slices; use the relationship between the coronal and cross-sectional coordinates of the ABVS image to project the position of the breast on the coronal plane to the cross-section, determine the position of the breast in the cross-section, and initially reduce the cross-section The search range of the surface area of interest can remove the influence of interference such as noise and artifacts outside the mammary gland;

3. For the 820 cross-sectional images, every ten images were merged into a new image through minimum value mapping, and 82 images were obtained after merging; for each image, the method based on the maximum direction phase was used, combined with the tumor location characteristics and texture features, extract the candidate region of interest, and obtain the binary image of the corresponding region;

4. From step 3, a series of binary images of continuous cross-sectional candidate regions of interest are obtained; according to the continuity of the tumor, if the connected region of a certain slice does not overlap with the adjacent slices, it is called an irrelevant slice and should be removed first ; For the remaining slices, group them according to their continuity and the central position of the connected area, and each group represents a suspicious tumor;

5. Classify each group of cross-sectional image slices processed in step 4, extract the shape, texture and other features of the binary image and the corresponding grayscale image, and input them into the logistic regression classifier to obtain the probability that each group may be a tumor. The one with the highest probability is the cross-sectional image slice corresponding to the real tumor area;

6. In the binary images of the continuous slices of the tumor area determined in step 5, find the largest connected region, and determine the smallest ellipse containing the connected region; according to the coordinate correspondence between the sagittal plane and the cross-sectional image, find the abscissa corresponding to the center of the ellipse The sagittal plane image slice, using the method of extracting the candidate region of interest in step 3 to obtain the candidate region of interest in the sagittal plane image slice, and then determine the smallest ellipse containing the candidate region of interest on the sagittal plane; use the sagittal plane and The two ellipsoids of the cross section determine the length of the three main axes of the smallest ellipsoid that contains the region of interest in the ABVS image, thereby obtaining the region of interest of the tumor.

The relevant technical details involved in each step of the method of the present invention will be further specifically described below.

About step 1. The ABVS image is reconstructed from continuous image slices scanned by the breast automatic full-volume imaging instrument along the transverse direction of the breast. The image format is generally in DICOM format, and there are three orthogonal planes: transverse plane, sagittal plane and Coronal plane; the ultrasound images on the transverse and sagittal planes are similar to the two-dimensional ultrasound images obtained by traditional hand-held ultrasound equipment, and the general outline of the breast can be observed on the ultrasound images on the coronal plane; the DICOM file contains all the images Information, including the gray value of the image pixel and the pixel pitch of the image. Therefore, the 3D ABVS image can be reconstructed by reading the DICOM file information. The scanning depth of the ABVS image is generally 20-60 mm. According to the different scanning depths, the size of the reconstructed ABVS image is (98-294)×750×820. The schematic diagram of the three slices of the original ABVS image is shown in Figure 1, and the images of the reconstructed transverse, sagittal, and coronal planes are shown in Figure 2.

According to the relationship between the distance between the pixels of the ABVS image read from the DICOM file and the actual breast image, the ABVS image is reconstructed, and the distance between adjacent pixels is 0.2 mm.

About step 2. This step is to obtain the template of the mammary gland and eliminate interference such as artifacts and noises caused by poor contact around the mammary gland. The implementation method is as follows: firstly take 10 coronal images representing the surface layer of breast skin, calculate the minimum value mapping images of the first 10 coronal images, and then combine the ten images according to the principle of minimum value mapping. For the merged image, OSTU’s threshold value processing ^[1] method is used to obtain a binary image; the morphological dilation and erosion method is used to obtain the approximate outline of the mammary gland, and then the Hough ellipse transform method is used to detect and find the position of the ellipse corresponding to the mammary gland to generate A binary template image, the value inside the ellipse is 1, and the value outside the ellipse is 0; for the ABVS image, the ordinate of the coronal image corresponds to the abscissa of the cross-sectional image, so the coronal image can be used to generate The template map is used to determine the mammary gland area on the cross-sectional image, and remove the influence of noise and artifact interference outside the mammary gland.

In order to reduce the calculation amount of Hough ellipse transformation, according to the prior knowledge, the major and minor axis, direction and the ratio of the major and minor axes of the ellipse are limited. Here, the length of the major axis of the ellipse is limited to be between 400 and 800, the ratio of the major axis to the minor axis is greater than 0.75, and the angle between the major axis of the ellipse and the positive direction of the x -axis is between (0, π). The equation of the ellipse is:

(1)

Where a represents the major axis of the ellipse, b represents the minor axis of the ellipse, θ represents the angle between the major axis of the ellipse and the positive direction of the x -axis, and ( x ₀ , y ₀ ) is the center of the ellipse.

After detecting the ellipse containing the mammary gland on the coronal plane, use the relationship between the coronal plane and the three-dimensional coordinates of the cross-section of the ABVS image to project the position of the mammary gland on the coronal plane to the cross-section, determine the position of the mammary gland in the cross-section, and initially narrow down the cross-section of interest The search range of the candidate area removes the influence of interference such as noise and artifacts outside the mammary gland. The results of this part of the template extraction are shown in Figure 3. Using the template of the coronal plane to act on the whole ABVS image, the obtained images of the three slices are shown in Figure 4.

About step 3. The purpose of this step is to extract the region of interest of the cross-sectional image slice. Firstly, every ten consecutive cross-sectional images are merged by minimum value mapping to reduce the amount of computation. Then, the candidate regions of interest are extracted for each combined cross-sectional image.

The process of extracting candidate regions of interest in two-dimensional breast ultrasound tumor images is as follows:

1) Image preprocessing.

a) The speckle noise of the ultrasound image is serious, so anisotropic speckle noise removal (SRAD) is performed first. For the image shown in Figure 5(a), the filtered result is shown in Figure 5(b).

b) Reduce the range of the gray value ^[2] , and use the linear normalization formula (2) to adjust the gray value of the pixels in the image, and the result is shown in Figure 5(c).

(2)

Where L _n is the number of gray levels, lbound and ubound are taken as Q (0.5) and Q (0.95) respectively, and Q is the quantile function of the cumulative distribution of the histogram.

c) Enhance the hypoechoic area ^[3] , and use the adaptive Z-type equation represented by formula (3) to perform grayscale transformation:

(3)

Among them z _a , z _b , z _c determine the shape of the Z-type function. z _a and z _c determine the range of the nonlinear transformation of the curve, z _a is generally set to 20, z _c is the mean value of the image, and z _b determines the slope of the curve, which is obtained according to the slope of the gray distribution:

(4)

(5)

Suppose the image obtained after preprocessing is , as shown in Figure 5(d).

2) Find The maximum orientation phase map PMO . According to the method mentioned in literature [4], the image is filtered along six directions (0°, 30°, 60°, 90°, 120°, 150°) in the frequency domain of the image, and the final extraction with the maximum Phase information of the energy direction. These six directions cover the entire frequency spectrum (0~360°).

For image I , in order to calculate the PMO matrix, a two-dimensional log-gabor filter is first required to filter the image, where the transfer function of the two-dimensional log-gabor filter is a Gaussian function on a logarithmic scale, which is decomposed into two parts : radial filter and angular filter. For each direction θ ₀ , the constructed filter is:

(6)

The first part in the curly brackets represents the radial filter, and the latter part is the angular filter. ω is the radial coordinate; θ is the angular coordinate; ω ₀ is the center frequency of the filter, set four center frequencies of different scales, and their values are 1/3, 1/(3×1.7), 1/(3×1.7 ² ), 1/(3×1.7 ³ ); σ _θ determines the angular bandwidth centered on θ , which is set to 30°; Determines the radial bandwidth, usually set it to 0.55 to ensure the filtering effect while avoiding aliasing.

After the image is filtered by a total of 24 log-gabor filters with different center frequencies and different angles, it is transformed into 24 image matrices by inverse fast Fourier transform (FFT) ( x is the direction, s represents the center frequency scale), and then calculate the phase matrix LPA on each angle:

(7)

in Represents the phase matrix at an angle θ ; n represents the number of different center frequency filters; is the enhanced phase matrix at scale s , and its calculation formula is:

(8)

(9)

where imag ( x , s ) is The imaginary part of , real ( x , s ) is the real part of .

After filtering at six angles, six LPA feature matrices are obtained. Then, these six matrices are combined into one feature matrix. Here, the phase eigenvalue with the largest energy angle at a certain point is extracted as the final phase eigenvalue of the point, so as to obtain a new phase eigenvalue matrix. Because the local energy characterizes the structural information, the angle with the largest local energy is closest to the boundary direction. Therefore, the definition of the phase feature PMO at each pixel point is:

(10)

in Represents the energy matrix on the angle θ , and its calculation formula is:

(11)

Here imag ( x , s ) and real ( x , s ) represent the imaginary part and real part of the image after the log-gabor filter, respectively. ρ is the angle with the maximum local energy calculated according to formula (5), represents the superposition of the phase matrices of the center frequencies of the four scales over the angle ρ .

Through the above process, the maximum directional phase map PMO of the image I can be obtained. For the preprocessed image shown in Fig. 6(a), its maximum orientation phase map is shown in Fig. 6(b). Then, in order to highlight the tumor area and remove background noise, PMO was combined with 256- Multiply, and use a 5×5 median filter to filter the multiplied image, and the result is shown in Figure 6(c). However, after the maximum directional phase extraction, the gray value range of most pixels of PMO is [0,0.5], so the whole image looks blurry and dark. Use the following two formulas for grayscale correction:

(12)

(13)

In summary, the maximum orientation phase map of the grayscale image is obtained. In order to make the area with a large gray value represent the tumor area, the phase map of the maximum direction is reversed, and the final image is shown in Figure 6(d).

3) Due to the variety of tumor types, the texture and shape of the ultrasound images are also quite different. Some tumors have uneven internal gray levels, but the edges of the tumors have continuous phases; some tumors have relatively uniform internal gray levels, but the phases of the edges Indistinct or aliased with phase in the background. Therefore, in the selection of the candidate region of interest, it is necessary to combine the information of the grayscale image and the maximum orientation phase map. The preprocessed image I , the normalized PMO and PMO+I are subjected to threshold segmentation processing with the OSTU threshold method, and a binary image containing many connected regions is obtained. Remove the noise area in the image with an area less than 300 pixels and a connection rate with the edge greater than 50%. Calculate the energy function E of all remaining connected regions:

(14)

Among them, Compact represents the compactness of the connected region, and its value is the ratio of the area of the connected region to the area of the smallest convex shape; Area is the area of the connected region; Centdis is the distance between the center of the connected region and the center of the image; Eccentricity is the smallest connected region The eccentricity of the circumscribing ellipse; is the area and edge coincidence ratio, and its value is the ratio of the number of pixels where the area and the edge overlap to the total number of image edge pixels. w ₁ ~ w ₄ respectively represent the impact of the regional compactness, area, center distance and eccentricity on the regional energy function value, which is generally taken as 1, and can also be adjusted according to the specific ultrasound tumor image database. The back range is 1±0.5. w is the weight of the edge coincidence rate, and its value range is generally (0,0.5).

After the above process, the candidate regions of interest of the preprocessed image I , the normalized PMO and the PMO+I three images are respectively obtained. Calculate the energy function E of the candidate regions of interest in the three images, and select the region with the largest E as the final candidate region of interest in the image. Image I , normalized PMO and PMO+I three images and the results of threshold segmentation and region selection are shown in Fig. 7.

About step 4. After the previous step, a series of binary images of the candidate regions of interest in continuous cross-sections are obtained. The image slices of the tumor in the cross section are continuous, and the connected regions in the corresponding binary image overlap, as shown in Figure 8. In FIG. 8 , the candidate regions of interest of two adjacent image slices overlap each other, and the area of the candidate regions of interest first changes from small to large, and then from large to small. According to the continuity of the tumor, irrelevant image slices in which the connected area of the binary image does not overlap with the adjacent slices can be removed. The basis for judging whether a certain image slice is an irrelevant image slice is:

(15)

Where BW ^k represents the binary image of the kth cross-sectional image, and m and n are the number of pixels in the axial and lateral directions of the image, respectively. When R ( k )<50, it is considered that the candidate region of interest of the kth image does not overlap with the adjacent slice, and the kth image is removed from the image set. As an example, Fig. 9 shows three consecutive cross-sectional images of the ABVS image and the binary images of the corresponding candidate regions of interest, wherein the candidate regions of interest in the second image do not overlap with the two images before and after, It is therefore removed from the image set for this ABVS.

For the remaining image slices, first divide the continuous image slices into one group according to the continuity of the tumor area in the cross-sectional image slices; then in each group of continuous image slices, calculate the candidate regions of interest for adjacent slices The distance between centers, grouping image slices within 50 pixels into one group. In this way, each group obtained by the final division represents a serial cross-sectional slice of a suspicious tumor. In these groups, in addition to the real tumor area, the shadow behind the nipple and some hypoechoic areas caused by poor probe contact were also included.

About step 5. The purpose of this step is to extract the shape and texture features of the binary image and grayscale image corresponding to each group of continuous 2D images, and use these features for classification.

The shape and texture features of each image slice can be obtained by calculation, but because each group contains a different number of image slices, each feature obtained by each set of image slices has multiple values corresponding to the number of image slices, Therefore, it needs to be summarized to obtain a feature value that can be used to describe the overall characteristics of each group of images.

1) Shape features:

The aspect ratio of the candidate region of interest in the cross-sectional image within the group: generally speaking, the aspect ratio of the tumor is less than 1, while the aspect ratio of some non-tumor areas such as the shadow area behind the breast is greater than 1, so the aspect ratio of the region of interest can be Aspect ratio as a feature to distinguish tumors from non-tumors. This parameter is usually expressed by the ratio of the longest axial diameter L _short and the longest transverse diameter L _long of the tumor:

(16)

Calculate the aspect ratio Mabr of the _ROI of each image in the group, and take the median value as the aspect ratio of this group.

Overlapping ratio of regions of interest between consecutive image slices: It can be seen from Figure 8 that the area of interest regions of continuous cross-sectional images of tumors changes from small to large, and then from large to small. Generally speaking, the area in the middle position The image slice has the largest area of interest.

Define the region of interest overlap rate of two adjacent images:

(17)

where BW ^k represents the binary image of the kth image in the group. Equation (17) defines the overlap rate as the ratio of the overlapping area of the candidate region of interest in the kth image and the k +1th image to the area of the candidate region of interest in the k +1th image. In this way, for a suspicious tumor image group containing l images, we can get l -1 overlapping ratios. In order to summarize these l -1 features, their mean, standard deviation, product and gradient mean are respectively calculated as the features of classification.

2) Texture features:

Texture features include gray-scale co-occurrence matrix ^[5] features and gray-scale features.

Gray level co-occurrence matrix ^[5] features include entropy, contrast, correlation, energy and homogeneity. For the simplicity of calculation, first find the minimum circumscribed rectangle of the candidate region of interest in each binary image, and then expand the minimum circumscribed rectangle by 20 pixels vertically and horizontally to obtain a rectangle containing the candidate region of interest. Calculate the texture features of the grayscale images located in the rectangle, and then calculate the mean value of the grayscale co-occurrence matrix features of all image slices in the group as the features for classification.

The grayscale features include: a) the grayscale mean of the candidate region of interest and standard deviation ; b) the ratio of the gray mean value of the candidate region of interest to the background region within the rectangle containing the candidate region of interest Ratio to Standard Deviation , generally speaking, the gray mean value inside the region of interest is lower than the mean gray value of the background region, and the gray standard deviation is also lower than the gray standard deviation of the surrounding background region; In order to distinguish the hypoechoic area such as the shadow behind the nipple, the calculation formula of the acoustic shadow feature MeAP before and after the area of interest is:

(18)

Including:

(19)

(20)

(twenty one)

Here m and n are the number of pixels in the axial direction and horizontal direction of the image respectively, I ^h and BW ^h respectively represent the hth image in the 82 ABVS cross-sectional image slices and the binary image of its candidate region of interest, MBW represents the area The largest candidate area of interest, MeA is the average gray value in front of the suspicious tumor corresponding to the candidate area of interest in MBW , MeP is the average gray value behind the suspicious tumor corresponding to the candidate area of interest in MBW , Me is all candidate areas of interest in the suspicious tumor gray value of .

In summary, the final extracted features have five shape features and ten texture features. These features are input into the logistic regression classifier for classification, and the probability that each group of suspicious tumor regions may be tumors is obtained, and the one with the highest probability is selected as the tumor. In order to verify the validity of the method, ten times cross-validation method is adopted.

About step 6. The purpose of this step is to determine the lengths of the three main axes of the ellipsoid containing the region of interest through the acquired binary images and grayscale images of the continuous slices of the tumor region, so as to obtain an ellipsoid-shaped region of interest.

First, find the candidate region of interest with the largest area in the middle of the tumor area, expand the candidate region of interest by 10 pixels using the morphological expansion method, and calculate the minimum circumscribed ellipse of the expanded region, and the major axis of the ellipse The length is a and the minor axis is b . It can be known from Figure 1 that the coordinates of the abscissa on the cross section correspond to the serial number of the sagittal plane image slice, and the coordinates of the center point of the ellipse circumscribing the region of interest are ( x ₀ , y ₀ ), then the x _0th sagittal plane image is taken Slice, in the region where the y -axis coordinates of the sagittal plane image are ( y ₀ - b /2, y ₀ + b /2), use the above-mentioned method of extracting the candidate region of interest based on the maximum direction phase to obtain the sagittal plane of interest Candidate area. Calculate the circumscribed ellipse and its major axis length c of the region after the candidate region of interest is extended outward by 10 pixels. In this way, the lengths a , b and c of the three main axes of the ellipsoid containing the region of interest are obtained.

The equation of the ellipsoid is:

(twenty two)

Where ( x ₀ , y ₀ , z ₀ ) are the coordinates of the center point of the ellipsoid, and a , b , and c are the axis lengths of the ellipsoid in three orthogonal directions. The schematic diagram of the ABVS image and its ellipsoid ROI is shown in Figure 10, and the final ROI of the actual ABVS image is shown in Figure 11.

Description of drawings

Figure 1. Schematic diagram of cross-section, sagittal plane and coronal plane of ABVS image, where Horizontal plane represents a schematic diagram of cross-section, Sagittal plane represents a schematic diagram of sagittal plane, and Coronal plane represents a schematic diagram of coronal plane.

Figure 2. Transverse, sagittal, and coronal views of the reconstructed ABVS image. A breast tumor is contained within a square wireframe.

Figure 3. Coronal template acquisition, the extracted template is inside the oval wireframe.

Figure 4. Applying the ellipse template to limit the search range of the region of interest: (a) Three slices of the original image, and inside the wireframe are images of each slice after using the ellipse template; (b) Three slices after applying the ellipse template.

Figure 5. Cross-sectional image preprocessing. Among them, (a) original image; (b) SRAD filtered image; (c) image with reduced gray scale; (d) enhanced image of hypoechoic area.

Figure 6. Maximum orientation phase map of cross-sectional images. Among them, (a) original image; (b) PMO image; (c) PMO × (256- I ) median-filtered image; (d) grayscale-adjusted image.

Figure 7. Extraction of regions of interest. Among them, (a), (b), (c) are I , PMO , I+PMO respectively; (d), (e), (f) are the OSTU binary image of three images respectively; (g), ( h), (i) are the candidate regions of interest obtained from the three images respectively; the final candidate region of interest is inside the rectangular wireframe.

Figure 8. Continuous images of tumor regions and their corresponding candidate regions of interest. The left side is the continuous original image slice, the right side is the corresponding binary image, and the white area is the candidate area of interest.

Figure 9. Candidate regions of interest for consecutive images of non-tumor regions. The left side is the continuous original image slice, the right side is the corresponding binary image, and the white area is the candidate area of interest.

Figure 10. A schematic diagram of an ABVS image and its ellipsoidal region of interest. Among them, (a) cross-sectional ROI map; (b) sagittal ROI map; (c) coronal ROI map.

Figure 11. Extraction results of the tumor region of interest: (a) images of three slices of the largest tumor region of interest; (b) images of three slices of the extracted region of interest.

Detailed ways

The method for automatically extracting the region of interest in the three-dimensional ultrasonic breast volume imaging (ABVS) proposed by the present invention is tested. The ABVS images were taken from an ACUSON S2000 ^TM ultrasound instrument from Siemens. The system is equipped with a broadband linear probe (14L5BV), which can obtain breast volume images of 15.4 cm×16.8 cm×(2~6) cm. A total of 15 ABVS images were collected in this experiment, each with 98-294 coronal images, 820 transverse images, and 750 sagittal images.

First, the original ABVS image is reconstructed, and the three sections (transverse, sagittal, and coronal) of the reconstructed image are shown in Figure 2, where the general outline of the breast can be seen on the coronal image. It can be seen from Fig. 2 that the breast contour is generally close to an ellipse, so the Hough transform is used to find the ellipse representing the breast on the coronal image, as shown in Fig. 3 . In the following processing, only the area inside the ellipse needs to be considered, which can reduce the amount of calculation and eliminate the interference of noise, artifacts, etc. outside the breast. Figure 4 shows the original image of each slice and the image of each slice after applying the ellipse template.

Figure 5 shows the process of preprocessing the cross-sectional images. Figure 5(a) shows the original image. In view of the characteristics of low contrast and severe speckle noise in ultrasound images, anisotropic speckle noise is firstly eliminated, and the result is shown in Fig. 5(b). Then, reducing the grayscale range actually stretches the grayscale value between 0 and 255, and the result is shown in Figure 5(c). Finally, the gray value of the hypoechoic region is enhanced, and the result is shown in Figure 5(d).

Figure 6 shows the maximum orientation phase map of the cross-sectional images. Since some tumors have large differences in the internal gray value, if only the gray image is used for threshold segmentation, the complete tumor area cannot be obtained. Therefore, the phase information is added. Fig. 7 shows the extraction of the ROI of the three images of I , PMO and I+PMO and the candidate ROI selected from the three images. Here, the last selected candidate region of interest is obtained by performing threshold segmentation and region selection on the image I+PMO .

Fig. 8 and Fig. 9 are cross-sectional continuous images of non-tumor area and tumor area and corresponding candidate regions of interest, respectively. It can be seen from this that the regions of interest in the image of the tumor region are continuous and overlap each other; the regions of interest in the image of the non-tumor region do not overlap with each other. According to this feature, a part of irrelevant image slices can be removed.

Figure 11 is the extraction result of the tumor region of interest. The images that remove irrelevant slices are grouped, and then the texture and shape features of each group of image slices are calculated, and the logistic regression classifier is used to classify, and the validity of the ten-fold cross-validation method is verified.

In the experiment, a total of 15 cases of ABVS images were processed, of which 14 cases could get the group containing the tumor area, and the accuracy rate reached 93.3%. Using the grayscale image and binary image of the image slices in the group containing the tumor area, combined with the sagittal plane image, calculate the length of the principal axes of the smallest ellipsoid containing the tumor area in three orthogonal directions to obtain the ellipsoid shape of the tumor The region of interest is shown in Figure 10. Figure 11(a) shows three slices of the ABVS image, and the finally obtained three slices of the actual region of interest are shown in Figure 11(b). The region of interest is represented by an ellipsoid, and the obtained ellipsoid is closer to the real tumor region than the cubic region.

In summary, the present invention is suitable for the extraction of tumor regions of interest in 3D breast full-volume images. The entire extraction process is completely automatic and accurate, and can realize the positioning of tumor regions without relying on the user's labeling, which is more objective and efficient.

references

[1] N. Otsu, "A threshold selection method from gray levelhistograms," IEEE Trans. Syst., Man Cybern., vol. 9, pp. 62–66, 1979.

[2] B. Liu, H.D. Cheng, J.H. Huang, J.W. Tian, X.L. Tang, J.F. Liu, “Fully automatic and segmentation-robust classification of breast tumors based on local texture analysis of ultrasound images,” Pattern Recognit., vol.43, pp. 280–298, 2010.

[3] M. Xian, Y.T. Zhang, H.D. Cheng, “Fully automatic segmentation of breast ultrasound images based on breast characteristics in space and frequency domains,” Pattern Recognit., vol.48, pp. 485–497, 2015.

[4] J. Shan, H.D. Cheng, Y. Wang, “Completely automated segmentation approach for breast ultrasound images using multiple-domain feature,” Ultrasound in Med. & Biol., vol.38, pp. 262–275, 2012.

[5] R. M. Haralick, K. Shanmugam, I. H. Dinstein, “Textural features for image classification,” IEEE Trans. Syst., Man Cybern., vol.3, pp. 10–621, 1973.

Claims (6)

1. an automatic extraction method of region of interest in three-dimensional ultrasonic breast full volume imaging (ABVS), it is characterized in that concrete steps are: (1) First, the ABVS image in DICOM format is reconstructed according to the distance between pixels in the three-dimensional direction, so that it corresponds to the size of the actual breast image; after reconstruction, 820 images of the transverse section and 750 images of the sagittal plane are obtained According to different scanning depths, 98-294 coronal image slices were obtained; (2) According to the feature that the mammary gland is generally elliptical in the coronal image, calculate the minimum value mapping image of the first ten coronal image slices of the reconstructed ABVS image, and use the threshold value processing method for the minimum value mapping image, and use Hough Ellipse transformation obtains the template of the breast on the coronal plane and applies it to all coronal plane images; using the relationship between the coronal plane and the three-dimensional coordinates of the cross-section of the ABVS image, the position of the breast on the coronal plane is projected to the cross-section to determine The location of the breast in the cross-section, initially narrowing the search range of the region of interest in the cross-section, and removing the influence of noise and artifacts outside the breast; (3) For the 820 images of the cross section, every ten images are merged into a new image through the minimum value mapping, and 82 images are obtained after the merger; image enhancement is performed on each image based on the maximum energy direction phase method, Combined with tumor location characteristics and texture characteristics, extract candidate regions of interest, and obtain binary images of corresponding regions; (4) From step (3), a series of binary images of continuous cross-sectional candidate regions of interest are obtained; according to the continuity of the tumor, if the connected region of a certain slice does not overlap with the adjacent slices, it is called an irrelevant slice. Remove them first; group the remaining slices according to their continuity and the central position of the connected area, and each group is represented as a suspicious tumor; (5) Classify each group of cross-sectional image slices in step (4); extract the shape and texture features of its binary image and corresponding grayscale image, and input it to the logistic regression classifier to obtain the probability that each group may be a tumor, the probability The largest one is the cross-sectional image slice corresponding to the real tumor area; (6) In the binary image of the continuous slice of the tumor area determined in step (5), find the largest connected region as the candidate region of interest, and use the morphological expansion method to expand the region outward by 10 pixels, and determine the region containing The minimum ellipse of the region after the expansion; find the sagittal plane image sheet corresponding to the abscissa of the ellipse center, adopt the algorithm of extracting the candidate region of interest in the step (3) to obtain the candidate region of interest of the sagittal plane image sheet, and use the same morphology The expansion method expands the area outward by 10 pixels, and then determines the smallest ellipse containing the expanded area; using the two ellipses of the sagittal plane and the cross section, determines the minimum ellipsoid of the ABVS image containing the region of interest The three main axes are long to obtain the tumor's region of interest. 2. the automatic extraction method of region of interest in three-dimensional ultrasonic breast full volume imaging (ABVS) according to claim 1, is characterized in that in step (1) and step (2), ABVS image is along mammary gland cross-sectional direction by instrument The serial image slices obtained by scanning are reconstructed in DICOM format and have three orthogonal planes: transverse, sagittal and coronal planes; the ultrasound images in the transverse and sagittal planes are compared with those acquired by traditional hand-held ultrasound equipment. Two-dimensional ultrasound images are relatively similar, and the general outline of the breast can be observed on the coronal ultrasound image. 3. the automatic extraction method of region of interest in the three-dimensional ultrasound breast full volume imaging (ABVS) according to claim 1, it is characterized in that adopting the interested candidate based on maximum energy direction phase and tumor prior knowledge in step (3) The region extraction method extracts the candidate region of interest from the cross-sectional image slice, and the specific process is as follows: (1) Perform image noise reduction and contrast enhancement preprocessing on ultrasound images: first use anisotropic noise reduction method to reduce the speckle noise of ultrasound images, and then improve image contrast by grayscale stretching transformation; (2) Calculate the maximum energy direction phase map PMO for the image obtained by the above preprocessing; since in the ultrasound breast image, the direction of the tumor boundary is constantly changing, along six directions in the frequency domain of the image: 0°, 30°, 60°, 90°, 120°, 150°, filter the image, and finally extract the phase information with the maximum energy direction; these six directions cover the entire frequency spectrum from 0 to 360°; For image I, calculate the PMO matrix. Here, a two-dimensional log-gabor filter is used to filter the image. The transfer function of the two-dimensional log-gabor filter is a Gaussian function on a logarithmic scale, which is decomposed into two parts: path direction filter and angle filter; for each filtering direction θ ₀ , the constructed filter is: Among them, the former part in the curly brackets represents the radial filter, the latter part is the angular filter, ω is the frequency of the radial filter, θ is the angular coordinate, ω ₀ is the center frequency of the radial filter, and the centers of four different scales are set Frequency, its values are 1/3, 1/(3×1.7), 1/(3×1.7 ² ), 1/(3×1.7 ³ ); σ _θ determines the angular bandwidth centered on θ, Determine the radial bandwidth to ensure the filtering effect while avoiding aliasing; After the image is filtered by a total of 24 log-gabor filters with different center frequencies and different angles, it is transformed into 24 image matrices I'(x,s) by inverse fast Fourier transform (FFT), where x is the direction, s represents the center frequency scale, and then calculate the phase matrix LPA at each angle: Among them, LPA _θ represents the phase matrix at the angle θ, n represents the number of filters with different center frequencies, phase'(s) is the enhanced phase matrix under the scale s, and its calculation formula is: phase'(s)=[π-πcos(2phase)]/2 (3) Where imag(x,s) is the imaginary part of I'(x,s), and real(x,s) is the real part of I'(x,s); After filtering through six angles, six LPA feature matrices are obtained; then, these six matrices are combined into one feature matrix; here, the phase eigenvalue with the largest energy angle at a certain point is extracted as the final phase eigenvalue of the point, Thus a new phase feature matrix is obtained; then, the phase feature PMO is defined at each pixel as: PMO(i,j)=LPA _ρ (i,j),ρ=argmax Eng _θ (i,j) (5) Where Eng _θ represents the energy matrix at angle θ, and its calculation formula is: ρ is the angle with the maximum local energy calculated according to formula (5), and LPA _ρ represents the mean value of the four-scale center frequency phase matrix at the ρ angle, thus forming the maximum direction phase map PMO; (3) In the selection of the candidate region of interest, combined with the information of the grayscale image and the maximum energy direction phase map, the preprocessed image I, the normalized PMO and the PMO+I three images were respectively processed by the OSTU threshold method. Perform threshold segmentation to obtain a binary image containing many connected regions; remove noise regions in the image whose area is less than 300 pixels and whose connectivity rate with the edge is greater than 50%, and calculate the energy function E of all remaining connected regions: Among them, Compact represents the compactness of the connected region, and its value is the ratio of the area of the connected region to the area of the smallest convex shape; Area is the area of the connected region; Centdis is the distance between the center of the connected region and the center of the image; Eccentricity is the smallest connected region The eccentricity of the circumscribed ellipse; OvlpRate is the area and edge coincidence rate, and its value is the ratio of the number of pixels overlapping the area and the edge to the total number of image edge pixels; w ₁ ~ w ₄ respectively represent the compactness of the area, The influence of area, center distance and eccentricity on the area energy function value, w is the weight of edge coincidence rate; Through the above process, the candidate regions of interest of the preprocessed image I, normalized PMO and PMO+I three images are respectively obtained; the energy function E of the candidate regions of interest of the three images is calculated, and the region with the largest E is selected as The final candidate region of interest for the image. 4. the automatic extraction method of region of interest in the three-dimensional ultrasonic breast full volume imaging (ABVS) according to claim 1, is characterized in that in the step (4), judges whether a certain image sheet is irrelevant image sheet basis is: Among them, BW ^k represents the binary image of the kth cross-sectional image, m and n are the number of pixels in the axial direction and lateral direction of the image respectively; when R(k)<50, it is considered that the candidate region of interest in the kth image is the same as There is no overlap between adjacent slices, and the kth image is removed from the image set; For the remaining image slices, first divide the continuous image slices into one group according to the continuity of the tumor area in the cross-sectional image slices; then in each group of continuous image slices, calculate the candidate regions of interest for adjacent slices The distance between the centers, the image slices with a distance within 50 pixels are grouped into one group; in this way, each group finally divided represents a continuous cross-sectional slice of a suspicious tumor. 5. the automatic extraction method of region of interest in three-dimensional ultrasonic breast full volume imaging (ABVS) according to claim 1, is characterized in that in step (5), extracts the texture and the shape feature of each group of image slices, input to Logistic regression classifier to obtain the probability that the group may be a tumor; the group with the highest probability represents the tumor area; Among the extracted features, texture features include gray-scale co-occurrence matrix features and gray-scale features; shape features include overlap rate features and aspect ratio features between slices within a group. 6. the method for automatically extracting the region of interest in three-dimensional ultrasonic breast volume imaging (ABVS) according to claim 1, is characterized in that in step (6), by the binary image of the continuous slice of the tumor region of acquisition and gray Degree image, determine the length of the three main axes of the ellipsoid containing the region of interest, so as to obtain the region of interest in the shape of the ellipsoid, the specific process is: First, find the candidate region of interest with the largest area in the middle of the tumor area, expand the candidate region of interest by 10 pixels using the morphological expansion method, and calculate the minimum circumscribed ellipse of the expanded region, and the major axis of the ellipse The length is a, and the length of the short axis is b; the coordinates of the abscissa on the cross section correspond to the serial number of the sagittal image slice, and the coordinates of the center point of the ellipse circumscribing the region of interest are (x ₀ , y ₀ ), then take the x _0th A slice of a sagittal plane image, in which the y-axis coordinates of the sagittal plane image are in the region of (y ₀ -b/2, y ₀ +b/2) using the above method for extracting candidate regions of interest based on the maximum energy direction phase Obtain the candidate region of interest on the sagittal plane; calculate the circumscribed ellipse and the major axis length c of the region after the candidate region of interest is extended outward by 10 pixels, so that the lengths of the three principal axes of the ellipsoid containing the region of interest are obtained a, b and c.