AU2021100714A4

AU2021100714A4 - Robust MRI Segmentation Using Background Noise Removal with Polyfit Surface Evolution

Info

Publication number: AU2021100714A4
Application number: AU2021100714A
Authority: AU
Inventors: Changjiang Liu
Original assignee: Sichuan University of Science and Engineering
Current assignee: Sichuan University of Science and Engineering
Priority date: 2021-02-04
Filing date: 2021-02-04
Publication date: 2021-05-06
Anticipated expiration: 2029-02-04

Abstract

The invention discloses a robust MRI segmentation method to outline potential abnormality blobs. Thresholding and boundary tracing strategies are employed to remove background noise, and hence the Rol in the whole process is set. Subsequently, a polyfit surface evolution is proposed to approximately estimate bias field, which makes segmentation robust to image noise. Simultaneously, customized initial level set functions are devised so as to detect subtle bright or dark blobs which are highly potential abnormality regions. The proposed method improves bias field estimation and level set method to acquire fine segmentation with low computational complexity. (a) (b) (c) (d) (e) (f (g) (h) (i) Figure 5

Description

(a) (b) (c)

(d) (e) (f

(g) (h) (i) Figure 5

BRIEF DESCRIPTION OF THE FIGURES

Figure 1 is the diagram of the Bias field polyfit estimation in embodiments of the invention;

Figure 2 is the diagram of customized initial contour of $',# in embodiments of the

invention;

Figure 3 is the diagram of energy functional values in the iterations in embodiments of the

invention;

Figure 4 is the effect diagram of background noises removal and Rol based segmentation in

embodiments of the invention;

Figure 5 is an image segmentation effect diagram of embodiments of the invention for bright and

dark blobs detection.

DESCRIPTION OF THE INVENTION

The technical scheme in embodiments of the invention is described clearly and completely

combining the appended figures in the novel embodiments, but the scope of protection of the

invention is not limited by the specific mode of implementation.

An improved multi-phase level set method for MRI segmentation with background noise

removal and polyfit surface.

<1> Multi-phase level set method image segmentation, that is, the model of an MR imaget

has been widely accepted as:

I(x) = b(x)J(x) + n(x) (1)

Among them, b(x) is bias field that accounts for the intensity inhomogeneity in the

observed image I(x), n(x) is the additive Gaussian noise with zero-mean, and J(x) is the true

image without bias field and noise. Assume J(x) is approximately piece wise constant, and ci

for the i-th region.

The level set formulation is given as follows: k F(4, c, , b) + v 2i L (0j) + y R (#P) (v, p 2 0) (2)

In (2), the energy functional E can be defined as:

E((D, c, b) = f (E 1 f K(y - x)|II(x) - b (y) ci l2 dy)Mi ((D(x))dx (3)

among them, (= (#1,...,#k) is a vector valued function, C = (c 1 ,..., CN), K(y - x) is a

non-negative window function, being K(y - x) = 0 for x V Oy, OY = {x:Ix - y| : p}. For the

case of N 3, two or more level set functions 01, ., are used to define membership

functions Mi(i = 1, . . , N), given by:

(1, x E El M0, otherwise

In (2), the energy term L(#) computes the length of the zero level set of $/ in the

conformal metric and it is defined by:

f(fj x) IVH(#I(x))Idx (4)

among them, V is the vector differential operator, namely gradient operator, H is the heaviside

function.

In (2), the energy term R(#1 (x)) is introduced to formulate level set evolution without

re-initializations. R(#j(x)) is given by:

R (#fI(x)) = f (IV| (x) I- 1)2 dx (5)

The energy functional F(, c, b) is minimized iteratively with respect to each of its

variables D, c, b, given that the other two have been updated in the previous iterations.

1) Optimization of $/(x):

(___ _ -_ N dMi(4b) a e + v(# )div( ____ )

+p[A#j - div(V)] (6) I2 Iv-PjI ei (x) (1*K)- 2ciI-(b*K)+c(b2 *K) (j = 1,2, . . , k)

among them, - calculates the per-element product of two matrices, * is the convolution operator,

1 is a matrix of ones with the same size asI, (.)is the derivative of the heaviside function and div(.) is the divergence.

2) Optimization of c:

c. f (b2*K)Mg(D(y))dy (b*K)IMi(D(y))dy (i = 1, N) (7) f

3) Optimization of b:

b- IJ(')*K (8) J(2)*K

among them,

N

J ~l = ciMi(<D(y))

N

j(2) cim,(<D(y)) i=1

<2> Background noise removal (Rol setting)

MRI is vulnerable to noise contamination. If noise in background pixels is severe, it will influence final curve evolution using multi-phase level set method. While the gray values of MRIs differ greatly in background and organs. Thus, thresholding is a simple but efficient method to remove the background from MRI. In our method, we first extract coarse regions of a target using Otsu's method. Subsequently, we retrieve connected components from the binary image and select the contour outlines in the image whose areas are greater than a predefined A value (Step 3 below). Finally, we fill the regions bounded by the preceding selected contours, which serve as the Rol in this paper. Assuming gray values of background are much lower than the foreground ones, the procedure can be outlined as follows:

1) Compute the optimal threshold value using Otsu's algorithm and denote it by T.

2) Apply a fixed-level threshold to each array element of image f and obtain the binary image g:

(255, f(i, j) > iT g(i,) D o0, otherwise

among them, constant 0 < ij 1 can guarantee that the regions of foreground with lower gray values are not ignored. Based on empirical observations this constant commonly takes values in the interval [0.7, 0.9].

3) Trace only the outer contours of image g. For the i-th contour, if the area is greater than

A, fill the area bounded by the contour with the gray value 255 and denote the corresponding mask

image by Bi. The constant A can skip small regions incurred by noise. The union set of such Bi,

denoted by B = UBj, is the final ROI image, with the notation R in which the value1 represents

the points to be further processed:

R (i,j) = 1, B(ij) = 255 (10) ) 0, otherwise

<3> Multi-phase level set method with bias field polyfit estimation

Local neighboring pixels are contributive with the definition of the window function K.

Confined to a local domain, the b can be approximately represented as a linear combination of

simple functions:

b(y) = w T G (y) (11)

In this invention, we employ the basis of polynomials with degree p that fits the bias field. It

is presented as:

G(y)= {xjy'-1|0 ! i O p, 0 j i}. (12) Substituting (11) into (3), the energy functional E is converted into:

E(W, c, D) = f Exdx (13)

among them, Ex = >N f K(y - x)|I(x) - WTG(y) 12 Mi((x))dy.

Therefore, - -2v + 2Aw, where OW

v f G(y)I(x)(I ciMi(D(x))K(y - x))dydx

Sf I(x) I ciMi(D(x))(f G(y)K(y - x)dy)dx A= f f G(y)G T(y)( c Mi((x))K(y - x))dydx

f c 2M(CD (x))(f G(y)G T(y)K(y - x)dy)dx

Applying the convolution operation, the equations above will be given by: v f (G * K)I(x)E ciMi(CD(x))dx i=1 (14) A = f (GG T * K) cN CM, (CD(x))dx

Note that v E RLx1, A E RLxL. Let O = 0, we get: aw

w A-1v (15)

and,

b(y) wT G(y) (16)

Bias field polyfit estimation can retain the smoothness, compared to existing literatures, one

of the differences in our approach is the additional convolution operations in Eq. (14), which help

localize and remove noise.

<4> Customized initial level set functions to detect abnormalities

One of the disadvantages of traditional level set methods is that it requires manual placement

of a closed curve near the desired boundary. Moreover, some abnormalities occur as negligible

patches, therefore no suitable initial level set functions can result in such abnormalities being

detected. As a result, we introduce level set functions which intersect the target region in MRIs as

much as possible. These initial level set functions help us find out the abnormalities automatically,

especially the negligible patches.

In order to segment more regions of interest, such as bright and dark blobs, we formulate

initial level set functions so as to make initial contours scatter linearly, equally spaced in image

domain. Classification of anatomical structure appearances in a MRI slice of the brain, suggests

that the Rol mentioned in this invention can be classified into light, gray and dark parts; namely,

two level set functions constructing 3 membership functions (k = 2, N = 3). The pseudo-code for

definition of 0'(x, y), 0'(x, y) is shown below, where the image dimension is W x H, p > 0 is

a constant with any value, and the Rol is bounded by a rectangle {(x,y)IXmin X X ! xmax,Ymin

Y 3 Ymax):

1) Initialization: #0(x, y) = p, 0 (x, y) = p, x = 0: 1: W - 1, y = 0: 1: H - 1.

Among them, the expression index = initVal:step:endVal means increment the index variable from initVal to endVal by step, and repeat execution of statements until index is greater than endVal.

2) Evaluation 1: #P(x,y) = -p, P (x + 1, y) = -p, x = xmin: 2: xmax, y=ymin: 2: Ymax.

3) Evaluation 2: 0 (x, y) = -p , 2(x + 1, y) = -p , X = Xmin: 2: Xmax, Y = ymin

+ 1: 2: Ymax.

4) Rol clipped: #1 = 1 - R, 02 = 02 - R.

In this case, the relationship between membership functions and level set functions is given by:

M 1(x) = 1 - H(#1(x)) M 2 (x) = H(#1(x))H(#2(x)) (17) M 3 (x)=)(1 - H(#2(x)))

among them, 0 Mi(x) 1 and M = 1.

<5> ROI based numerical implementation

In numerical practice, the Dirac function 6(x) in (7) is smoothed as:

S (x) = er 2+x2 (18) (8

Iteration for #P is proposed by:

#j+1 = #j + r (19) at

In this paper, the window function K is taken as a Gaussian function with a standard deviation a, denoted by K,.Our proposed Rol based numerical implementation is presented as follows:

1) ROI mask. I = I -R.

2) Initialization. Set the following parameters: N, il, p, o, v, y, e and -. Initialize b and

01,..,(Pk. Due to their invariability during iterative procedures, we prepare G, GGT, 1**K,

(1 * K,) . I 2 , (G*K,).I and GGT * K,.

3) Iterative procedure. For the i-th iteration, first update c via (7). Second, update <D via

(18), (6) and (19). Finally, update b via (14), (15) and (16).

4) ROI mask. #j = #j . R for next iteration.

5) Redo Step 3, until i exceeds the predefined maximal iteration.

The results of the experiment are as follows:

The invention mainly considers the segmentation of the MRI, and we analyzed 112 frames of

MR brain images, detailed information can be seen in Table 1. The algorithm was implemented

using the Microsoft Visual Studio development platform with the Open CV (Computer Vision)

library, and Matlab. In the experiments, the number for the maximal iteration needed is 50. If not

specified, the parameters used in the experiment are as follows: N = 3, rj = 0.7, p = 4,

4.0, v = 4.0, p = 1.0, E = 1.0, T = 0.1.

Table 1 Parameters for MRI of experiment Width 96 Modality MR Height 192 Repetition Time 36s Bit Depth 12 Echo Time 9.2s Slice Thickness 1 mm Imaging Frequency 63.6250

Pixel Spacing (1.0417mm, Protocol Name COR3DPreSc 1.0417mm) Norm OFF ' S'

Embodiment 1, referring to Figure 1, picture (a) is bias estimation by traditional method and

picture (b) is bias field polyfit estimation proposed. Due to the additional convolution operations of

our proposed algorithm compared to existing literatures, it helps localize and remove noises of

MRI.

Embodiment 2, referring to Figure 2, picture (a) is Rol based initial contour, and picture (b) is

the partial enlarged view. Here, red and blue contours are a visual presentation of the level set

functions #0 and 0 respectively. The initial level set functions guarantee that initial contours

intersect the tissues extremely, merge or split to evolve to the final desirable contour.

Embodiment 3, referring to Figure 4, picture (a) is inaccurate curve evolution by traditional

method with spurious contours incurred by background noise, picture (b) is open contours with

some faulty contours based on background noise removal and no Rol localization, and picture (c)

is the accurate contour by proposed algorithm in this invention. Therefore, the invention can not

only obtain accurate final contours, but also can speed up calculation.

Embodiment 4, referring to Figure 5, picture (a) is source image of frame 88, picture (b) is

Rol, picture (c) is the curve evolved after 10 iterations, picture (d) to picture (f) are final contours

and segmentation results, and picture (g) to picture (i) are three split parts of segmentation results.

Among them, picture (h) is bright blob detected and picture (i) is dark blob detected in MRI, which

are high likelihood of potential abnormalities. Therefore, our method can extract accurate contours

of MRI, specifically bright blobs and dark blobs simultaneously.

To sum up, the invention can not only segment MRI, but also detect subtle bright or dark

blobs automatically. Background removal helped us to exclude background noise, which also

provided Rol to speed up successive processing. Polyfit level set method with customized initial

level set functions in this invention does segment MRIs into individual parts, including tiny tissues.

It has good performance in extracting all the blobs in MRIs, in which there are potential lesions.

The above disclosure is only a specific embodiment of the invention, however, the

embodiment of the invention is not limited to thereof, and any change which the technician in this

field can think of shall fall within the scope of protection of the invention.

EDITORIAL NOTE 2021100714

There are 6 pages of claims only.

Claims

1. An improved multi-phase level set method for MRI abnormality detection with background

noise removal and polyfit surface is characterized by:

has been widely accepted as:

I(x) = b(x)J(x) + n(x) (1)

for the i-th region.

The level set formulation is given as follows:

k 'F(4, c, b) = E(D, c, b) + v E _1 L (0j) + y R (j) (v, p 2 0) (2) ,j=1

In (2), the energy functional E can be defined as:

E((D, c, b) = f (E 1 f K(y - x)|II(x) - b (y) ci l2 dy)Mi ((D(x)) dx (3)

among them, (= (#1, .. ,# ) is a vector valued function, C = (c 1 ,..., CN), K(y - x) is a

functions Mi(i = 1, . . , N), given by:

(1, x E El M0, otherwise

conformal metric and it is defined by:

f(fj x) IVH(#I(x))Idx (4)

function.

re-initializations. R(#P(x)) isgivenby:

R (f(x) (fV(|7 (x)| - 1)2 dx (5)

The energy functional F(,c,b) is minimized iteratively with respect to each of its

1) Optimization of $/(x):

N asbj- aM(ej + v5(#j)div( j at 1 = (j '0-1(j 6 +p[A#j - div(v~J)] (6) Iv$pjI ej(x) =2 - (1* K) - 2cil - (b * K) + c?(b2 * K) (j = 1, 2, ... , k)

1 is a matrix of ones with the same size asI, (.)is the derivative of the heaviside function and

div(.) is the divergence.

2) Optimization of c:

c-= f (b*K)IMi((y))dy (i = 1, ... ,N) (7) f (b2*K)Mg(D(y))dy

3) Optimization of b:

b= IJ(')*K (8) J(2)*K

among them,

N

JCi)- = ciMi(D(y))

N

j(2) ciz2m(D(y)) i=1

<2> Background noise removal (Rol setting)

MRI is vulnerable to noise contamination. If noise in background pixels is severe, it will

influence final curve evolution using multi-phase level set method. While the gray values of MRIs differ greatly in background and organs. Thus, thresholding is a simple but efficient method to remove the background from MRI. In our method, we first extract coarse regions of a target using

Otsu's method. Subsequently, we retrieve connected components from the binary image and select

the contour outlines in the image whose areas are greater than a predefined A value (Step 3

below). Finally, we fill the regions bounded by the preceding selected contours, which serve as the

Rol in this paper. Assuming gray values of background are much lower than the foreground ones,

the procedure can be outlined as follows:

2) Apply a fixed-level threshold to each array element of image f and obtain the binary

image g:

(255, f(i, j) > iT g(i,) D o0, otherwise

among them, constant 0 < i ! 1 can guarantee that the regions of foreground with lower gray

values are not ignored. Based on empirical observations this constant commonly takes values in

the interval [0.7, 0.9].

denoted by B = UBj, is the final ROI image, with the notation R in which the value 1 represents

the points to be further processed:

R1, B(ij) = 255 ) 0, otherwise

<3> Multi-phase level set method with bias field polyfit estimation

simple functions:

b(y) = w T G (y) (11)

is presented as:

G(y) {xjy''-10 ! i O p, 0 j ! i}. (12) Substituting (11) into (3), the energy functional E is converted into:

E(W, c, D) = f Edx (13)

among them, Ex = >N f K(y - x)|I(x) - WTG(y) 2 Mi(D(x))dy.

Therefore, = -2v + 2Aw, where OW

v = f f G(y)I(x)( 1 ciMi(D(x))K(y - x))dydx

Sf I(x) N ciMi(D(x))(f G(y)K(y - x)dy)dx A- ~ (13) A = ff G (y)G T(y)(EN Jc2Mi ((D(x)) K(y - x)) dydx

= f N 2cMi (CD(x)) (f G(y) G T (y) K(y - x) dy) dx

Applying the convolution operation, the equations above will be given by:

v f (G * K)I(x) N ciMi (D(x))dx

A= f (GGT * K) 1c2Mi ((x))dx

Note that v E RLx, A E RLxL. Let OW = 0, we get:

w A-1v (15)

and,

b(y) wT G(y) (16)

localize and remove noise.

<4> Customized initial level set functions to detect abnormalities

especially the negligible patches.

In order to segment more regions of interest, such as bright and dark blobs, we formulate initial level set functions so as to make initial contours scatter linearly, equally spaced in image domain. Classification of anatomical structure appearances in a MRI slice of the brain, suggests that the Rol mentioned in this invention can be classified into light, gray and dark parts; namely, two level set functions constructing 3 membership functions (k = 2, N = 3). The pseudo-code for definition of PO(x, y), P (x, y) is shown below, where the image dimension is W x H, p > 0 is a constant with any value, and the Rol is bounded by a rectangle {(x, y)|xmin x X max, ymin 5

Y 5 ymax):

1) Initialization: #P(x,y) = p, P (x, y) = p, x = 0: 1: W - 1, y = 0: 1: H - 1.

2) Evaluation 1: 00(x, y) = -p, 2 (x + 1, y) = -P, X =Xmin: 2:Xmax, Y = Ymin: 2: Ymax.

3) Evaluation 2: 0 (x, y) = -p , 2(x + 1, y) = -p , X = Xmin: 2 Xmax, : Y = Ymin

+ 1: 2: Ymax.

4) Rol clipped: #1 = #1 - R, 02 = 02 -R.

M 1(x) = 1 - H(q# 1 (x)) M 2 (x) = H(#1 (x))H(# 2 (x)) (17) M 3 (x) = H(#P1 (x))(1 - H(#2(x)))

among them, 0 Mi(x) 1 and 2 1 MZ (x)= 1.

<5> ROI based numerical implementation

In numerical practice, the Dirac function 6(x) in (7) is smoothed as:

S (x) = Ia f2+x2( (18)

Iteration for (jis proposed by:

#l+ = #P + T A (19)

In this paper, the window function K is taken as a Gaussian function with a standard

deviation o, denoted by K,.Our proposed Rol based numerical implementation is presented as

follows:

1) ROI mask. I = I - R.

2) Initialization. Set the following parameters: N, ij, p, u, v, y, E and T. Initialize b and

01,...,#. Due to their invariability during iterative procedures, we prepare G, GGT, 1* K',

(1 * K') . 1 2 , (G*K,)- Iand GGT * Ka.

(18), (6) and (19). Finally, update b via (14), (15) and (16).

4) ROI mask. #j = #j - R for next iteration. ) Redo Step 3, until i exceeds the predefined maximal iteration.

FIGURES 2021100714

(a) (b) Figure 1

(a) (b) Figure 2

Figure 3

(a) (b) (c) Figure 4

(a) (b) (c)

(d) (e) (f)

(g) (h) (i) Figure 5