Flattening of the Lung Surface with Temporal Consistency for the
Follow-Up Assessment of Pleural Mesothelioma
Peter Faltin
1
, Thomas Kraus
2
, Marcin Kopaczka
1
and Dorit Merhof
1
1
Institute of Imaging and Computer Vision, RWTH Aachen University, Aachen, Germany
2
Institute and Out-patient Clinic of Occupational Medicine, Uniklinik RWTH Aachen, Aachen, Germany
Keywords:
Surface Representation, Planar Visualization, Multidimensional Scaling, Flattening, Temporal Consistency.
Abstract:
Malignant pleural mesothelioma is an aggressive tumor of the lung surrounding membrane. The standardized
workflow for the assessment comprises an inspection of 3D CT images to detect pleural thickenings which
act as indicators for this tumor. Up to now, the visualization of relevant information from the pleura has only
been superficially addressed. Current approaches still utilize a slice-wise visualization which does not allow
a global assessment of the lung surface. In this publication, we present an approach which enables a planar
2D visualization of the pleura by flattening its surface. A distortion free mapping to a planar representation
is generally not possible. The present method determines a planar representation with low distortions directly
from a voxel-based surface. For a meaningful follow-up assessment, the consistent representation of a lung
from different points in time is highly important. Therefore, the main focus in this publication is to guarantee
a consistent representation of the pleura from the same patient extracted from images taken at two different
points in time. This temporal consistency is achieved by our newly proposed link of both surfaces during the
flattening process. Additionally, a new initialization method which utilizes a flattened lung prototype speeds
up the flattening process.
1 INTRODUCTION
Asbestos is a silicate mineral which was commonly
used in many countries because of its advantageous
properties, e.g. resistance to acid and heat. During
mechanical processing, asbestos fibers can get into
the lung by inhalation. A deposition of these fibers in
the tissue can cause different kinds of cancer such as
the malignant pleural mesothelioma (MPM) which is
an aggressive tumor of the pleura. The pleura is a lung
surrounding membrane consisting of two layers. Be-
cause of the long latency of about 38 years and a peak
usage of asbestos in the 1980s, an increasing number
of MPM cases is expected in Europe till 2018 (Pis-
tolesi and Rusthoven, 2004). In some other countries,
especially in Asia, the asbestos consumption is still
rising (Virta, 2012). Therefore, an increasing number
of MPM cases is expected in these countries. Because
of the aggressive tumor growth, early stage detection
is absolutely essential. For this purpose, high risk
patients undergo a regular medical check-up. These
checks also include CT image acquisition to identify
pleural thickenings which are indicators for the MPM.
Especially the identification of their growth is an im-
Figure 1: CT slice with thickenings indicated by orange cir-
cles.
portant aspect in the diagnosis. The growth can be es-
timated by a follow-up assessment, based on images
from two different points in time.
As shown in Figure 1, pleural thickenings are
difficult to spot because of the low image con-
trast between their tissue and the surrounding tis-
sue. At an early stage, they have a small thick-
Faltin, P., Kraus, T., Kopaczka, M. and Merhof, D.
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma.
DOI: 10.5220/0005718702330243
In Proceedings of the 11th Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2016) - Volume 2: IVAPP, pages 235-245
ISBN: 978-989-758-175-5
Copyright
c
2016 by SCITEPRESS Science and Technology Publications, Lda. All rights reserved
235
ness in terms of the image resolution of typically
1-3 voxels. Therefore, a manual analysis of the im-
ages is time consuming and subject to strong inter-
and intra-reader variability (Ochsmann et al., 2010).
Several (semi-)automated methods (Rudrapatna et al.,
2005; Chen et al., 2011; Frauenfelder et al., 2011;
Sensakovic et al., 2011; Faltin et al., 2013; Kengne
Nzegne et al., 2013) have been developed to sup-
port the physician in the decision process. Never-
theless, the final decision is made by the physician.
For this purpose, a visualization of the relevant data
is highly important. Two publications (Faltin et al.,
2013; Kengne Nzegne et al., 2013) propose a global
3D visualization of the thickening locations, while
other publications only consider 2D slice-wise visual-
ization. With the 3D visualization a global overview
is possible. However, due to the 3D shaped surface,
thickenings might be occluded and the affected areas
on the surface can only be qualitatively recognized. A
slice-wise view facilitates a more precise assessment
of the affected surface, but it is time consuming, does
not visualize the complex thickening morphologies,
and does not give an overview. To address these is-
sues, a planar visualization of the lung surface is sug-
gested in this publication. An exemplary mapping of a
3D lung surface into a planar visualization is shown in
Figure 2. This kind of visualization enables a global
overview and a precise recognition of the affected ar-
eas on the surface. Different features, e.g. tissue at-
tenuation or surface roughness are utilized in a man-
ual or automated thickening identification. These fea-
tures can be visualized for the whole lung at once by
utilizing the planar representation. Beside the detec-
tion of thickenings, it is important to compare the dif-
ferent points in time for a growth assessment. For this
purpose the planar representations of the lung at both
points in time need to be consistent regarding their
spatial arrangement.
The concept of mapping a 2D surface from a 3D
space to a planar 2D space is often referred to as
flattening. To provide a meaningful interpretation of
the flattened surface, the mapping should preserve an-
gles and be equiareal, as far as possible. If both an-
gle and area are preserved the mapping is called iso-
metric or distance-preserving. However, an isomet-
ric mapping of arbitrary shaped surfaces is not gen-
erally possible (Floater and Hormann, 2005). The
aim of flattening methods is to approximate this iso-
metric mapping. Most of these methods (Zigelman
et al., 2002; Bronstein et al., 2003; Tosun et al., 2004;
Zhao et al., 2013; Hurdal and Stephenson, 2009) are
based on triangular meshes. In our case, the rele-
vant feature information is available for each of the
350 000 lung surface voxels. For a mesh-based flat-
(a) 3D lung surface. (b) 2D planar visualization of
lung surface from Figure 2a.
Figure 2: Exemplary mapping of the 3D surface to a planar
2D visualization. The 3D coordinates of each surface point
are taken as color components for a RGB coding.
tening a very highly resolved mesh needs to be pro-
cessed or resampled. This results either in high com-
putational demands or the need of additional mesh
simplification methods which might induce inaccu-
racies. An efficient alternative is the method from
Saroul et al. (Saroul et al., 2006) which can be di-
rectly applied on discrete voxel data and focuses on an
interactive visualization with adjustable planes and a
center of low distortion. However, it does not include
a global optimization of the chosen plane. There-
fore, our approach utilizes multidimensional scaling
(MDS) which has been applied for the flattening of
discrete surfaces (Grossmann et al., 2002) and pro-
vides a global optimization on a discrete voxel sur-
face. With this approach, the preserved surface dis-
tances are directly calculated on the high resolution
data set and the computation can still be performed on
a down-sampled set of surface points. In this publica-
tion, we suggest a new method to link both points in
time globally to achieve consistent representations of
the surfaces. Existing approaches to compare similar
shaped surface (Tosun et al., 2004; Zhao et al., 2013;
Hurdal and Stephenson, 2009) often use a standard-
ized mapping target. The approach from Bronstein
et al. (Bronstein et al., 2003) was especially designed
for an inter- and intra-subject comparison. However,
none of these approaches uses a mapping target which
is individually optimized for each subject. In con-
trast to these methods, our method results in subject
adapted representations and facilitates the comparison
between different points in time. Additionally, we de-
duce a method to efficiently initialize the MDS pro-
cess to speed up the processing of a similar shaped
surface. For this purpose, the results of a flattening
from a prototype lung are utilized. A closed-form so-
lution to directly compute the initial points from this
prototype is derived.
IVAPP 2016 - International Conference on Information Visualization Theory and Applications
236
adjacency matrix
dissimilarity matrices
2D points after MDS
considered interpolated planar representation
MDS
t = t
1
:
BFS
NMS
Tr
Tr
p
t = p:
BFS
t = t
2
:
BFS
lung surface I
I
I
of set I
I
I
of set I
I
I
sub,r
of set I
I
I
sub,r
of set I
I
I
Figure 3: Overview of the presented method: First, an adjacency matrix is extracted from the lung surface of the two points in
time t
1
and t
2
and the prototype surface t = p. By a non-maxima suppression (NMS) the set I
I
I of all surface points is reduced to
the set I
I
I
sub,r
I
I
I. From the full adjacency matrices and the reduced point set, the dissimilarities are extracted by a breadth first
search (BFS). Then, the multidimensional scaling (MDS) extracts a subsampled 2D representation. Finally, an interpolation
by radial basis functions (RBF) is used to determine a full planar visualization. The prototype t = p can be used to initialize
the MDS. To connect the different lung surfaces, the transformations Tr and Tr
p
are required.
2 METHODS
In the first part of this section, the extraction of the
discrete lung surface is shortly described. This is fol-
lowed by a general introduction to MDS. The idea
of mapping a discrete 3D surface to a 2D plane by
MDS was introduced by Grossmann et al. (Gross-
mann et al., 2002) and is described in Section 2.3, 2.4
and 2.6. We embed this method into a follow-up as-
sessment which requires a temporal consistent flatten-
ing at both points in time. This methodical extension
is described in Section 2.5. To speed up the flattening
processing, we also derive an initialization approach
to map the data from a flattened prototype lung to an
actual lung, described in Section 2.7. An overview of
the complete work flow is shown in Figure 3.
2.1 Discrete Lung Surface
The extraction of the discrete lung surface R
R
R from
the chest CT image data g is performed in multiple
steps, as shown in Figure 4. Firstly, the body is seg-
mented by a simple thresholding, keeping the largest
connected component, followed by a morphological
filling. The chosen threshold of -250 HU is not critical
regarding the exact value. Within this body segment,
the lungs are identified. In this case the exact thresh-
old of -250 HU is relevant and conforms to the lower
boundary of the standard thickening window. This
window is used to visualize the CT data g in a stan-
dardized manual assessment of pleural thickenings.
Additionally, the mediastinal region R
R
R
medias
in be-
tween left and right lung is identified. It is the largest
concave region of the lung surfaces in each slice. Fi-
nally, a discrete surface which guarantees connectiv-
ity within a 26-connected neighborhood N
26
is ex-
tracted. For this purpose, a lung mask with removed
6-neighborhood N
6
is subtracted from the segmented
lung R
R
R itself
R
R
R
full
= R
R
R \ (R
R
R N
6
), (1)
with the morphological erosion . The surface part
surrounded from the mediastinal region R
R
R
medias
is ex-
cluded, because it is not relevant for the thickening
assessment. The resulting lung surface is given by
R
R
R = R
R
R
full
\ (R
R
R
medias
N
26
), (2)
with the morphological dilation .
2.2 Multidimensional Scaling
MDS is an approach to map a set of N points into
a space of any dimension M, by considering only
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma
237
Figure 4: 2D visualization of the lung surface segmentation for an exemplary slice. First, the body (green) is detected, then
lungs (yellow and orange) and mediastinum (blue) are segmented. Finally, the pleura surface is extracted (yellow and orange).
the dissimilarities δ
i, j
of the data points, with i, j
I
I
I =
{
1. . . N
}
. This section gives a short idea of the
method. Details can be found in the relevant litera-
ture (Cox and Cox, 2001). The results are only mean-
ingful if the dissimilarities can be, at least approxi-
mately, described in the target space of dimension M.
In the target space, the positions p
p
p
MD,i
are estimated
depending on their distances d
i, j
=
p
p
p
MD,i
p
p
p
MD, j
.
These distances should be as close as possible to the
associated dissimilarities in the original space. This is
achieved by solving
R
R
R
MD
= arg min
p
p
p
MD,1
...p
p
p
MD,N
loss
(p
p
p
MD,1
, . . . , p
p
p
MD,N
),
(3)
with a certain loss measure and the matrix
contain-
ing the dissimilarities δ
i, j
. The MDS methods in lit-
erature are differentiated into
classical vs. non-classical and
metric vs. non-metric
approaches. These terms and the resulting loss crite-
ria named stress or strain are sometimes ambiguously
used. In this publication the terminology from Buja
et al. (Buja et al., 2008) is chosen and non-classical
MDS with the loss criterion
loss
(p
p
p
MD,1
, . . . , p
p
p
MD,N
)
=Sstress
(p
p
p
MD,1
, . . . , p
p
p
MD,N
)
=
i, j=1...N
i6= j
d
2
i, j
δ
2
i, j
2
i, j=1...N
i6= j
δ
4
i, j
(4)
is utilized. This criterion is differentiable even if two
points are coincident. The solution is found by a
gradient-based minimization, utilizing the conjugate-
gradient method (Navon and Legler, 1987). Weight-
ing terms w
i, j
for each pair of points can be chosen
individually, resulting in the new criterion
Sstress
,W
W
W
(p
p
p
MD,1
, . . . , p
p
p
MD,N
)
=
i, j=1...N
i6= j
w
i, j
d
2
i, j
δ
2
i, j
2
i, j=1...N
i6= j
w
i, j
δ
4
i, j
, (5)
with the weights denoted by the matrix W
W
W of same
size as
. In the case of a planar visualization, the
target dimension is chosen as M = 2.
2.3 Distances on Discrete Surface
As described in the last section, the optimization for
the 2D representation aims for distances d
i, j
identical
to the dissimilarities δ
i, j
. For a surface flattening, the
dissimilarities are chosen as the distances on the lung
surface between all surface points p
p
p
3D
R
R
R (Gross-
mann et al., 2002). To calculate the distances between
the discrete points on the surface, first a weighted ad-
jacency matrix with distances between the connected
surface points p
p
p
3D,i
with i I
I
I, I
I
I =
{
1. . .
|
R
R
R
|}
of a
26-connected neighborhood is extracted. Based on
this information, a breadth first search (BFS), based
on the Boost C++ Library (Gleich, 2009), is used to
determine the distances on the surface. A lung sur-
face has approximately
|
R
R
R
|
= 350 000 surface points.
Just to store all pair distances, approximately 450 GB
of memory would be necessary. So, the distances
are only extracted for a subsampled number of sur-
face points. A regular subsampling of the unparam-
eterized 3D surface is not easily possible. There-
fore, non-maxima suppression (Neubeck and Gool,
2006) (NMS) with a suppression radius r is applied to
the binary surface R
R
R. The free parameter r represents
a compromise between exact surface description and
computational complexity. In this case, the NMS on a
binary image is only utilized to generate a uniformly
distributed subset of points R
R
R
sub,r
R
R
R, with
p
p
p
3D,i
p
p
p
3D, j
r, i, j I
I
I, i 6= j. (6)
The indices of the subsampled surface points are con-
tained in the set I
I
I
sub,r
I
I
I. Also other point distribu-
IVAPP 2016 - International Conference on Information Visualization Theory and Applications
238
tion algorithms could be applied. The NMS was cho-
sen to be conforming to the procedure in a subsequent
step.
2.4 Surface Mapping for a Single Point
in Time
The sample points describe a 2D surface embed-
ded in a 3D space, therefore a meaningful 2D rep-
resentation can be obtained, but of course a perfect
isometric mapping cannot be achieved. All sample
points are mapped depending on their dissimilari-
ties δ
i, j
= δ(p
p
p
3D,i
, p
p
p
0
3D, j
), i, j I
I
I
sub,r
and indepen-
dent on their actual 3D position p
p
p
3D,i
. These dis-
similarity measures δ
i, j
are given by the surface dis-
tances from Section 2.3. Solving the minimization
problem in Equation 3 leads to the 2D representa-
tion R
R
R
2D
. This representation is not unique regard-
ing the degrees of freedom (DOF) rotation, transla-
tion and reflection, as these properties are not observ-
able from the dissimilarities. For a single surface all
weights w
i, j
are kept equal to one, if all regions are
of similar importance. Focusing on some important
regions which are, e.g. subject to thickenings can be
achieved by increasing the weights in these regions.
2.5 Consistent Mapping for Two Points
in Time
The process described in the previous section, could
be independently applied to the lungs’ surfaces of
the same patient from images, taken at two different
points in time. In most cases, this leads, beside the
DOFs, to a similar spatial arrangement of the points
in 2D which is however not guaranteed. The problem
described in Equation 3 might have multiple solutions
and also local minima. To avoid different solutions at
both points in time, we propose to join both indepen-
dent problems into a single problem. For this pur-
pose, identical points on the lung surfaces R
R
R
t
1
and
R
R
R
t
2
from both points in time t
1
and t
2
are linked by
enforcing identical target positions in the 2D visual-
ization.
Firstly, suitable points to establish a link are se-
lected. They should be located at the same 3D sur-
face position for both points in time, after a registra-
tion Tr : p
p
p
3D,t
1
,i
7→ p
p
p
3D,t
2
,i
0
. This registration is usu-
ally not perfect and afterwards the points Tr
p
p
p
3D,t
1
,i
and p
p
p
3D,t
2
,i
0
are not identical. For each surface
point p
p
p
3D,t
1
,i
, i I
I
I the closest point at the other point
in time and their temporal distance is extracted by uti-
lizing the k-d tree algorithm (Bentley, 1975).
conn(i) = arg min
i
0
=1...
|
R
R
R
t
2
|
Tr
p
p
p
3D,t
1
,i
p
p
p
3D,t
2
,i
0
2
(7)
d
t
(i) = p
p
p
3D,t
1
,i
p
p
p
3D,t
2
,conn(i)
, (8)
In the subsampling process, point pairs which can be
well aligned by the registration should be locally pre-
ferred to pairs which can only be poorly aligned by the
registration. This quality can be measured by the tem-
poral distance d
t
(i). Therefore, the NMS (Neubeck
and Gool, 2006) with radius of r is now applied to
the inverse distances
1
d
t
(i)
to obtain the subsampled
set I
I
I
sub,r
for the linking process. The flattening in re-
gions of high importance could be improved by keep-
ing additional points in these regions. Points within
growing thickening regions will result in inconsistent
surface distances at both points in time and should be
avoided. Strong growth also results in large temporal
distances. Hence, thickening points which are critical
for the flattening are suppressed automatically.
In the second step, the minimization problems are
independently formulated as in the previous section
and then connected. The set of linked points does not
have to be identical with the set of subsampled points
in Section 2.3. However, the problem statement is
given by the matrices W
W
W and
, whose size depends
on the number of involved points. Therefore, the point
sets are chosen to be identical to reduce the number
of variables in the optimization problem. The matri-
ces
t
1
and
t
2
are the dissimilarities for the different
points with identical order of the points at both points
in time. So, the entries δ
δ
δ
t
1
,i, j
, δ
δ
δ
t
2
,conn(i),conn( j)
corre-
spond to the same point combination at both points
in time, because point i corresponds to point i
0
and j
corresponds to j
0
. Then, the matrices describing the
linked problem are
conn
=
t
1
0
0
0
0
0
0
t
2
, (9)
W
W
W
conn
=
1
1
1 κ diag(1, ..., 1)
κ diag(1, ..., 1) 1
1
1
,(10)
with 1
1
1 being matrices of size
|
I
I
I
sub,r
|
×
|
I
I
I
sub,r
|
with
ones in each entry, diag(1, ..., 1) being matrices of
same size with ones on the diagonal and zeros at each
other entry, and 0
0
0 being a matrix of same size with
zeros at each entry. The parameter κ is a free weight-
ing parameter to balance the influence of inter- and
intra-temporal distances given in
conn
. Only the di-
agonal entries of the matrices 0
0
0 are relevant to set the
temporal distances of connected points. The other en-
tries are of no importance and are ignored because of
their zero weights in W
W
W
conn
. Finally, the connected
problem can be solved as described in Equation 3.
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma
239
2.6 Interpolation
Only the subset of points I
I
I
sub,r
is mapped to
points p
p
p
2D,c
, c I
I
I
sub,r
on the 2D plane by MDS. For
all other discrete image points p
p
p
3D,i
, i (I
I
I \ I
I
I
sub,r
)
their matching 2D positions p
p
p
2D,i
are estimated
by an interpolation based on radial basis functions
(RBF) (Morse et al., 2005). RBF approximate a func-
tion f by a sum of radially symmetric functions Φ,
around the subsampled points p
p
p
3D,c
with c I
I
I
sub,r
by
f
ζ
(p
p
p
3D,i
) =poly
ζ
(p
p
p
3D,i
)+
cI
I
I
sub,r
λ
ζ
(c)Φ
p
p
p
3D
p
p
p
3D,c
, (11)
for any point p
p
p
3D,i
, where λ
ζ
(c) are individ-
ual weighting terms for each center point p
p
p
3D,c
,
poly
ζ
(p
p
p
3D,i
) represents the polynomial part of the
function, and ζ allows the separate interpolation for
both planar dimensions ζ {x, y}. A resulting planar
point p
p
p
2D,i
is given by p
p
p
2D,i
= ( f
x
(p
p
p
3D,i
), f
y
(p
p
p
3D,i
))
T
.
The weights λ
ζ
(c) and coefficients for the polyno-
mial part poly
ζ
are determined by a least squares ap-
proach (Morse et al., 2005), considering all mapped
points I
I
I
sub,r
. This method has been successfully ap-
plied to various medical interpolation problems (Carr
et al., 1997; Morse et al., 2005). Different RBFs with
similar interpolation results have been applied in pre-
liminary experiments. In this publication, only radi-
ally symmetric function Φ(ρ) = ρ is regarded in detail
which showed a slightly superior performance.
2.7 Initialization
MDS is an optimization process, thus its calculation
time depends on the initialization with the 2D point
positions. To reduce the calculation time a prototype
lung can be utilized. For this prototype, a planar rep-
resentation must be calculated once. Then, it can be
utilized for all future flattening processes. A registra-
tion Tr
p
: p
p
p
3D,i
7→ p
p
p
3D,p,i
between points p
p
p
3D,p,i
R
R
R
p
of the prototype lung surface and the correspond-
ing points p
p
p
i
R
R
R of the actual lung surface R
R
R
is required. For the subsampled points with in-
dices i I
I
I
sub,r
, the initial 2D positions can be approx-
imated from the prototype positions p
p
p
2D,p,i
.
Different physical dimensions of the actual and
prototype lung in 3D cannot be directly transferred to
the planar 2D representation. The 2D offset between
two points is a mixture of the 3D offset, obtained from
the optimization process. However, we deduce an ap-
proximation for a deformation of the flattened proto-
type which considers different 3D dimensions of the
lung. As a result of Section 2.3, the surface distance
between two points p
p
p
i
and p
p
p
j
is calculated by a sum
of distances between neighboring surface points, on
the shortest path E
E
E between those points. This path is
described by the edges (µ, ν)
T
E
E
E and results in the
surface distance
δ
i, j
=
(µ,ν)
T
E
E
E
p
p
p
3D,µ
p
p
p
3D,ν
. (12)
An isometric mapping to 2D is not possible for lung
surfaces and therefore, the distances between the
points p
p
p
2D,µ
and p
p
p
2D,ν
of each edge (µ, ν)
T
are sub-
ject to an error caused by the mapping. The planar
distances can therefore be described by
d
i, j
=
(µ,ν)
T
E
E
E
e
µ,ν
p
p
p
3D,µ
p
p
p
3D,ν
| {z }
δ
µ,ν
, (13)
with the mapping error e
µ,ν
. The prototype and ac-
tual lungs have similar shapes and the aim is a sim-
ilar 2D mapping. Therefore, it is assumed that the
distance d
p,i, j
, of the corresponding points p
p
p
2D,p,i
and
p
p
p
2D,p, j
on the prototype surface, is subject to the same
mapping errors e
µ,ν
and given by the same shortest
path E
E
E. Then, the prototype distance of two neighbor
points is given by
d
p,µ,ν
= e
µ,ν
Tr
p
(p
p
p
3D,µ
) Tr
p
(p
p
p
3D,ν
)
| {z }
δ
0
µ,ν
. (14)
As a result of Equation 13 and 14, the planar distances
of the actual lung can be approximated by
d
i, j
=
(µ,ν)
T
E
E
E
δ
µ,ν
δ
0
µ,ν
d
p,µ,ν
, (15)
with the planar distances d
p,µ,ν
known from the proto-
type and the shortest path E
E
E, which can be efficiently
calculated during the BFS (Gleich, 2009), discussed
in Section 2.3.
By choosing a fixed point p
p
p
2D, f
, with f I
I
I
sub
, the
distances d
f ,i
from Equation 15 to this point can be
used for an initialization. The position p
p
p
2D,i
of all
other points i I
I
I
sub
\ { f } can be described in polar
coordinates, as shown in Figure 5, by
p
p
p
2D,i
= p
p
p
2D, f
+ d
f ,i
cos(φ
f ,i
)
sin(φ
f ,i
)
, (16)
and for the prototype by
p
p
p
2D,p,i
= p
p
p
2D,p, f
+ d
p, f ,i
cos(φ
p, f ,i
)
sin(φ
p, f ,i
)
, (17)
with the unknown angles φ
f ,i
and φ
p, f ,i
. Without loss
of generality these angles and also the fixed point can
be chosen identically
p
p
p
2D, f
= p
p
p
2D,p, f
results in identical translation, and
IVAPP 2016 - International Conference on Information Visualization Theory and Applications
240
(a) Ground truth. Red
indicates thickenings, yel-
low uncertain findings and
green the healthy pleura
surface.
0
1
(b) Predicted thickening con-
fidence in an exemplary clas-
sification.
0 mm
41 mm
(c) Distance to the closest
bone. The ribs are located in
between the line-shaped re-
gions with high distances.
-760 HU
1300 HU
(d) Average attenuation val-
ues of tissue surrounding the
lung. The ribs are located in
the line-shaped regions with
high attenuation.
Figure 6: Planar visualization of lung surface with ground truth of thickening locations, automatic prediction, and exemplary
features considered for the prediction.
d
f ,i
φ
f ,i
p
p
p
2D, f
p
p
p
2D,i
Figure 5: Sketch of the planar lung surface, with description
of the point p
p
p
2D,i
, relative to the fixed point p
p
p
2D, f
in polar
coordinates d
f ,i
and φ
p, f ,i
.
φ
f ,i
= φ
p, f ,i
results in identical rotation and reflec-
tion
of prototype and actual lung. Finally, the initial point
position is approximated by
p
p
p
2D,i
= p
p
p
2D,p, f
+
d
f ,i
d
p, f ,i
p
p
p
2D,p,i
p
p
p
2D,p, f
. (18)
The derived initialization does not compensate po-
tential changes in the angle φ
f ,i
, caused by the differ-
ent spacings. Therefore, they are subject to a mapping
error φ
err, f ,i
and the correct condition is
φ
f ,i
= φ
p, f ,i
+ φ
err, f ,i
, (19)
resulting in the point positions
˜
p
p
p
2D,i
. So the error
caused by incorrect angles is
ε
i
=
p
p
p
2D,i
˜
p
p
p
2D,i
2
. (20)
With the help of trigonometric functions and Equa-
tions 16 and 20 it can be shown, that the error is
ε
i
= 2d
2
f ,i
sin
2
φ
err, f ,i
2
. (21)
Since the term sin
2
φ
err, f ,i
2
cannot be determined
without the solution, it is replaced by an unknown ex-
pected value which is assumed to be identical for all
combinations of f and i. Under this condition, the
fixed point, which minimizes the overall error
i
ε
i
, is
given by the mean value of all points
p
p
p
2D,opt, f
= |I
I
I
sub
|
1
|I
I
I
sub
|
i=1
p
p
p
2D,i
. (22)
A prototype lung should have a similar shape to
the processed lung so the surfaces can be well aligned
by the registration Tr
p
. That is why creating a mean
shape as a prototype could improve the initialization.
If the prototype and actual lung are too different, the
assumptions in this section might not hold. However,
the resulting imperfect initialization would not have
a negative influence on the mapping quality, but only
on the computation time.
3 RESULTS
To assess potential thickenings at an early stage,
high resolution chest CT data with a resolution
of 512 × 512 × 600 to 512 × 512 × 700 voxels is uti-
lized. The anisotropic image resolution varies be-
tween 0.63 mm and 1 mm in the transverse plane and
between 0.5 mm and 0.8 mm orthogonal to this plane.
Altogether, twenty data sets from ten different pa-
tients at two different points in time have been evalu-
ated. All calculations have been performed on an Intel
Core i5 with 16 GB of memory.
Some visual results for an exemplary patient are
given in Figure 6. The ground truth, extracted on the
volumetric data from a physician by labeling approx-
imately 550 2D slices, is shown in the first image 6a.
The following images 6c–6d show relevant scalar fea-
tures which have been extracted on the lung surface
for each surface point and Figure 6b shows the con-
fidence of an exemplary thickening classification, not
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma
241
0
20
40
60
80
100
120
140
160
0
0,05
0,1
0,15
0,2
0,25
0,3
0,35
0,4
0,45
16 32 48 64 80 96 112 128
duration/s
STRESS
r
STRESS1 STRESS2 duration
Figure 7: Distortion and computation time, depending on
radius r for the non-maxima suppression. The results are
given as mean values for all patients of left and right lung.
The error bars indicate the standard deviation.
(a) Differences with align-
ment (κ = 500).
-10 mm
10 mm
(b) Differences without
alignment (κ = 0).
Figure 8: Planar visualization of a lung at two points in time
from the same patient. The colors represent the differences
of the rib distances clipped to the interval [-10, 10] mm.
discussed in this publication. For a quantitative as-
sessment, the duration and the error of the flattening
as stress (Kruskal, 1964)
STRESS1 =
i> j
i, jI
I
I
(d
i, j
δ
i, j
)
2
i> j
i, jI
I
I
δ
2
i, j
0.5
, (23)
STRESS2 =
i> j
i, jI
I
I
(d
i, j
δ
i, j
)
2
i> j
i, jI
I
I
(d
i, j
mean (d
i, j
))
2
0.5
(24)
is plotted against the radius r of the NMS. The stress
measures indicate if the distances δ
i, j
on the 3D sur-
face are similar to the distances d
i, j
in the planar rep-
resentation. Low stress indicates similar distances
and high stress different distances. A vague verbal
scale of the STRESS1 measures was published by
Kruskal (Kruskal, 1964): starting with a values of 0
for perfect matching and ending with a value of 0.2
for poor matching. The average results for the left
and right lung of 10 patients at one point in time are
given in Figure 7.
For the temporal consistent flattening, a non-
rigid registration, based on B-splines (Rueckert et al.,
1999), was applied. This registration has already been
utilized for the automated assessment of the thicken-
ings and is available at no extra computational cost.
An exemplary visual comparison of two points in time
showing the effect of temporal consistency is shown
in Figure 8. The relation of slice-wise visualization
and planar visualization for a growing thickening is
shown in Figure 11. Quantitative results for the con-
sistent planar visualization are given in Figure 9. The
radii values (visualized by the different colors), lead-
ing to stable results for the planar visualization of a
single point in time, are investigated for the prob-
lem with two points in time with different connec-
tion strengths κ
{
0, 1, 500, 1000
}
(indicated by the
grouping). Additionally to the average STRESS1 cri-
terion, shown in Figure 9a, the average temporal dis-
tance between both data sets on the 2D plane
¯
d
2D
=
1
|
R
R
R
t
1
|
i=1...
|
R
R
R
t
1
|
min
i
0
=1...
|
R
R
R
t
2
|
d
2D
(p
p
p
2D,t
1
,i
, p
p
p
2D,t
2
,i
0
)
(25)
with
d
2D
(p
p
p
i
, p
p
p
j
) =
p
p
p
2D,t
1
,i
p
p
p
2D,t
2
, j
2
(26)
is given in Figure 9b to evaluate the consistency.
For the MDS initialization, a rough mapping to the
prototype is sufficient. Therefore, only a simple rigid
registration was applied which can be calculated in a
couple of seconds, while a non-rigid registration takes
typically a couple of minutes (Glocker et al., 2010).
An additional lung, not utilized in the evaluation, has
been chosen as the prototype. The calculation time
for the different components: adjacency matrix, dis-
similarity calculation by BFS and the MDS (with and
without initialization) are given in Figure 10. The
duration of the initialized MDS covers all additional
computations required, namely, the mapping of the
points Tr : p
p
p
3D,i
7→ p
p
p
3D,p,i
and the initialization itself.
4 DISCUSSION
Comparing the resulting computation time and the
quality of the flattening process depending on the ra-
dius, as shown in Figure 7, reveals a good compro-
mise of calculation time and error for r = 40. On aver-
age, the smallest evaluated radius of r = 16 results in
approximately 2350 surface points, the radius r = 40
results in approx. 410 surface points and a radius of
r = 128 results in only 14 surface points. For smaller
radii the computation time is significantly larger and
for larger radii the computation time is only slightly
decreasing, while the error is strongly increasing. The
stress stays relatively low even for large radii, which
indicates a good performance of the RBFs in interpo-
lating the lung shape. The relatively large acceptable
IVAPP 2016 - International Conference on Information Visualization Theory and Applications
242
0
0,02
0,04
0,06
0,08
0,1
0,12
0 1 500 1000
average STRESS1
κ
r=16 r=24 r=32 r=40
(a) The average STRESS1 results.
0
20
40
60
80
0 1 500 1000
κ
r=16 r=24 r=32 r=40
avgerage
d
2D
/pixel
130
110
90
(b) The average distances between the 2D representations
of both points in time.
Figure 9: Results for consistent mapping for various radii r
and various weightings κ.
radius of r = 40 might be surprising at first glance.
However the experiments showed that it is important
to map a few point optimally to the planar plane, and
all points in between can be interpolated with a simi-
lar quality as the mapping itself. An analogy would
be to thinly cover a 3D surface with liquid rubber.
After cooling down, the resulting rubber sheet can be
pinned to a planar surface at a few fixed points. If
the fixed points are chosen properly, the rubber sheet
arranges the surface in between these points auto-
matically with low stress. The stress asymptotically
approaches a minimum value larger zero, which is
caused by the shape of the lung itself which does not
allow a perfect isometric mapping.
As shown in Figure 9b, the temporal distance
between two points in time is significantly reduced
for κ = 1, but decreases only slightly for larger val-
ues κ {500, 1000}. As expected, both points in
time can be represented by a similar planar represen-
tation. Therefore, the problem defined by the matrices
in Equations 9 and 10 need only a small connection
0
10
20
30
40
50
60
70
80
90
16 20 24 28 32 36
duration/s
r
adjacency
matrix
dissimilarities
by BFS
MDS with
initialization
MDS without
initialization
Figure 10: Computation time for the different parts of the
algorithm, depending on the non-maxima suppression ra-
dius r.
(a) Slice view at point in
time t
1
.
(b) Slice view at point in
time t
2
.
(c) Planar representation
of thickening probability
at point in time t
1
.
(d) Planar representation
of thickening probability at
point in time t
2
.
Figure 11: Exemplary comparison of grown thickening (in-
dicated by orange circle) between two points in time t
1
and
t
2
utilizing the result of a pleura classification.
weight κ to reach a similar solution at both points in
time. The stable stress values for increasing connec-
tion strength κ, shown in Figure 9a, indicate that tem-
poral consistency has only a minor influence on the
individual flattening quality. The suggested method
results in planar representations of both points in time
which are much better comparable then before, with
minimal influence on the individual flattening. In
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma
243
cases of κ > 0, the stress slightly increases with in-
creasing radius r and at the same time, the temporal
distance is decreasing. An increasing radius r results
in a reduced number of constrains and therefore a less
accurate but more consistent flattening at both points
in time.
In the visualization of the non-aligned (κ = 0) sur-
faces in Figure 8b, a rotation between the rib struc-
tures is visible whereas the aligned (κ = 500) vi-
sualization reveals only minor differences between
the surfaces. The different alignment for the non-
connected MDS is mainly caused by the rotational
DOF. An exemplary comparison of a single thicken-
ing at two points in time is shown in Figure 11.
The initialization with the prototype reduces the
calculation time for the MDS approximately by a fac-
tor of 2, as shown by the blue and purple curve in Fig-
ure 10. For larger radii, the duration to calculate the
dissimilarities by BFS is getting more dominant. The
calculation of the adjacency matrix is independent of
the chosen radius and therefore its computation time
is constant.
5 CONCLUSIONS AND
OUTLOOK
The presented approach results in a planar visual-
ization of the lung surface, which enables a global
assessment of the pleura. Potential findings can be
exactly located and assessed regarding their extent.
With the newly introduced temporal link, which im-
proves the consistency of the visualizations from dif-
ferent points in time, the planar representation can
also be utilized for a follow-up assessment. An im-
portant result is that this consistency does not have a
considerable negative influence on the flattening qual-
ity. Additionally, with a closed-form initialization by
a prototype lung, the MDS can be sped up. The new
visualization can present simple scalar features. For a
correct judgment of pleural thickenings, the 3D view
is still required. However, the one-to-one mapping
of 2D points and 3D points enables an synchronized
3D volume and planar visualization. With this syn-
chronization, the user can benefit from both visual-
ization types, as motivated in Figure 11. For docu-
mentation purposes, the planar representation itself is
highly useful to record the final decision of the physi-
cian about the thickening locations.
Because of the rigid registration, only a rough
mapping of surface points to the prototype is possi-
ble. This limits the use of the prototype lung to the
initialization. With the presented approach, a perma-
nent link of all prototype points during optimization,
similar to the temporal link, would have a negative
impact on the flattening quality. However, a perma-
nent link without negative influence on the flattening
is desirable for inter-patient consistency and might be
an interesting topic for future research.
REFERENCES
Bentley, J. L. (1975). Multidimensional binary search
trees used for associative searching. Commun. ACM,
18(9):509–517.
Bronstein, A. M., Bronstein, M. M., and Kimmel, R.
(2003). Expression-invariant 3D face recognition. In
Audio-and Video-Based Biometric Person Authentica-
tion, pages 62–70. Springer.
Buja, A., Swayne, D. F., Littman, M. L., Dean, N., Hof-
mann, H., and Chen, L. (2008). Data visualization
with multidimensional scaling. Journal of Computa-
tional and Graphical Statistics, 17(2):444–472.
Carr, J. C., Fright, W. R., and Beatson, R. K. (1997). Sur-
face interpolation with radial basis functions for med-
ical imaging. IEEE Transactions on Medical Imaging,
16(1):96–107.
Chen, M., Helm, E., Joshi, N., and Brady, M. (2011).
Random walk-based automated segmentation for the
prognosis of malignant pleural mesothelioma. In
IEEE International Symposium on Biomedical Imag-
ing: From Nano to Macro, pages 1978–1981. IEEE.
Cox, T. F. and Cox, M. A. A. (2001). Multidimensional
Scaling. CRC Press.
Faltin, P., Chaisaowong, K., Ochsmann, E., and Kraus,
T. (2013). Vollautomatische Detektion und Ver-
laufskontrolle der pleuralen Verdickungen. In
53. Wissenschaftliche Jahrestagung der Deutschen
Gesellschaft f
¨
ur Arbeitsmedizin und Umweltmedizin
e.V. (DGAUM), Arbeitsmedizin in Europa, Muskel-
Skelett-Erkrankungen und Beruf, page 144, Bregenz
(Austria).
Floater, M. S. and Hormann, K. (2005). Surface parame-
terization: a tutorial and survey. In Advances in mul-
tiresolution for geometric modelling, pages 157–186.
Springer.
Frauenfelder, T., Tutic, M., Weder, W., G
¨
otti, R., Stahel,
R., Seifert, B., and Opitz, I. (2011). Volumetry: an
alternative to assess therapy response for malignant
pleural mesothelioma? European Respiratory Jour-
nal, 38(1):162–168.
Gleich, D. F. (2009). Models and Algorithms for PageRank
Sensitivity. PhD thesis, Stanford University. Chapter
7 on MatlabBGL.
Glocker, B., Komodakis, N., Paragios, N., and Navab, N.
(2010). Non-rigid registration using discrete mrfs:
Application to thoracic ct images. In Workshop Eval-
uation of Methods for Pulmonary Image Registration,
MICCAI 2010.
Grossmann, R., Kiryati, N., and Kimmel, R. (2002). Com-
putational surface flattening: a voxel-based approach.
IVAPP 2016 - International Conference on Information Visualization Theory and Applications
244
IEEE Transactions on Pattern Analysis and Machine
Intelligence, 24(4):433–441.
Hurdal, M. K. and Stephenson, K. (2009). Discrete confor-
mal methods for cortical brain flattening. Neuroimage,
45(1):S86–S98.
Kengne Nzegne, P., Faltin, P., Kraus, T., and Chaisaowong,
K. (2013). 3D lung surface analysis towards segmen-
tation of pleural thickenings. In Bildverarbeitung f
¨
ur
die Medizin 2013, pages 253–258. Springer.
Kruskal, J. B. (1964). Multidimensional scaling by opti-
mizing goodness of fit to a nonmetric hypothesis. Psy-
chometrika, 29(1):1–27.
Morse, B. S., Yoo, T. S., Rheingans, P., Chen, D. T., and
Subramanian, K. R. (2005). Interpolating implicit
surfaces from scattered surface data using compactly
supported radial basis functions. In Proc. ACM SIG-
GRAPH, pages 78–87. ACM.
Navon, I. and Legler, D. M. (1987). Conjugate-gradient
methods for large-scale minimization in meteorology.
Monthly Weather Review, 115(8):1479–1502.
Neubeck, A. and Gool, L. J. V. (2006). Efficient non-
maximum suppression. In Proc. of International Con-
ference on Pattern Recognition, pages 850–855.
Ochsmann, E., Carl, T., Brand, P., Raithel, H., and Kraus,
T. (2010). Inter-reader variability in chest radiogra-
phy and HRCT for the early detection of asbestos-
related lung and pleural abnormalities in a cohort of
636 asbestos-exposed subjects. INT ARCH OCC ENV
HEA, 83:39–46.
Pistolesi, M. and Rusthoven, J. (2004). Malignant pleu-
ral mesothelioma: Update, current management, and
newer therapeutic strategies. Chest, 126(4):1318–
1329.
Rudrapatna, M., Mai, V., Sowmya, A., and Wilson, P.
(2005). Knowledge-driven automated detection of
pleural plaques and thickening in high resolution CT
of the lung. In Information Processing in Medical
Imaging, pages 147–167. Springer.
Rueckert, D., Sonoda, L. I., Hayes, C., Hill, D. L. G., Leach,
M. O., and Hawkes, D. J. (1999). Nonrigid regis-
tration using free-form deformations: Application to
breast MR images. IEEE Transactions on Medical
Imaging, 18:712–721.
Saroul, L., Figueiredo, O., and Hersch, R. D. (2006). Dis-
tance preserving flattening of surface sections. IEEE
Transactions on Visualization and Computer Graph-
ics, 12(1):26–35.
Sensakovic, W. F., Armato III, S. G., Straus, C., Roberts,
R. Y., Caligiuri, P., Starkey, A., and Kindler, H. L.
(2011). Computerized segmentation and measurement
of malignant pleural mesothelioma. Medical physics,
38(1):238–244.
Tosun, D., Rettmann, M. E., and Prince, J. L. (2004).
Mapping techniques for aligning sulci across multiple
brains. Medical image analysis, 8(3):295–309.
Virta, R. L. (2012). 2012 Minerals Yearbook. U.S. Depart-
ment of the Interior, U.S. Geological Survey.
Zhao, X., Su, Z., Gu, X. D., Kaufman, A., Sun, J., Gao, J.,
and Luo, F. (2013). Area-preservation mapping using
optimal mass transport. Visualization and Computer
Graphics, IEEE Transactions on, 19(12):2838–2847.
Zigelman, G., Kimmel, R., and Kiryati, N. (2002). Tex-
ture mapping using surface flattening via multidimen-
sional scaling. Visualization and Computer Graphics,
IEEE Transactions on, 8(2):198–207.
Flattening of the Lung Surface with Temporal Consistency for the Follow-Up Assessment of Pleural Mesothelioma
245