CN108876769B - Left auricle CT image segmentation method - Google Patents

Abstract

A left auricle CT image segmentation method relates to medical image processing, and comprises the steps of firstly carrying out noise reduction processing on a left auricle CT image, constructing missing boundary information of the left auricle by using priori knowledge, and obtaining a single initial segmentation image as an Atlas single marker by using a region growing algorithm; performing improved Atlas segmentation by using the information relationship between the single mark and the adjacent sequences of the CT of the left auricle to obtain an initial segmentation sequence; and finally, obtaining an accurate segmentation result sequence through a rapid level set algorithm so as to achieve the purpose of segmenting the CT image of the left auricle.

Description

Left auricle CT image segmentation method

Technical Field

The invention relates to medical image processing, in particular to a left atrial appendage CT image segmentation method.

Background

The left atrial appendage is a component of the left atrium of the body and functions to assist in pumping blood. The multi-branch structure of the left auricle is easy to block the blood backflow, when the left auricle can not produce contraction, the atrial pressure is increased, atrial fibrillation is caused, and the structure of the left auricle has important influence on cardiovascular and cerebrovascular diseases. When the left atrial appendage is subjected to an occlusion operation, how to select a proper occlusion device according to the shape of the left atrial appendage is the basis of successful operation. Therefore, the left atrial appendage CT image segmentation is one of the bases for analyzing and processing the left atrial appendage image, and has an important meaning. However, the left atrial appendage CT image segmentation has a certain difficulty due to the characteristics of no boundary between the left atrial appendage and the left atrial adhesion part, no texture feature in the left atrial appendage CT image, and the like. The research on the existing left atrial appendage image segmentation is less, most relevant to the method of the invention is Huang 38891, Gardenia et al (Huang 38891; Gardenia, Korean Louis, Liuqi, white-silk, Denghua, He-Ling.) provide a left atrial appendage ultrasound image segmentation algorithm, which automatically extracts the outline of the ultrasound left atrial appendage image [ J ]. Sichuan university newspaper (engineering science edition), 2016,48(03):87-93.), realizes the segmentation function of the left atrial appendage ultrasound image, realizes the accurate segmentation of the left atrial appendage from a single ultrasound image, and the algorithm firstly carries out two-step position identification and then uses an active outline model to obtain accurate segmentation through the initial outline close to the real edge, the position determination of the first two steps needs to rely on priori knowledge, and the algorithm has low robustness. In addition, the method carries out left atrial appendage image segmentation processing on the ultrasonic image, and can not be directly used for left atrial appendage CT image segmentation.

The left atrial appendage has anatomical features of variable shape and complex structure. The left atrial appendage CT image has the following characteristics: the left auricle and the left atrium adhesive part have no boundary, the left auricle has no texture characteristics substantially, and the like, so the boundary of the left auricle CT image can not be directly obtained, and the effect of using the traditional image segmentation algorithm is poor.

Disclosure of Invention

The present invention is directed to provide a left atrial appendage CT image segmentation method for overcoming the above-mentioned shortcomings of the prior art.

The method comprises the steps of firstly, carrying out noise reduction on a CT image of the left auricle, constructing missing boundary information of the left auricle by using priori knowledge, and obtaining a single initial segmentation image as an Atlas single marking image by using a region growing algorithm; performing improved Atlas segmentation by using the information relationship between the single mark and the adjacent sequences of the CT of the left auricle to obtain an initial segmentation sequence; and finally, obtaining an accurate segmentation result sequence through a rapid level set algorithm so as to achieve the purpose of segmenting the CT image of the left auricle.

The invention comprises the following steps:

1) inputting a left auricle CT sequence image A;

2) carrying out noise reduction processing on the sequence image A in the step 1) to obtain a sequence image B;

3) selecting a single CT image B in the sequence image B in the step 2)_iThe remaining sequence images are denoted as B';

4) according to the single CT image B obtained in the step 3)_iManually selecting two characteristic points according to the position characteristic information of the middle left auricle, constructing the boundary information of the left auricle missing, and obtaining an image C;

5) obtaining an image L by using a region growing algorithm for the image C obtained in the step 4)_i；

6) With the image L obtained in step 5)_iFor the initial marker image, Atlas segmentation of the image B' in adjacent registration is performed to obtain an initial segmentation sequenceL；

7) Calculating approximate Euclidean distances from all foreground points to the background points of the sequence image B in the step 2) to obtain a distance image sequence D;

8) and taking the contour of the initial segmentation sequence L obtained in the step 6) as an initial contour, taking the distance image sequence D obtained in the step 7) as a characteristic image, and performing accurate segmentation on the left atrial appendage CT sequence image on the sequence image B by using a fast level set algorithm.

In step 2), the sequence image a is subjected to noise reduction processing, and a specific method for obtaining the sequence image B may be:

(1) the image of the ith image in the sequence image A after being subjected to the iterative filtering for t times is marked as I_t,iRecalculated using the curvature anisotropy diffusion equation:

wherein c (| x |) represents the conductance,

v represents taking the gradient, t represents the number of iterations, and I (x, y, I) represents the ith image in the sequence image a;

