System and method for automated detection of mucus plugs within bronchial tree in MSCT images
View Patent ↗A method for detecting and localizing mucus plugs in digitized lung images, includes providing a digitized lung image volume comprising a plurality of intensities corresponding to a 3-dimensional grid of points, extracting a bronchial tree from said lung image, said bronchial tree comprising a plurality of branching airways terminating at terminal points, providing a model of a 2-dimensional cross section of an airway, selecting an extended point beyond a terminal point of an airway branch in a direction of said airway branch, obtaining a 2-dimensional cross section I of size m×n points from said lung image about said selected point, processing said 2D cross section I by calculating a local neighborhood function for each point in the cross section and forming a union of all local neighborhood functions, and calculating a correlation between processed 2D cross section and said airway model, wherein said correlation is indicative of the presence of a mucus plug within said airway.
1. A method for detecting and localizing mucus plugs in digitized lung images comprising the steps of:
providing a digitized lung image volume comprising a plurality of intensities corresponding to a 3-dimensional grid of points;
extracting a bronchial tree from said lung image, said bronchial tree comprising a plurality of branching airways terminating at terminal points;
providing a model of a 2-dimensional cross section of an airway;
selecting an extended point beyond a terminal point of an airway branch in a direction of said airway branch;
obtaining a 2-dimensional cross section I of size m×n points from said lung image about said selected point;
processing said 2D cross section I by calculating a local neighborhood function for each point in the cross section and forming a union of all local neighborhood functions; and
calculating a correlation between processed 2D cross section and said airway model, wherein said correlation is indicative of the presence of a mucus plug within said airway.
2. The method of claim 1 , wherein extracting a bronchial tree from said lung image comprises segmenting said image volume using an adaptive threshold method and applying an iterative region-growing until terminal points of the airway branches are reached, and skeletonizing to obtain a bronchial tree structure comprising a set of branches, each branch comprising a set of ordered sites.
3. The method of claim 1 , further comprising computing a distance map along said bronchial tree to obtain a lumen diameter along the bronchial tree.
4. The method of claim 3 , wherein said airway model is a 2-dimensional p×p image template defined as
K
(
x
,
y
)
=
(
k
0
,
0
⋯
k
0
,
p
-
1
⋮
⋰
⋮
k
p
-
1
,
0
⋯
k
p
-
1
,
p
-
1
)
,
wherein x, yεZ, xε[0, . . . , p−1], and yε[0, . . . p−1], and wherein
k
x
,
y
=
{
Φ
2
if
h
<
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
≤
p
/
2
,
Φ
1
if
h
≥
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
,
Φ
0
otherwise
,
wherein Φ 0 , Φ 1 , Φ 2 are distinct integers with Φ 0 <Φ 1 <Φ 2 , h is an airway lumen diameter, and (p/2)−h is an airway wall thickness.
5. The method of claim 4 , wherein Φ 0 =0, Φ 1 =1, and Φ 2 =2.
6. The method of claim 4 , wherein calculating the lumen diameter of said airway model for the airway branch of the extended point comprises selecting a portion of said airway branch, selecting a plurality of cross-sections along said branch portion, obtaining a lumen diameter for each cross section from said distance map, and calculating an average of each cross section lumen diameter.
7. The method of claim 6 , wherein calculating the wall thickness of said airway model for the airway branch comprises selecting a cross-section at a mid point of said branch portion, projecting a plurality of rays from a skeleton mid point of said cross section, finding a wall thickness from an intensity profile along each ray using a full-width-half-max criteria, and averaging the wall thickness of each ray.
8. The method of claim 6 , wherein the wall thickness is estimated as 1/(lumen diameter).
9. The method of claim 7 , further comprising estimating lower and upper intensity thresholds for healthy airways by selecting a cross section at the mid point of said branch, setting said lower intensity threshold t 1 to the maximum intensity of the lumen in said cross section, estimating the upper intensity threshold t 2 by averaging intensity values along said plurality of rays within an estimated wall thickness interval, and estimating a mucus minimum density d muc as a maximum intensity value over said branch portion.
10. The method of claim 7 , further comprising estimating lower and upper intensity thresholds for healthy airways by setting said thresholds to pre-determined values, and estimating a mucus minimum density as being approximately equal to said lower threshold.
11. The method of claim 4 , wherein calculating a local neighborhood function for each point in the cross section comprises defining a local p-neighborhood W ij of each point (x, y) of said cross-section as
W
ij
(
x
,
y
)
=
(
w
0
,
0
ij
⋯
w
0
,
p
-
1
ij
⋮
⋰
⋮
w
p
-
1
,
0
ij
⋯
w
p
-
1
,
p
-
1
ij
)
,
∀
i
∈
[
0
,
…
,
m
-
1
]
,
∀
j
∈
[
0
,
…
,
n
-
1
]
,
wherein said union of all local neighborhood functions is
⋃
i
,
j
W
ij
for iε[0, . . . , m−1] and jε[0, . . . , n−1], wherein, ∀xε[0, . . . , p−1], ∀yε[0, . . . , p−1], for a clean airway,
w
x
,
y
ij
=
{
Φ
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
2
,
Φ
1
if
t
1
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
1
,
Φ
0
otherwise
,
wherein t 1 and t 2 are respectively lower and upper intensity bounds for a healthy lumen, and for a mucus-filled airway,
w
x
,
y
ij
=
{
Φ
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
2
,
Φ
1
if
d
muc
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
1
,
Φ
0
otherwise
,
wherein d muc is a mucus minimum density.
12. The method of claim 11 , wherein said correlation C between processed 2D cross section and said airway model is calculated as
C
ij
=
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
w
x
,
y
ij
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
k
x
,
y
wherein 0≦C ij ≦1, and wherein a correlation greater a predefined first tolerance value is indicative of a mucus plug and a correlation greater a predefined second tolerance value less than said first tolerance value is indicative of a clean airway.
13. A method for detecting and localizing mucus plugs in digitized lung images comprising the steps of:
providing a digitized lung image volume comprising a plurality of intensities corresponding to a 3-dimensional grid of points;
correlating a p×p airway model to an m×n 2D airway cross section I extracted from said lung image from
C
ij
=
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
w
x
,
y
ij
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
k
x
,
y
,
wherein 0≦C ij ≦1, x, yεZ, xε[0, . . . , p−1] and yε[0, . . . , p−1],
k
x
,
y
=
{
2
if
h
<
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
≤
p
/
2
,
1
if
h
≥
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
,
0
otherwise
,
wherein h is a lumen diameter of said airway, and (p/2)−h is a wall thickness of said airway, wherein for a clean airway,
w
x
,
y
ij
{
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
2
,
1
if
t
1
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
1
,
0
otherwise
,
wherein t 1 and t 2 are respectively lower and upper intensity bounds for a healthy lumen, wherein for a mucus-filled airway,
w
x
,
y
ij
{
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
2
,
1
if
d
muc
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
1
,
0
otherwise
,
wherein d muc is a mucus minimum density, and wherein a correlation greater a predefined high tolerance value is indicative of a mucus plug and a correlation greater a predefined low tolerance value less than said first tolerance value is indicative of a clean airway.
14. The method of claim 13 , further comprising extracting a bronchial tree from said lung image, said bronchial tree comprising a plurality of branching airways terminating at terminal points; and computing a distance map along said bronchial tree to obtain a lumen diameter along the bronchial tree.
15. The method of claim 14 , wherein extracting said m×n 2D airway cross section I comprises selecting an extended point beyond a terminal point of an airway branch in a direction of said airway branch; extracting said 2-dimensional cross section I of size m×n points from said lung image about said extended point.
16. The method of claim 13 , further comprising forming a union of all local neighborhood functions
⋃
W
ij
i
,
j
for iε[0, . . . , m−1] and jε[0, . . . , n−1] wherein
W
ij
(
x
,
y
)
=
(
w
0
,
0
ij
…
w
0
,
p
-
1
ij
⋮
⋰
⋮
w
p
-
1
,
0
ij
…
w
p
-
1
,
p
-
1
ij
)
,
∀
i
∈
[
0
,
…
,
m
-
1
]
,
∀
j
∈
[
0
,
…
,
n
-
1
]
.
17. A non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executable by the computer to perform the method steps for detecting and localizing mucus plugs in digitized lung images, said method comprising the steps of:
providing a digitized lung image volume comprising a plurality of intensities corresponding to a 3-dimensional grid of points;
extracting a bronchial tree from said lung image, said bronchial tree comprising a plurality of branching airways terminating at terminal points;
providing a model of a 2-dimensional cross section of an airway;
selecting an extended point beyond a terminal point of an airway branch in a direction of said airway branch;
obtaining a 2-dimensional cross section I of size m×n points from said lung image about said selected point;
processing said 2D cross section I by calculating a local neighborhood function for each point in the cross section and forming a union of all local neighborhood functions; and
calculating a correlation between processed 2D cross section and said airway model, wherein said correlation is indicative of the presence of a mucus plug within said airway.
18. The computer readable program storage device of claim 17 , wherein extracting a bronchial tree from said lung image comprises segmenting said image volume using an adaptive threshold method and applying an iterative region-growing until terminal points of the airway branches are reached, and skeletonizing to obtain a bronchial tree structure comprising a set of branches, each branch comprising a set of ordered sites.
19. The computer readable program storage device of claim 17 , the method further comprising computing a distance map along said bronchial tree to obtain a lumen diameter along the bronchial tree.
20. The computer readable program storage device of claim 19 , wherein said airway model is a 2-dimensional p×p image template defined as
K
(
x
,
y
)
=
(
k
0
,
0
…
k
0
,
p
-
1
⋮
⋰
⋮
k
p
-
1
,
0
…
k
p
-
1
,
p
-
1
)
,
wherein x, yεZ, xε[0, . . . , p−1], and yε[0, . . . , p−1], and wherein
k
x
,
y
=
{
Φ
2
if
h
<
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
≤
p
/
2
,
Φ
1
if
h
≥
(
k
x
,
y
,
k
p
/
2
,
p
/
2
)
,
Φ
0
otherwise
,
wherein Φ 0 , Φ 1 , Φ 2 are distinct integers with Φ 0 <Φ 1 <Φ 2 , h is an airway lumen diameter, and (p/2)−h is an airway wall thickness.
21. The computer readable program storage device of claim 20 , wherein Φ 0 =0, Φ 1 =1, and Φ 2 =2.
22. The computer readable program storage device of claim 20 , wherein calculating the lumen diameter of said airway model for the airway branch of the extended point comprises selecting a portion of said airway branch, selecting a plurality of cross-sections along said branch portion, obtaining a lumen diameter for each cross section from said distance map, and calculating an average of each cross section lumen diameter.
23. The computer readable program storage device of claim 22 , wherein calculating the wall thickness of said airway model for the airway branch comprises selecting a cross-section at a mid point of said branch portion, projecting a plurality of rays from a skeleton mid point of said cross section, finding a wall thickness from an intensity profile along each ray using a full-width-half-max criteria, and averaging the wall thickness of each ray.
24. The computer readable program storage device of claim 22 , wherein the wall thickness is estimated as 1/(lumen diameter).
25. The computer readable program storage device of claim 23 , the method further comprising estimating lower and upper intensity thresholds for healthy airways by selecting a cross section at the mid point of said branch, setting said lower intensity threshold t 1 to the maximum intensity of the lumen in said cross section, estimating the upper intensity threshold t 2 by averaging intensity values along said plurality of rays within an estimated wall thickness interval, and estimating a mucus Minimum density d muc as a maximum intensity value over said branch portion.
26. The computer readable program storage device of claim 23 , the method further comprising estimating lower and upper intensity thresholds for healthy airways by setting said thresholds to pre-determined values, and estimating a mucus minimum density as being approximately equal to said lower threshold.
27. The computer readable program storage device of claim 20 , wherein calculating a local neighborhood function for each point in the cross section comprises defining a local p-neighborhood W ij of each point (x, y) of said cross-section as
W
ij
(
x
,
y
)
=
(
w
0
,
0
ij
…
w
0
,
p
-
1
ij
⋮
⋰
⋮
w
p
-
1
,
0
ij
…
w
p
-
1
,
p
-
1
ij
)
,
∀
i
∈
[
0
,
…
,
m
-
1
]
,
∀
j
∈
[
0
,
…
,
n
-
1
]
,
wherein said union of all local neighborhood functions is
⋃
i
,
j
W
ij
for iε[0, . . . , m−1] and jε[0, . . . , n−1], wherein, ∀xε[0, . . . , p−1], ∀yε[0, . . . , p−1], for a clean airway,
w
x
,
y
ij
=
{
Φ
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
2
,
Φ
1
if
t
1
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
1
,
Φ
0
otherwise
,
wherein t 1 and t 2 are respectively lower and upper intensity bounds for a healthy lumen, and for a mucus-filled airway,
w
x
,
y
ij
=
{
Φ
2
if
t
2
≤
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
2
,
Φ
1
if
d
muc
≥
I
(
i
-
(
p
-
1
2
)
+
x
,
j
-
(
p
-
1
2
)
+
y
)
and
k
i
,
j
≥
Φ
1
,
Φ
0
otherwise
,
wherein d muc is a mucus minimum density.
28. The computer readable program storage device of claim 27 , wherein said correlation C between processed 2D cross section and said airway model is calculated as
C
ij
=
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
w
x
,
y
ij
∑
x
=
0
p
-
1
∑
y
=
0
p
-
1
k
x
,
y
wherein 0≦C ij <1, and wherein a correlation greater a predefined first tolerance value is indicative of a mucus plug and a correlation greater a predefined second tolerance value less than said first tolerance value is indicative of a clean airway.