(2) iterating the process of the step (1) for n times to obtain an ith image subjected to noise reduction;

(3) increasing i by an increment, and enabling i to be i +1, executing the step (1), and obtaining an image sequence B after preserving edges and reducing noise;

in the step 3), selecting a single CT image B in the sequence images B in the step 2)_iThe specific method of (3) may be:

(1) calculating the number n of images in the CT sequence image A of the left auricle;

(2) selecting the ith image as a picture B_iWherein i is the largest integer no greater than n/2.

In step 4), the single CT image B obtained in step 3) is used_iThe specific method for acquiring the image C by manually selecting two feature points and constructing the boundary information of the left auricle missing from the position feature information of the middle left auricle comprises the following steps:

(1) let two manually selected seed points be p₁(x₁，y₁)、p₂(x₂，y₂) These two points can determine a straight line in a logical coordinate system

In the physical coordinate system, the coordinate system is,

can pass through from point p₁To point p₂Fitting and forming pixel points extending one by one, namely point p₁The next pixel point position of (a) is expressed as:

wherein N is₀Indicating the position of the next pixel, N₁，N₂Representing candidate positions, d₁，d₂Respectively represent N₁，N₂Distance to a point above and below the logical coordinate;

(2) the straight line can be known according to the position relation on the left auricle anatomical structure

The slope k of (a) satisfies: k e (-45 °, 0 °) U (-90 °, -45 °), then its corresponding N₁，N₂Expressed as:

(3) when the slope k of the straight line satisfies k ∈ (-45 °, 0 °):

wherein p is_i+1Denotes the selection criterion of the next location, dx denotes the logical horizontal offset value, dy denotes the logical vertical offset value, in terms of p_i+1The value of (a) can determine the result of the next coordinate:

when straight line

When the slope k of (a) satisfies k e (-90 °, -45 °):

wherein p is_i+1Indicating the selection criteria for the next location, dx indicating the logical horizontal offset value and dy indicating the logical vertical offset value. According to p_i+1The value of (a) can determine the result of the next coordinate:

(4) modifying image B_iMiddle coordinate N₀Gray value of the pixel:

I(N₀)＝0

an image C is obtained.

In step 5), obtaining an image L by using a region growing algorithm for the image C obtained in step 4)_iThe specific method of (3) may be:

(1) selecting a pixel point in the image C as a seed point s according to the position of the left auricle in the image C, setting a threshold, and adding the point s into a seed point set seed when the input seed point meets the following threshold:

I₀(s)∈[lower,upper]

wherein, I₀(s) represents the gray value of the seed point s in the image C, lower represents the lower limit of the gray value, and upper represents the upper limit of the gray value;

(2) for each s_iE is determined, calculate s_iAll elements U in the eight connected domain_iBelongs to U, all the points meeting the threshold value condition set in the step (1) are added into a seed point set seed,repeating the step (2) until no new point is added into the seed point set seed to obtain an image L_i。

In step 6), with the image L obtained from step 5)_iFor the initial marker image, Atlas segmentation of adjacent registration is performed on the image sequence B', and a specific method for obtaining the initial segmentation sequence L may be:

(1) the serial number of the image Bi in the sequence image B is marked as i, and the next image in the sequence image B is marked as B_i+1A 1 to B_iAnd the image L obtained in the step 5)_iIs marked as a mark pair (B)_i，L_i) A 1 to B_iAs floating images, B_i+1Image registration for fixed images, obtaining transformation parameters f_i+1. Image registration was achieved using an elastix tool (s.klein, m.staring, k.murphy, m.a.viargever, j.p.w.pluim, "elastix: atomic for intensive based imaging," IEEE transaction on medical imaging, vol.29, No.1, pp.196-205, January2010) with the following specific steps:

(a) image B using a pyramid model_i，B_i+1Dividing into three precision grades;

(b) image B_i，B_i+1Mapping the physical coordinates into logical coordinates by an interpolation method;

(c) from picture B separately by means of random sampling_i，B_i+1Selecting a plurality of sample points as a sample point set S representing the whole_i，S_i+1；

(d)S_i，S_i+1The similarity between the two is expressed by mutual information;

(e) the optimization algorithm of the registration uses an adaptive gradient descent algorithm;

(f) using a multilevel registration strategy, three transformation modes of rigid transformation, affine transformation and B spline transformation (D.Rueckert, L.I.Sonoda, C.Hayes, D.L.G.Hill, M.O.Leach, and D.J.Hawks.Nonrigig-transforming from-to-form are used in sequence to obtain a transformation parameter f._i+1。

(2) Converting the parameter f_i+1Acting on the marker image L_iObtaining B_i+1Is marked with an image L_i+1；

(3) B is to be_i+1And L_i+1As a new label pair (B)_i+1，L_i+1) Repeating the steps (1) and (2) until the end;

(4) image B_iAs a floating image, with the last image B of the sequence image B_i-1Image registration according to the method described in step (1), wherein B_i-1For fixing the image, transformation parameters f are obtained_i-1；

(5) Converting the parameter f_i-1Acting on the marker image L_iObtaining B_i-1Is marked with an image L_i-1；

(6) B is to be_i-1And L_i-1As a new label pair (B)_i-1，L_i-1) And (5) repeating the steps (4) and (5) until the end.

In step 7), the step of calculating the approximate euclidean distances from all foreground points to the background points for the sequence image B in step 2) may be performed by: (see Akmal B M, Maragos P. optimal design of a chain distance transformations [ J ]. IEEE Transactions on imaging Processing A Publication of the IEEE Signal Processing Society,1998,7(10): 1477-84.):

(1) each image in the sequence image B is processed as follows: sliding the distance transformation template from left to right and from top to bottom, sequentially calculating the minimum Euclidean distance from each pixel to the background, and obtaining an intermediate distance image sequence B_d；

(2) For intermediate range image sequence B_dEach image of (a) is processed as follows: and sliding the distance transformation template from right to left and from bottom to top, and sequentially calculating the minimum Euclidean distance from each pixel to the background to obtain a distance characteristic image sequence D.

In step 8), the contour of the initial segmentation sequence L obtained in step 6) is used as an initial contour, the range image sequence D obtained in step 7) is used as a feature image, and the sequence image B is subjected to accurate segmentation of the left atrial appendage CT sequence image by using a fast level set algorithm, which may be specifically performed by (see document: Anon.Insight segmentation and integration Toolkit [ EB/OL ] [2008-07-06]. http:// www.itk.org.):

(1) calculating each image Di (i is 1,2,3.. times.n) in the distance characteristic image sequence D by using a sigmoid function to form an input characteristic DS of the rapid level set algorithm_i(i ═ 1,2,3.. n), where the sigmoid function is:

wherein I represents an image D_iMax and Min are images D_iMiddle maximum and minimum gray value, alpha, represents image D_iThe width of the luminance range (window width), β the luminance at the center of the range (window level), I' the output image, i.e. DS_i(ii) a All DS_i(i 1,2,3.. n) is a fast level set algorithm feature sequence and is marked as DS;

(2) taking the sequence L obtained in step 6) as an initial contour sequence of a fast level set algorithm, taking the DS obtained in step 8) as a feature sequence of the fast level set algorithm, and for each image Li, i being 1,2,3.. n in the sequence L, the level set algorithm sets Li as a zero level set of a high-dimensional function, which is called a level set function: y (X, t), then the level set function moves into a differential equation, the profile of the motion is obtained by extracting the zero level set G ((X), t) { Y (X, t) } 0} from the output, and the update of the differential equation solution ψ is calculated using a level set equation:

wherein, a (x) represents a horizontal convection coefficient, p (x) represents a propagation coefficient (also called an expansion coefficient), z (x) represents a spatial modifier coefficient of a curvature mean, and α, β, and γ represent weights of respective parameters, and are constants.

According to the method, the left auricle image is segmented according to the characteristics of the left auricle image with missing boundary information, and the image sequence is subjected to noise reduction processing to reduce the influence caused by noise; then, in order to obtain an initial contour model sequence, a single initial segmentation image containing contour model information is obtained by supplementing missing boundary information to a single image; then, the single initial segmentation image and the corresponding gray level image form an initial Atlas mark pair, and a segmentation sequence is obtained through an Atlas algorithm which is firstly registered and then segmented; the Atlas mark sequence is an initial segmentation sequence containing an initial contour model, and contour information is input into a level set algorithm in sequence to obtain a corresponding segmentation result sequence. The invention finally achieves the purpose of segmenting the left auricle of the image, namely separating the left auricle with partial missing boundaries from the left ventricle image connected with the left auricle.

Constructing boundary deletion of image deletion of a sticky part in the initial segmentation process of a single left atrial appendage through priori knowledge, and obtaining single initial segmentation by using a segmentation algorithm; using an improved Atlas algorithm to convert a single initial segmentation image into an initial segmentation sequence image, using similarity of adjacent images in the sequence image to meet the requirement of the Atlas algorithm on a gray scale image in a mark pair, using the obtained single initial segmentation gray scale image to register with a gray scale image adjacent to the single initial segmentation gray scale image in the sequence, applying a transformation result of the registration to the Atlas mark (initial segmentation image), obtaining Atlas marks of the adjacent gray scale images, combining the Atlas marks into a new Atlas mark pair, and iterating the steps to form a complete Atlas mark sequence; the accurate segmentation uses a distance map to replace a gray map so as to obtain distinguishing information of the left auricle and the left atrium, and then the initial contour information is input into a rapid level set algorithm, so that the purpose of accurate segmentation is finally achieved.

Compared with the prior art, the invention has the beneficial effects that:

1. aiming at the characteristics of the CT image of the left auricle, the invention utilizes the priori knowledge to supplement the information of the boundary missing part, makes up the overgrowth defect of the region growing algorithm connected with the threshold, and successfully obtains the first initial segmentation containing the contour information.

2. The invention uses an improved Atlas algorithm, combines a single initial segmentation image and a corresponding gray map into an initial mark pair, then uses the gray map in the initial mark pair to register with the gray map adjacent to the sequence position, acts the registered transformation result on the Atlas mark (initial segmentation image) to obtain the Atlas marks of the adjacent gray images, combines the Atlas mark pair into a new Atlas mark pair, and iterates the steps to form a complete Atlas mark sequence, thereby achieving the purpose of automatically obtaining the initial segmentation sequence containing the outline information.

3. The improved Atlas algorithm used by the invention utilizes the correlation between sequences, solves the problem of difficult acquisition of similar label pairs, not only improves the segmentation precision of the traditional single Atlas algorithm, but also greatly reduces the time complexity of registration and fusion in multiple Atlas algorithms.

4. According to the method, the distance image is used for replacing the gray level image to increase the position information, so that the boundary identification capability of a level set algorithm in the weak boundary segmentation problem is improved; and the contour map obtained by the improved Atlas algorithm is combined to improve the segmentation precision, so that the precise segmentation of the left atrial appendage sequence is realized.

5. The method increases the missing boundary information in the CT sequence image of the left auricle and improves the segmentation precision by three steps, and the methods supplement each other and are buckled with each other in a ring mode, and the left auricle image can be segmented more accurately by combining the methods, so that the method has the advantages of high segmentation precision and strong robustness.

The method aims at the characteristics of the left auricle image with missing boundary information, and carries out processing and analysis in multiple steps, so as to achieve the purpose of accurately segmenting the left auricle image from the CT image. The steps of single mark acquisition, Atlas acquisition of an initial contour sequence, distance feature acquisition and rapid level set segmentation can supplement each other in function, the single mark acquisition mainly utilizes prior information to supplement missing boundary information by combining with prior knowledge, and then utilizes a region growing algorithm connected with a threshold to acquire segmentation, the algorithm is simple and the time complexity is low; the improved Atlas algorithm obtains an initial contour sequence, provides an initial contour which is as close to a real boundary as possible for a level set algorithm, reduces the time complexity of the level set algorithm, prompts missing boundary position information and improves the segmentation accuracy; the distance characteristics can solve the problem of insufficient boundary information in the gray level image, improve the accuracy of the contour of a level set algorithm at a position without a boundary, and improve the robustness; the fast level set algorithm can actively generate a closed contour and generate an accurate segmentation result sequence image with the help of initial contour information and distance characteristics. Compared with the closest prior art, the method has obvious distinguishing characteristics and has obvious advantages: the method can generate the sequence segmentation result, has strong practicability, does not need excessive manual intervention, can generate the accurate left atrial appendage sequence segmentation result by manually selecting three points, and is superior in accuracy and robustness.

Drawings

FIG. 1 is a basic flow diagram of the present invention;

FIG. 2 is a left atrial appendage CT image sequence according to an embodiment of the present invention;

FIG. 3 is a distance transform template according to an embodiment of the present invention;

fig. 4 is one of the final segmentation result image sequences according to the embodiment of the present invention.

Detailed Description

The invention is described in further detail below with reference to the figures and examples, which are explained in detail in the flow chart shown in fig. 1.

The invention comprises the following steps:

1) inputting a left auricle CT sequence image A;

2) carrying out noise reduction processing on the sequence image A in the step 1) to obtain a sequence image B;

3) selecting a single CT image B in the sequence image B in the step 2)_iThe remaining sequence images are denoted as B';

4) according to the image B obtained in the step 3)_iManually selecting two characteristic points according to the position characteristic information of the middle left auricle, constructing the boundary information of the left auricle missing, and obtaining an image C;

5) obtaining an image L by using a region growing algorithm for the image C obtained in the step 4)_i；

6) With the image L obtained in step 5)_iFor the initial marker image, Atlas segmentation of image B' in adjacent registration is performed to obtain an initial segmentation sequence L.

7) Calculating approximate Euclidean distances from all foreground points to the background points of the sequence image B in the step 2) to obtain a distance image sequence D;

8) and taking the contour of the initial segmentation sequence L obtained in the step 6) as an initial contour, taking the distance image sequence D obtained in the step 7) as a characteristic image, and performing accurate segmentation on the left atrial appendage CT sequence image on the sequence image B by using a fast level set algorithm.

One of the left atrial appendage CT image sequences of an embodiment of the present invention is shown in fig. 2.

In step 2), the sequence image a is subjected to noise reduction processing, and a specific method for obtaining the sequence image B may be:

(1) the image of the ith image in the sequence image A after being subjected to the iterative filtering for t times is marked as I_t,iRecalculated using the curvature anisotropy diffusion equation:

wherein c (| x |) represents the conductance, wherein,

v represents taking the gradient, t represents the number of iterations, and I (x, y, I) represents the ith image in the sequence image a;

in this embodiment, the values of the parameters are as follows: c is 5 and t is 4.

(2) Iterating the process of the step (1) for n times to obtain an ith image subjected to noise reduction;

(3) increasing i by an increment, and enabling i to be i +1, executing the step (1), and obtaining an image sequence B after preserving edges and reducing noise;

in step 3), selecting a single CT image B in the sequence images B_iThe specific method comprises the following steps:

(1) calculating the number n of images in the CT sequence image A of the left auricle;

(2) selecting the ith image as a picture B_iWherein i is the largest integer no greater than n/2;

in this embodiment, the values of the parameters are as follows: the 30 th image is selected when n is 60 and i is 30.

In step 4), two feature points are manually selected from the image Bi obtained in step 3), and boundary information of the left atrial appendage missing is constructed, wherein the specific method for obtaining the image C comprises the following steps:

(1) let two manually selected seed points be p₁(x₁，y₁)、p₂(x₂，y₂) These two points can determine a straight line in a logical coordinate system

In the physical coordinate system, the coordinate system is,

can pass through from point p₁To point p₂Fitting and forming pixel points extending one by one, namely point p₁The next pixel point position of (a) is expressed as:

wherein N is₀Indicating the position of the next pixel, N₁，N₂Representing candidate positions, d₁，d₂Respectively represent N₁，N₂Distance to a point above and below the logical coordinate;

(2) the straight line can be known according to the position relation on the left auricle anatomical structure

The slope k of (a) satisfies: k e (-45 °, 0 °) U (-90 °, -45 °), then its corresponding N₁，N₂Expressed as:

(3) when the slope k of the straight line satisfies k ∈ (-45 °, 0 °):

wherein p is_i+1Indicating the selection criteria for the next location, dx indicating the logical horizontal offset value and dy indicating the logical vertical offset value. According to p_i+1The value of (a) can determine the result of the next coordinate:

when straight line

When the slope k of (a) satisfies k e (-90 °, -45 °):

wherein p is_i+1Indicating the selection criteria for the next location, dx indicating the logical horizontal offset value and dy indicating the logical vertical offset value. According to p_i+1The value of (a) can determine the result of the next coordinate:

(4) modifying image B_iMiddle coordinate N₀Gray value of the pixel:

I(N₀)＝0

obtaining an image C;

in this embodiment, the values of the parameters are as follows:

p1(x1，x2)＝(115,66)，p2(x2，y2)＝(166，110)

in step 5), obtaining an image L by using a region growing algorithm for the image C obtained in step 4)_iThe specific method comprises the following steps:

(1) selecting a pixel point in the image C as a seed point s according to the position of the left auricle in the image C, setting a threshold, and adding the point s into a seed point set seed when the input seed point meets the following threshold:

I₀(s)∈[lower,upper]

wherein, I₀(s) represents the coordinate value of the seed point s in the image C, lower represents the lower limit of the gray value, and upper represents the upper limit of the gray value.

(2) For each s_iE is determined, calculate s_iAll elements U in the eight connected domain_iE.g. U, adding the seed point set seed to the points meeting the set threshold condition in the step (1), repeating the step (2) until no new points are added to the seed point set seed, and obtaining an image L_i；

In this embodiment, the values of the parameters are as follows: lower 137, upper 255, and s (152, 85).

In step 6), with the image L obtained from step 5)_iFor the initial marker image, Atlas segmentation of adjacent registration is performed on the image sequence B', and a specific method for obtaining an initial segmentation sequence L is as follows:

(1) the serial number of the image Bi in the sequence image B is marked as i, and the next image in the sequence image B is marked as B_i+1. B is to be_iAnd the image L obtained in the step 5)_iIs marked as a mark pair (B)_i，L_i) A 1 to B_iAs floating images, B_i+1Image registration for fixed images, obtaining transformation parameters f_i+1. Image registration was achieved using an elastix tool (s.klein, m.staring, k.murphy, m.a.viargever, j.p.w.pluim, "elastix: atomic for intensive based imaging," IEEE transaction on medical imaging, vol.29, No.1, pp.196-205, January2010) with the following specific steps:

(a) image B using a pyramid model_i，B_i+1Dividing into three precision grades;

(b) image B_i，B_i+1Mapping the physical coordinates into logical coordinates by an interpolation method;

(c) from respective images B by means of random sampling_i，B_i+1Selecting a plurality of sample points as a sample point set S representing the whole_i，S_i+1；

(d)S_i，S_i+1The similarity between them is expressed by mutual information(see, P.Th' evenazandM.Unser. optimization of multiresolution imaging registration. IEEETranss. image processing, 9(12): 2083-2099, 2000.);

mutual information can be expressed as:

wherein L is_FAnd L_MIs a collection of regularly spaced intensity distribution centers, p is a discrete joint probability, and p is_FAnd p_MThe edge dispersion probabilities of the stationary and moving images are obtained by summing m and f, respectively. The joint probability is estimated using a B-spline Parzen window:

wherein, ω is_FAnd ω_MRepresenting a fixed and moving B-spline Parzen window. Scale constant σ_FAnd σ_MIs equal to L_FAnd L_MA defined intensity bin width. These are directly according to I_FAnd I_MAnd the number | L of histogram bins specified by the user_FL and L_M|；

(e) The optimization algorithm for registration uses an adaptive gradient descent algorithm (references: s.klein, j.p.w.pluim, m.staring, and m.a.viewer.adaptive storage gradient specific optimization for image registration. international Journal of computer Vision,81(3): 227-;

(f) the transformation parameter f is obtained by using a multilevel registration strategy and sequentially using three transformation modes of rigid transformation, affine transformation and B spline transformation (reference documents: D.Rueckert, L.I.Sonoda, C.Hayes, D.L.G.Hill, M.O.Leach, and D.J.Hawks.Nonrigig transformation from and for formation: Application to break MR images. IEEETrans.Med.Imag, 18(8): 712-721, 1999)_i+1。

The position between the floating image and the fixed image is first adjusted using a rigid transformation:

wherein R (x) represents a rotation matrix, C represents a rotation center, t represents an original transformation parameter set, μ represents a parameter vector, and is composed of an Euler angle (rad) and a translation vector, i.e., a vector μ with a length of 3, wherein θ_zRepresenting rotation about an axis perpendicular to the image, t_x、t_yRespectively representing translations in x and y directions;

and then, using affine transformation to adjust the overall shape:

wherein C represents the rotation center, t represents the original transformation parameter group, mu represents the parameter vector, A (x) represents the matrix capable of realizing translation, rotation, scaling and shearing, and the parameter vector mu is composed of matrix element a_ijAnd a translation vector t_x、t_yComposition t_x、t_yRespectively representing translations in x and y directions;

finally, local deformation is carried out by using B spline transformation:

wherein x is_kDenotes a control point, β³(x) Represents a cubic multi-dimensional B-spline polynomial (M.Unser. Splines: Aperfect for Signal and image processing. IEEE Signal processing. Mag.,16(6): 22-38,1999.), p_kRepresents the B-spline coefficient vector (control point displacement), sigma represents the B-spline control point spacing, N_xRepresenting a set of points constrained by control of the x points in the B-spline; the control point grid is defined by the amount of space between control points (σ ═ g [ ([ sigma ])₁，...，σ_d) (where d represents the image size), which may be different for each direction; the number of control points P ═ (P1., Pd) the number of parameters M is determined by M ═ P1 ×. Pd) × d; the parameter vectors are as follows: mu ═ p (p)_1x，p_2x，...，pP₁，p_1y，p_2y，...，pP₂)^T(ii) a Thereby obtaining a transformation parameter f_i+1。

(2) Converting the parameter f_i+1Acting on the marker image L_iObtaining B_i+1Is marked with an image L_i+1；

(3) B is to be_i+1And L_i+1As a new label pair (B)_i+1，L_i+1) Repeating the steps (1) and (2) until the end;

(4) image B_iAs a floating image, with the last image B of the sequence image B_i-1Image registration according to the method described in step (1), wherein B_i-1For fixing the image, transformation parameters f are obtained_i-1；

(5) Converting the parameter f_i-1Acting on the marker image L_iObtaining B_i-1Is marked with an image L_i-1；

(6) B is to be_i-1And L_i-1As a new label pair (B)_i-1，L_i-1) Repeating the steps (4) and (5) until the end;

in this embodiment, the values of the parameters are as follows: 2048, 3, Iter (Iter)_r，Iter_a，Iter_B)＝(800，1200，1000)，|L_F|＝|L_M|＝(16，32，64)

In step 7), the specific method for calculating the approximate euclidean distances from all foreground points to the background points for the sequence image B in step 2) to obtain the distance image sequence D includes: (see Akmal B M, Maragos P. optimal design of a chain distance transformations [ J ]. IEEE Transactions on imaging Processing A Publication of the IEEE Signal Processing Society,1998,7(10): 1477-84.):

(1) each image in the sequence image B is processed as follows: sliding the distance transformation template from left to right and from top to bottom, sequentially calculating the minimum Euclidean distance from each pixel to the background, and obtaining an intermediate distance image sequence B_d；

(2) For intermediate range image sequence B_dEach image of (a) is processed as follows: sliding the distance transformation template from right to left and from bottom to top, and sequentially calculating the minimum Euclidean distance from each pixel to the background to obtain a distance characteristic image sequence D;

in this embodiment, the distance transformation template is shown in fig. 3, where values of each parameter are as follows: a is 5, b is 7, c is 11;

in step 8), the contour of the initial segmentation sequence L obtained in step 6) is used as an initial contour, the range image sequence D obtained in step 7) is used as a feature image, and a fast level set algorithm is used to perform accurate segmentation of the left atrial appendage CT sequence image on the sequence image B, which is specifically described in the following steps (see document: Anon.Insight segmentation and integration Toolkit [ EB/OL ] [2008-07-06]. http:// www.itk.org.):

(1) calculating each image Di (i is 1,2,3.. times.n) in the distance characteristic image sequence D by using a sigmoid function to form an input characteristic DS of the rapid level set algorithm_i(i ═ 1,2,3.. n), where the sigmoid function is:

wherein I represents an image D_iMax and Min are images D_iMiddle maximum and minimum gray value, alpha, represents image D_iThe width of the luminance range (window width), β the luminance at the center of the range (window level), I' the output image, i.e. DS_i(ii) a All DS_i(i 1,2,3.. n) is a fast level set algorithm feature sequence and is marked as DS;

(2) taking the sequence L obtained in step 6) as an initial contour sequence of a fast level set algorithm, taking the DS obtained in step 8) as a feature sequence of the fast level set algorithm, and for each image Li (i ═ 1,2,3.... n) in the sequence L, setting Li to a zero level set of a high-dimensional function, where the high-dimensional function is called a level set function: y (X, t), and then the level set function motion becomes a differential equation, and the profile of the motion is obtained by extracting the zero level set G ((X), t) = { Y (X, t) ═ 0} from the output. The update of the differential equation solution ψ is calculated using a level set equation:

wherein, a (x) represents a horizontal convection coefficient, p (x) represents a propagation coefficient (also called an expansion coefficient), z (x) represents a spatial regulator coefficient of a curvature mean, and α, β, and γ represent weights of each parameter and are constants;

in this embodiment, the values of the parameters are as follows: p (x) the fast level set algorithm signature sequence DS obtained using the procedure of (1) in step 8). Alpha, beta, gamma are set as default values in the ITK development kit (see literature: Anon. Insightsegmentation and registration Toolkit [ EB/OL ] [2008-07-06]. http:// www.itk.org.). Finally, a result image of the left atrial appendage CT image segmentation is obtained, as shown in fig. 4.

Claims (8)

1. A left atrial appendage CT image segmentation method is characterized by comprising the following steps:

1) inputting a left auricle CT sequence image A;

2) carrying out noise reduction processing on the sequence image A in the step 1) to obtain a sequence image B;

3) selecting a single CT image B in the sequence image B in the step 2)_iThe remaining sequence images are denoted as B';

4) according to the single CT image B obtained in the step 3)_iManually selecting two characteristic points according to the position characteristic information of the middle left auricle, constructing the boundary information of the left auricle missing, and obtaining an image C;

5) obtaining an image L by using a region growing algorithm for the image C obtained in the step 4)_i；

6) With the image L obtained in step 5)_iPerforming Atlas segmentation of adjacent registration on the image B' for an initial marker image to obtain an initial segmentation sequence L;

7) calculating approximate Euclidean distances from all foreground points to the background points of the sequence image B in the step 2) to obtain a distance image sequence D;

8) and taking the contour of the initial segmentation sequence L obtained in the step 6) as an initial contour, taking the distance image sequence D obtained in the step 7) as a characteristic image, and performing accurate segmentation on the left atrial appendage CT sequence image on the sequence image B by using a fast level set algorithm.

2. The method for segmenting the CT image of the left atrial appendage according to claim 1, wherein in the step 2), the sequential image a is subjected to noise reduction processing, and the specific method for obtaining the sequential image B is as follows:

(1) the image of the ith image in the sequence image A after being subjected to the iterative filtering for t times is marked as I_t,iRecalculated using the curvature anisotropy diffusion equation:

wherein c (| x |) represents the conductance,

representing gradient taking, t representing iteration number, and I (x, y, I) representing ith image in the sequence image A;

(2) iterating the process of the step (1) for n times to obtain an ith image subjected to noise reduction;

(3) and (4) incrementing i, and making i equal to i +1, executing the step (1) to obtain the image sequence B after the edge-preserving noise reduction.

3. The method as claimed in claim 1, wherein in step 3), the single CT image B in the sequence image B in step 2) is selected_iThe specific method comprises the following steps:

(1) calculating the number n of images in the CT sequence image A of the left auricle;

(2) selecting the ith image as a picture B_iWherein i is the largest integer no greater than n/2.

4. The method for segmenting the CT image of the left atrial appendage of claim 1, wherein in the step 4), the single CT image B obtained in the step 3) is used as the basis_iThe specific method for acquiring the image C comprises the following steps of manually selecting two characteristic points from the position characteristic information of the middle left auricle, constructing the boundary information of the left auricle missing, and acquiring the image C:

(1) let two manually selected seed points be p₁(x₁，y₁)、p₂(x₂，y₂) These two points define a straight line l in a logical coordinate system, in the physical coordinate system, l passing from the point p₁To point p₂Fitting and forming pixel points extending one by one, namely point p₁The next pixel point position of (a) is expressed as:

wherein N is₀Indicating the position of the next pixel, N₁，N₂Representing candidate positions, d₁，d₂Respectively represent N₁，N₂Distance to a point above and below the logical coordinate;

(2) according to the position relation on the left auricle anatomical structure, the slope k of the straight line l satisfies the following conditions: k e (-45 °, 0 °) U (-90 °, -45 °), then its corresponding N₁，N₂Expressed as:

(3) when the slope k of the straight line satisfies k ∈ (-45 °, 0 °):

wherein p is_i+1Denotes the selection criterion of the next location, dx denotes the logical horizontal offset value, dy denotes the logical vertical offset value, in terms of p_i+1The value of (a) determines the result of the next coordinate:

when the slope k of the straight line l satisfies k ∈ (-90 °, -45 °):

wherein p is_i+1Indicating the selection criteria of the next location, dx indicating the logical horizontal offset value and dy indicating the logical vertical offset value; according to p_i+1The value of (a) determines the result of the next coordinate:

(4) modifying image B_iMiddle coordinate N₀Gray value of the pixel:

I(N₀)＝0

an image C is obtained.

5. The method for segmenting the CT image of the left atrial appendage of claim 1, wherein in the step 5), the image L is obtained by using a region growing algorithm for the image C obtained in the step 4)_iThe specific method comprises the following steps:

(1) selecting a pixel point in the image C as a seed point s according to the position of the left auricle in the image C, setting a threshold, and adding the point s into a seed point set seed when the input seed point meets the following threshold:

I₀(s)∈[lower,upper]

wherein, I₀(s) represents the gray value of the seed point s in the image C, lower represents the lower limit of the gray value, and upper represents the upper limit of the gray value;

(2) for each s_iE is determined, calculate s_iAll elements U in the eight connected domain_iE.g. U, adding seed point set seed to the points meeting the set threshold condition in the step (1), and repeating the step (2) until no new points are added to the seed pointsSeed is collected to obtain an image L_i。

6. The method as claimed in claim 1, wherein in step 6), the image L obtained from step 5) is used for segmentation_iFor the initial marker image, Atlas segmentation of adjacent registration is performed on the image sequence B', and a specific method for obtaining an initial segmentation sequence L is as follows:

(1) the serial number of the image Bi in the sequence image B is marked as i, and the next image in the sequence image B is marked as B_i+1A 1 to B_iAnd the image L obtained in the step 5)_iIs marked as a mark pair (B)_i，L_i) A 1 to B_iAs floating images, B_i+1Image registration for fixed images, obtaining transformation parameters f_i+1(ii) a The image registration is realized by using an elastix tool, and the specific steps are as follows:

(a) image B using a pyramid model_i，B_i+1Dividing into three precision grades;

(b) image B_i，B_i+1Mapping the physical coordinates into logical coordinates by an interpolation method;

(c) from picture B separately by means of random sampling_i，B_i+1Selecting a plurality of sample points as a sample point set S representing the whole_i，S_i+1；

(d)S_i，S_i+1The similarity between the two is expressed by mutual information;

(e) the optimization algorithm of the registration uses an adaptive gradient descent algorithm;

(f) using a multilevel registration strategy, and sequentially using three transformation modes of rigid transformation, affine transformation and B spline transformation to obtain a transformation parameter f_i+1；

(2) Converting the parameter f_i+1Acting on the marker image L_iObtaining B_i+1Is marked with an image L_i+1；

(3) B is to be_i+1And L_i+1As a new label pair (B)_i+1，L_i+1) Repeating the steps (1) and (2) until the end;

(4) image B_iAs a floating image, with the last image B of the sequence image B_i-1Image registration according to the method described in step (1), wherein B_i-1For fixing the image, transformation parameters f are obtained_i-1；

(5) Converting the parameter f_i-1Acting on the marker image L_iObtaining B_i-1Is marked with an image L_i-1；

(6) B is to be_i-1And L_i-1As a new label pair (B)_i-1，L_i-1) And (5) repeating the steps (4) and (5) until the end.

7. The method as claimed in claim 1, wherein in step 7), the approximate euclidean distances between all foreground points and background points are calculated for the sequence image B in step 2), and the specific method for obtaining the distance image sequence D is as follows:

(1) each image in the sequence image B is processed as follows: sliding the distance transformation template from left to right and from top to bottom, sequentially calculating the minimum Euclidean distance from each pixel to the background, and obtaining an intermediate distance image sequence B_d；

(2) For intermediate range image sequence B_dEach image of (a) is processed as follows: and sliding the distance transformation template from right to left and from bottom to top, and sequentially calculating the minimum Euclidean distance from each pixel to the background to obtain a distance characteristic image sequence D.

8. The left atrial appendage CT image segmentation method according to claim 1, wherein in the step 8), the contour of the initial segmentation sequence L obtained in the step 6) is used as an initial contour, the range image sequence D obtained in the step 7) is used as a feature image, and a fast level set algorithm is used to perform accurate segmentation on the left atrial appendage CT sequence image on the sequence image B, and the specific method is as follows:

(1) calculating each image Di, i in the distance characteristic image sequence D as 1,2,3.. n, wherein n is the total number of the distance characteristic image sequence D by using a sigmoid function, and forming input characteristics of a fast level set algorithmDS_i1,2,3.. n, wherein the sigmoid function is:

wherein I represents an image D_iMax and Min are images D_iMiddle maximum and minimum gray value, alpha, represents image D_iThe width of the luminance range, β the luminance in the center of the range, I' the output image, i.e. DS_i(ii) a All DS_iA fast level set algorithm signature sequence is marked as DS, wherein i is 1,2,3.. n;

(2) taking the sequence L obtained in step 6) as an initial contour sequence of a fast level set algorithm, taking the DS obtained in step 8) as a feature sequence of the fast level set algorithm, and for each image Li, i being 1,2,3.. n in the sequence L, the level set algorithm sets Li as a zero level set of a high-dimensional function, which is called a level set function: y (X, t), then the level set function moves into a differential equation, the profile of the motion is obtained by extracting the zero level set G ((X), t) { Y (X, t) } 0} from the output, and the update of the differential equation solution ψ is calculated using a level set equation:

wherein, a (x) represents a horizontal convection coefficient, p (x) represents a propagation coefficient, z (x) represents a spatial modifier coefficient of a curvature mean, and α, β, and γ represent weights of respective parameters, which are constants.

Images