2.25cm2.25cm2.25cm2.25cm

# Perfectly packing a cube by cubes of nearly harmonic sidelength

Rory McClenagan  Department of Mathematics and Statistics
University of Northern British Columbia
Prince George, BC V2N4Z9
Canada

###### Abstract.

Let $d$ be an integer greater than $1$, and let $t$ be fixed such that
$\frac{1}{d}<t<\frac{1}{d-1}$. We prove that for any $n_{0}$ chosen
sufficiently large depending upon $t$, the $d$-dimensional cubes of sidelength
$n^{-t}$ for $n\geq n_{0}$ can perfectly pack a cube of volume
$\sum_{n=n_{0}}^{\infty}\frac{1}{n^{dt}}$. Our work improves upon a previously
known result in the three-dimensional case for when $1/3<t\leq 4/11$ and
$n_{0}=1$ and builds upon recent work of Terence Tao in the two-dimensional
case.

Research of the author is partially supported by NSERC CGS-M

##  1\. Introduction

Let $d$ be an integer greater than $1$. We define a brick to be a closed
$d$-dimensional hyperrectangle and use the term cube to refer to a brick with
equal sidelengths. We define a packing of a finite or infinite collection of
bricks $\mathcal{B}$ to be a particular configuration of the bricks in
$\mathbb{R}^{d}$ such that the interior of the bricks are disjoint and the
facets of the bricks are parallel to the coordinate hyperplanes. A packing of
$\mathcal{B}$ in a solid $\Omega\subset\mathbb{R}^{d}$ is a packing of
$\mathcal{B}$ such that every brick is contained in $\Omega$. The packing is
perfect if the measure $m(\Omega\setminus\mathcal{B})$ is $0$. In this case
the sum of the volumes of the bricks must be equal to $m(\Omega)$.

A famous question posed by Meir and Moser [meir_moser-original] asks whether
rectangles of dimensions $\frac{1}{n}\times\frac{1}{n+1}$ for $n\geq 1$ can
perfectly pack a square of area $1$. They also ask whether squares of
dimensions $\frac{1}{n}\times\frac{1}{n}$ for $n\geq 2$ can perfectly pack a
square of area $\frac{\pi^{2}}{6}-1$. While both of these problems remain
open, there are two directions in which partial results have been obtained.

First, one can instead try to pack the same squares into a slightly larger
square. For instance, Paulhus [paulhus-squares] showed that the squares of
sidelength $1/n$ for $n\geq 2$ could be packed into a square of area
$\frac{\pi^{2}}{6}-1+\frac{1}{1244918662}$. However, it was pointed out in
[joos-rectangles] that the proof contained some errors; these errors were
corrected in [januszewski-correction].

Second, one can instead consider the problem of trying to perfectly pack the
squares of sidelength $1/n^{t}$ for $n\geq 1$ and some fixed $t>1/2$ into a
square of area $\sum_{n=1}^{\infty}\frac{1}{n^{2t}}$. This becomes harder as
$t\rightarrow 1^{-}$, and is obviously equivalent to the original problem when
$t=1$. Januszewski and Zielonka [januszewski-pack] verified this for
$1/2<t\leq 2/3$. At the expense of dropping the first few squares, Tao [tao-
pack] has recently proved that one could perfectly pack the collection of
squares of sidelength $1/n^{t}$ for $n\geq n_{0}$ (where $n_{0}$ depends on
$t$) into a square of area $\sum_{n=n_{0}}^{\infty}\frac{1}{n^{2t}}$ in the
entire range $1/2<t<1$.

The analogous problem for cubes in $\mathbb{R}^{3}$ has also been studied. For
instance, building on the methods presented in [joos-cubes], Januszewski and
Zielonka [januszewski-pack] have shown that cubes of sidelength $1/n^{t}$ for
$n\geq 1$ can be packed perfectly into a cube of volume
$\sum_{n=1}^{\infty}\frac{1}{n^{3t}}$ provided that $1/3<t\leq 4/11$. In the
general $d$-dimensional case, Joós showed in [joos-d-cubes] that $d$-cubes of
sidelength $n^{-t}$ for $n\geq 1$ can be packed into a cube of volume
$\zeta(dt)$ as long as $t$ is in the interval $[d_{0},2^{d-1}/(d2^{d-1}-1)]$,
where $d_{0}$ depends only on $d$. The lower bound, $d_{0}$, which Joós
implicitly defined, was later improved by Januszewski and Zielonka in
[januszewski-d-cubes] to $1/d$.

In this paper, we extend Tao’s work [tao-pack] in the $2$-dimensional case to
the $d$-dimensional case of cubes and prove the following result:

######  Theorem 1.1.

If $\frac{1}{d}<t<\frac{1}{d-1}$, and $n_{0}$ is sufficiently large depending
on $t$, then the cubes of sidelength $n^{-t}$ for $n\geq n_{0}$ can perfectly
pack a cube of volume $\sum_{n=n_{0}}^{\infty}\frac{1}{n^{dt}}$.

To prove Theorem 1.1, we apply an inductive-type argument similar to that used
by Tao in [tao-pack]. Initially, we suppose that we can pack a finite set of
cubes $\mathcal{C}$ of sidelength $n^{-t}$ for $n_{0}\leq n<n_{0}^{\prime}$
into our single cube $S$. As long as $S\setminus\mathcal{C}$ can be
partitioned into bricks $\mathcal{B}$ with small enough total surface area,
then we can find a brick $B\in\mathcal{B}$ which is wide enough to pack the
next cube of sidelength $(n_{0}^{\prime})^{-t}$. We pack $B$ by cubes
$\mathcal{C}^{\prime}$ of sidelength $n^{-t}$ for $n_{0}^{\prime}\leq
n<n_{0}^{{}^{\prime\prime}}$ in some efficient manner. By efficient, we mean
that the remaining space $B\setminus\mathcal{C}^{\prime}$ can be partitioned
into bricks $\mathcal{B}^{\prime}$ with small enough total surface area. In
the next iteration, we choose a wide brick from
$\mathcal{B}\setminus\\{B\\}\cup\mathcal{B}^{\prime}$ and pack it efficiently.
We proceed recursively until we have packed an arbitrarily large finite number
of cubes into $S$. Theorem 1.1 would then follow from a compactness argument.

This type of argument reduces the problem to finding a general technique for
packing cubes efficiently into some brick, in essence, forming the inductive
step in the above argument. Up until now, we have followed Tao’s argument in
[tao-pack] closely. However, while it is fairly straightforward to generalize
the standard two-dimensional packing algorithm used in [tao-pack] to the
higher-dimensional case, Tao’s method of explicitly verifying that this
packing is legal and efficient becomes much more difficult in three dimensions
or more. Our innovation is to introduce the notion of “snugness” (see Section
2); this allows us to perform this portion of the argument in an elegant
fashion which does not become too complex in the higher-dimensional setting.

In Section 2, we introduce our notation and prove some simple lemmas. In
Section 3, we reduce the proof of Theorem 1.1 to a more general result,
Proposition 3.1, which can be proved via induction. The inductive step of this
argument is furnished by Theorem 3.2 which provides a general and efficient
method for packing a brick by cubes. This result is proved in Section 4.

### Acknowledgements

We would like to thank Professor Terence Tao for the discussion on his blog
post describing his work in [tao-pack] and for answering our questions. We
would also like to thank Dr. Alia Hamieh for helpful comments and proof-
reading this paper, as well as Jaume de Dios, Dr. Rachel Greenfeld, and Dr.
Jose Madrid for helpful comments. Finally, we would like to thank the referee
for useful suggestions.

##  2\. Preliminary Lemmas and Notation

Throughout this paper, we will use the standard asymptotic notation $X=O(Y)$,
$X\ll Y$, and $Y\gg X$ to refer to the relation $X\leq C|Y|$. The constant $C$
will only be allowed to depend on the parameters $t$ and $\delta$, which we
will introduce shortly, but will be independent of all other parameters. If,
instead, we use the notation $X=O_{M}(Y)$ or $X\ll_{M}Y$, then the
corresponding constant $C$ is allowed to depend on the parameter $M$. We use
$X\asymp Y$ if $X\ll Y$ and $Y\ll X$. Finally, we use the notation $X=o(Y)$ if
$X/Y\rightarrow 0$ with respect to some explicit or implicit limiting
behaviour defined in context. We will implicitly assume that all of our
constants are allowed to depend upon the dimension $d$.

Note that we also use some non-standard asymptotic notation. If
$\bm{x}=(x_{1},x_{2},\dots,x_{d})\in\mathbb{R}^{d}$, then we use $\bm{x}+O(X)$
to refer to a vector $(x_{1}+O(X),\dots,x_{d}+O(X))$, and analogously for
little-$o$ notation. Similarly, if
$B=[B_{1},B_{1}^{\prime}]\times\dots\times[B_{d},B_{d}^{\prime}]$ is a brick
positioned in $\mathbb{R}^{d}$, then we use $B+O(X)$ to refer to a brick
$[B_{1}+O(X),B_{1}^{\prime}+O(X)]\times\dots\times[B_{d}+O(X),B_{d}^{\prime}+O(X)]$,
and analogously for little-$o$ notation.

Let $i,j\in\\{1,2,\dots,d\\}$. Given a brick $B$, we will denote its
sidelengths by $w_{i}(B)$, ordered so that $w_{i}(B)\leq w_{j}(B)$ for any
$i\leq j$. We say that the width of $B$ is the smallest sidelength, and denote
it by $w(B):=w_{1}(B)$. Clearly, $w_{i}(B)=w_{j}(B)$ for every $i$ and $j$ if
and only if $B$ is a cube. We define the volume of a single brick $B$ to be

$\mathrm{vol}(B):=w_{1}(B)w_{2}(B)\dots w_{d}(B).$

We define the eccentricity of a brick as

$\mathrm{ecc}(B):=\frac{\mathrm{vol}(B)}{w(B)^{d}}\geq 1.$

Note that $\mathrm{ecc}(B)=1$ if and only if $B$ is a cube.

Let $\mathcal{B}$ be a collection of bricks. Define the volume of
$\mathcal{B}$ to be

$\mathrm{vol}(\mathcal{B}):=\sum_{B\in\mathcal{B}}w_{1}(B)w_{2}(B)\dots
w_{d}(B).$

Define the unweighted surface area of $\mathcal{B}$ to be

$\mathrm{surf}(\mathcal{B}):=2\sum_{B\in\mathcal{B}}\sum_{1\leq
i_{1}<i_{2}<\dots<i_{d-1}\leq d}w_{i_{1}}(B)w_{i_{2}}(B)\dots
w_{i_{d-1}}(B)\asymp\sum_{B\in\mathcal{B}}w_{2}(B)w_{3}(B)\dots w_{d}(B).$

For $0\leq\delta<1$, define the weighted surface area of $\mathcal{B}$ to be

$\mathrm{surf}_{\delta}(\mathcal{B}):=\sum_{B\in\mathcal{B}}w_{1}(B)^{\delta}w_{2}(B)\dots
w_{d}(B).$

Clearly, $\mathrm{surf}_{0}(\mathcal{B})\asymp\mathrm{surf}(\mathcal{B})$.
Weighted surface area is roughly speaking a version of unweighted surface area
which weights high eccentricity bricks a little less than low eccentricity
bricks. We can use the inequality $w(B)\leq(w_{2}(B)w_{3}(B)\dots
w_{d}(B))^{\frac{1}{d-1}}$ for any brick $B$, to derive the crude bound

$\mathrm{surf}_{\delta}(\mathcal{B})\ll(\mathrm{surf}(\mathcal{B}))^{1+\frac{\delta}{d-1}},$
(2.1)

for a finite collection of bricks $\mathcal{B}$.

A solid $S\subset\mathbb{R}^{d}$ is called simple if it is connected and can
be formed as a union of a finite collection of bricks. A packing of a finite
collection of bricks $\mathcal{B}$ in a simple solid $S$ is called
$\varepsilon$-snug for some $\varepsilon>0$ if, for every brick
$B\in\mathcal{B}$, the portion of $\partial B$ which does not intersect the
boundary of another brick or the boundary of $S$ has surface area
$\ll(w\varepsilon)^{d-1}$ and the portion of $\partial S$ which does not
intersect the boundary of any brick in $\mathcal{B}$ also has surface area
$\ll(w\varepsilon)^{d-1}$. Here $w$ is the width of the widest brick in
$\mathcal{B}$. The size discrepancy of a finite collection of bricks
$\mathcal{B}$ is

$\mathrm{sd}(\mathcal{B})=\frac{\max_{B\in\mathcal{B}}w(B)}{\min_{B\in\mathcal{B}}w(B)}-1.$

The following lemma gives a criterion for the existence of a brick of a
certain minimum width in terms of an elegant relationship between volume and
weighted surface area.

######  Lemma 2.1.

Let $0\leq\delta<1$. For any finite collection of bricks $\mathcal{B}$, there
exists a brick with width at least
$\big{(}\frac{\mathrm{vol}(\mathcal{B})}{\mathrm{surf}_{\delta}(\mathcal{B})}\big{)}^{\frac{1}{1-\delta}}$.

###### Proof.

By definition,

$\mathrm{vol}(\mathcal{B})=\sum_{B\in\mathcal{B}}w_{1}(B)w_{2}(B)\dots
w_{d}(B)\leq\Big{(}\sup_{B\in\mathcal{B}}w(B)^{1-\delta}\Big{)}\mathrm{surf}_{\delta}(\mathcal{B}).$

This implies that
$(\mathrm{sup}_{B\in\mathcal{B}}w(B))^{1-\delta}\geq\frac{\mathrm{vol}(\mathcal{B})}{\mathrm{surf}_{\delta}(\mathcal{B})}$,
giving the desired result. ∎

This illustrates the principal behind using weighted surface area. If we were
to use unweighted surface area, namely setting $\delta=0$, then to guarantee
the existence of a brick of width $w$, we would need an upper bound on the
surface area of the form $w^{-1}\mathrm{vol}(\mathcal{B})$. However, if
$\delta>0$, then we only need a weaker bound of the form
$w^{-(1-\delta)}\mathrm{vol}(\mathcal{B})$.

The following result follows from a compactness argument (see, for example,
[martin-compactness]).

######  Lemma 2.2.

Let $\mathcal{B}$ be an, at most, countable collection of bricks and let
$\Omega\subset\mathbb{R}^{d}$ be compact. Suppose that an arbitrarily large,
but finite, number of bricks from $\mathcal{B}$ can be packed into $\Omega$.
Then $\mathcal{B}$ in its entirety can be packed into $\Omega$. Furthermore,
if $\mathrm{vol}(\mathcal{B})=\mathrm{m}(\Omega)$, then this packing is
perfect.

The following lemma states that if a packing of bricks is sufficiently snug,
then the region between the bricks has negligible surface area.

######  Lemma 2.3.

Suppose that a finite collection of bricks $\mathcal{B}$, where the widest
brick has width $w$, has a $\varepsilon$-snug packing in a brick $B$, for some
$\varepsilon>0$. Then, $B\setminus\mathcal{B}$ can be partitioned into bricks
with weighted surface area $\ll C_{|\mathcal{B}|}\nu w^{d-1+\delta}$, where
$C_{|\mathcal{B}|}$ is a constant that depends on $|\mathcal{B}|$ and
$\nu\rightarrow 0$ as $\varepsilon\rightarrow 0$.

###### Proof.

Partition $B\setminus\mathcal{B}$ into a finite number of bricks
$\mathcal{B}^{\prime}$. The maximum number of bricks in $\mathcal{B}^{\prime}$
can be bounded by a constant dependent upon $|\mathcal{B}|$. By the definition
of snugness, we know that the true surface area, $A$, (in the sense of the
($d-1$)-dimensional Lebesgue measure) of the solid $\cup\mathcal{B}^{\prime}$
could not exceed $(\varepsilon w)^{d-1}(|\mathcal{B}|+1)$. The result follows
from the crude bound $\mathrm{surf}(\mathcal{B^{\prime}})\ll
A|\mathcal{B}^{\prime}|$ and (2.1). ∎

##  3\. Initial Reductions

In this section we prove the higher-dimensional analogue of Proposition 2.1 in
[tao-pack] which will allow us to deduce Theorem 1.1.

######  Proposition 3.1.

Fix $\frac{1}{d}<t<\frac{1}{d-1}$ and $\delta$ depending on $t$, such that
$0<\delta<1$ and $(d-1)t+\delta t<1$. Choose a scale $M$ sufficiently large
and choose $N_{0}$ sufficiently large depending on $M$. Let
$n_{\mathrm{max}}\geq n_{0}\geq N_{0}$, and suppose that $\mathcal{B}$ is a
family of bricks with volume

$\mathrm{vol}(\mathcal{B})=\sum_{n=n_{0}}^{\infty}\frac{1}{n^{dt}},$ (3.1)

weighted surface area bound

$\mathrm{surf}_{\delta}(\mathcal{B})\ll\frac{1}{M^{1-\delta/2}}\sum_{n=1}^{n_{0}-1}\frac{1}{n^{(d-1)t+\delta
t}},$ (3.2)

and height bound

$\sup_{B\in\mathcal{B}}w_{d}(B)\ll 1.$ (3.3)

Then one can pack $\bigcup_{B\in\mathcal{B}}B$ by cubes of sidelength $n^{-t}$
for $n_{0}\leq n<n_{\mathrm{max}}$.

First we see how we can derive our main result from Proposition 3.1.

###### Proof of Theorem 1.1.

Fix $\delta=\frac{1}{d-1}-t$, and note that it easily satisfies the necessary
conditions. Take $\mathcal{B}=\\{C\\}$ where $C$ is the cube of volume
$\sum_{n=n_{0}}^{\infty}\frac{1}{n^{dt}}$, having sidelength $\ll
n_{0}^{1/d-t}$ (since $t>1/d$). Observe that (3.2) is satisfied, since

$\mathrm{surf}_{\delta}(\mathcal{B})\ll n_{0}^{(1/d-t)(d-1+\delta)}\ll
n_{0}^{1-dt}\ll\frac{1}{M^{1-\delta/2}}n_{0}^{1-(d-1)t-\delta
t}\ll\frac{1}{M^{1-\delta/2}}\sum_{n=1}^{n_{0}-1}\frac{1}{n^{(d-1)t+\delta
t}},$

recalling that $(d-1)t+\delta t<1$ and $0<\delta<1$. We also have used the
fact that $\frac{n_{0}^{t-\delta t}}{M^{1-\delta/2}}\gg 1$ since $n_{0}\geq
N_{0}$, which is sufficiently large depending upon $M$. Since (3.1) and (3.3)
are trivially satisfied, we can then apply Proposition 3.1 to conclude that
$C$ can be packed by cubes of sidelength $n^{-t}$ for $n_{0}\leq
n<n_{\mathrm{max}}$, and the result follows from Lemma 2.2. ∎

The inductive step in the proof of Proposition 3.1 requires us to pack a brick
by a collection of cubes. We isolate this result as a corollary to the
following more general theorem:

######  Theorem 3.2.

Fix $0\leq\delta<1$. Let $M=M_{1}\leq M_{2}\leq\dots\leq M_{d}$ be natural
numbers, and $\mathcal{C}$ be a family of $M_{*}=M_{1}M_{2}\dots M_{d}$ cubes
with maximum width $w$ and with size discrepancy $\varepsilon$, for some
$\varepsilon>0$. Let $S$ be a brick with dimensions $S_{1}\times
S_{2}\times\dots\times S_{d}$ satisfying $M_{i}w\leq S_{i}\leq M_{i}w+O(w)$
for $i\in\\{1,2,\dots,d\\}$. Then, there exists a packing of $\mathcal{C}$ in
$S$ such that $S\setminus\mathcal{C}$ can be partitioned into bricks
$\mathcal{B}$ satisfying

$\mathrm{surf}_{\delta}(\mathcal{B})\ll\frac{M_{*}}{M}w^{d-1+\delta}+C_{M}\nu
w^{d-1+\delta},$

where $C_{M}$ is some constant dependent upon $M$ and $\nu\rightarrow 0$ as
$\varepsilon\rightarrow 0$.

We will prove this theorem in Section 4. For now, we use it to derive the
corollary we need in the proof of Proposition 3.1:

######  Corollary 3.3.

Fix $\frac{1}{d}<t<\frac{1}{d-1}$ and $\delta$ depending on $t$, such that
$0<\delta<1$ and $(d-1)t+\delta t<1$. Choose a scale $M$ sufficiently large
and choose $N_{0}$ sufficiently large depending on $M$. Suppose that $S$ is a
brick satisfying the width bound $Mn_{0}^{-t}\leq w(S)\leq
Mn_{0}^{-t}+O(n_{0}^{-t})$ for some $n_{0}\geq N_{0}$ and satisfying the
eccentricity bound $\mathrm{ecc}(S)=o(n_{0})$. Then we can find
$n_{0}^{\prime}\geq n_{0}$ with
$n_{0}^{\prime}-n_{0}\asymp\mathrm{ecc}(S)M^{d}$, such that $S$ can be
perfectly packed by cubes of sidelength $n^{-t}$ for $n_{0}\leq
n<n_{0}^{\prime}$ and a collection of bricks $\mathcal{B}$ satisfying the
weighted surface area bound

$\mathrm{surf}_{\delta}(\mathcal{B})\ll\frac{1}{M}\sum_{n=n_{0}}^{n_{0}^{\prime}-1}\frac{1}{n^{(d-1)t+\delta
t}}.$

###### Proof.

Let the parameters be chosen as in the corollary, and use the notation
$S=S_{1}\times S_{2}\times\dots\times S_{d}$ such that $S_{1}\leq
S_{2}\leq\dots\leq S_{d}$. Thus, $Mn_{0}^{-t}\leq S_{1}\leq
Mn_{0}^{-t}+O(n_{0}^{-t})$. Define $M_{i}=\lfloor S_{i}/n_{0}^{-t}\rfloor$ for
$i\in\\{1,2,\dots,d\\}$, so that $M_{1}\asymp M$. Choose
$n_{0}^{\prime}=n_{0}+M_{*}$. Note that

$M_{*}\asymp\mathrm{ecc}(S)M^{d},$ (3.4)

as required. Let $\mathcal{C}$ be the collection of cubes of sidelength
$n^{-t}$ for $n_{0}\leq n<n_{0}^{\prime}$. The size discrepancy is
$\frac{n_{0}^{\prime t}}{n_{0}^{t}}-1$. This can be made arbitrarily small as
long as $N_{0}$ is chosen to be sufficiently large compared with $M$. This
makes the second term in the bound of Theorem 3.2 negligible with respect to
the first. Thus, we can apply Theorem 3.2 to get a packing of $\mathcal{C}$ in
$S$ such that $S\setminus\mathcal{C}$ can be partitioned into bricks
$\mathcal{B}$ satisfying

$\mathrm{surf}_{\delta}(\mathcal{B})\ll\frac{M_{*}}{M}(n_{0}^{-t})^{d-1+\delta}\ll\frac{M_{*}}{M}\frac{(1+\mathrm{sd}(\mathcal{C}))^{d-1+\delta}}{(n_{0}^{\prime})^{(d-1)t+\delta
t}}\ll\frac{1}{M}\sum_{n=n_{0}}^{n_{0}^{\prime}-1}\frac{1}{n^{(d-1)t+\delta
t}},$

since $\mathrm{sd}(\mathcal{C})\rightarrow 0$. This completes the proof. ∎

Observe that the power of $M$ in the weighted surface area bound of the
corollary is independent of $\delta$. This fact allows us to loosen our
weighted surface area bound (3.2) by a factor of $M^{\delta/2}$, which is
enough to let us complete the inductive step of the proof of Proposition 3.1,
illustrating the advantage of working with weighted surface area (see also the
discussion after Lemma 2.1).

We now use this corollary to prove Proposition 3.1. Our proof closely mirrors
the proof of Proposition 2.1 in [tao-pack] except for higher dimensions.
However, for the reader’s convenience, we include it here.

###### Proof of Proposition 3.1.

We prove this via downward induction on $n_{0}$. Fix $n_{\mathrm{max}}\geq
N_{0}$. Clearly, the result holds if $n_{0}=n_{\mathrm{max}}$. Fix some
$n_{0}\leq n_{\mathrm{max}}$, and assume the result holds with $n_{0}$
replaced by any strictly larger integer up to $n_{\mathrm{max}}$. We show that
the result will then hold for $n_{0}$.

Since $t>1/d$, (3.1) implies that $\mathrm{vol}(\mathcal{B})\asymp
n_{0}^{1-dt}$. Furthermore, since $(d-1)t+\delta t<1$, we have
$\mathrm{surf}_{\delta}(\mathcal{B})\ll
M^{-(1-\delta/2)}n_{0}^{1-(d-1)t-\delta t}$. Thus, Lemma 2.1 implies the
existence of a brick $B^{\prime}\in\mathcal{B}$ satisfying

$w(B^{\prime})\gg\bigg{(}\frac{n_{0}^{1-dt}}{M^{-(1-\delta/2)}n_{0}^{1-(d-1)t-\delta
t}}\bigg{)}^{\frac{1}{1-\delta}}=M^{\frac{1-\delta/2}{1-\delta}}n_{0}^{-t}.$

Since $\frac{1-\delta/2}{1-\delta}>1$ for $0<\delta<1$, then as long as we
take $M$ sufficiently large, we can drop the implied constant and conclude
that $w(B^{\prime})\geq Mn_{0}^{-t}$. Partition $B^{\prime}$ into two bricks
$B$ and $B^{\prime}\setminus B$, so that $Mn_{0}^{-t}\leq w(B)\leq
Mn_{0}^{-t}+O(n_{0}^{-t})$. By the height bound, (3.3), $\mathrm{ecc}(B)\ll
M^{-(d-1)}n_{0}^{(d-1)t}=o(n_{0})$, which means that we can apply Corollary
3.3 to pack $B$ by cubes of sidelengths $n^{-t}$ for $n_{0}\leq
n<n_{0}^{\prime}$ with $n_{0}^{\prime}-n_{0}\gg M^{d}$ and a collection of
bricks $\mathcal{B}_{0}$ satisfying

$\mathrm{surf}_{\delta}(\mathcal{B}_{0})\ll\frac{1}{M}\sum_{n=n_{0}}^{n_{0}^{\prime}-1}\frac{1}{n^{(d-1)t+\delta
t}}.$ (3.5)

Now, if $n_{0}^{\prime}\geq n_{\mathrm{max}}$, then we are done, as we have
packed every cube of sidelengths $n^{-t}$ for $n_{0}\leq n<n_{\mathrm{max}}$.
Otherwise, suppose that $n_{0}^{\prime}<n_{\mathrm{max}}$. Since
$n_{0}^{\prime}$ is strictly larger than $n_{0}$, it makes sense to now apply
our inductive hypothesis, replacing $n_{0}$ by $n_{0}^{\prime}$ (which is
strictly larger than $n_{0}$), and replacing $\mathcal{B}$ by
$\mathcal{B}^{\prime}=(\mathcal{B}\setminus\\{B^{\prime}\\})\cup\\{B^{\prime}\setminus
B\\}\cup\mathcal{B}_{0}$. First, however we have to check to assure that the
conditions of the proposition are met. Observe,

$\mathrm{vol}(\mathcal{B}^{\prime})=\sum_{n=n_{0}}^{\infty}\frac{1}{n^{dt}}-\sum_{n=n_{0}}^{n_{0}^{\prime}-1}\frac{1}{n^{dt}}=\sum_{n=n_{0}^{\prime}}^{\infty}\frac{1}{n^{dt}},$

and so $\mathcal{B}^{\prime}$ has the required total volume (3.1). Clearly,
$\mathcal{B}^{\prime}$ satisfies the height bound (3.3). Finally, by (3.5),
(3.2), and the fact that $\mathrm{surf}_{\delta}\\{B^{\prime}\setminus
B\\}\leq\mathrm{surf}_{\delta}\\{B^{\prime}\\}$, we have

$\mathrm{surf}_{\delta}(\mathcal{B}^{\prime})\leq\mathrm{surf}_{\delta}(\mathcal{B}\cup\mathcal{B}_{0})\ll\frac{1}{M^{1-\delta/2}}\sum_{n=1}^{n_{0}-1}\frac{1}{n^{(d-1)t+\delta
t}}+\frac{1}{M}\sum_{n=n_{0}}^{n_{0}^{\prime}-1}\frac{1}{n^{(d-1)t+\delta
t}}\ll\frac{1}{M^{1-\delta/2}}\sum_{n=1}^{n_{0}^{\prime}-1}\frac{1}{n^{(d-1)t+\delta
t}},$

and thus $\mathcal{B}^{\prime}$ satisfies the weighted surface bound (3.2).
Thus, we can apply the inductive hypothesis for $n_{0}^{\prime}$ and pack
$\bigcup_{B\in\mathcal{B}^{\prime}}B$ by the remaining cubes of sidelength
$n^{-t}$ for $n_{0}^{\prime}\leq n<n_{\mathrm{max}}$, which in turn implies
that we can pack $\bigcup_{B\in\mathcal{B}}B$ by cubes of sidelength $n^{-t}$
for $n_{0}\leq n<n_{\mathrm{max}}$. ∎

All that remains is proving Theorem 3.2 which provides a general and efficient
brick-packing algorithm.

##  4\. Efficient Brick-Packing Algorithm

###### Proof of Theorem 3.2.

Note that, without loss of generality, we can assume that $S_{i}=M_{i}w$ for
$i\in\\{1,2,\dots,d\\}$. To see this, suppose that $S^{\prime}$ is a cube that
contains $S$ and instead satisfies $M_{i}w\leq S_{i}^{\prime}\leq
M_{i}w+O(w)$. Then, $S^{\prime}\setminus S$ can be partitioned into a $O(1)$
bricks, each of which contributes an allowable weighted surface area
$\ll\frac{M_{*}}{M}w^{d-1+\delta}$.

To explicitly define our packing, we position $S$ in $\mathbb{R}^{d}$ as

$[0,M_{1}w]\times[0,M_{2}w]\times\dots\times[0,M_{d}w].$

Index $\mathcal{C}$ as $\\{C_{n}\\}_{n=0}^{M_{*}-1}$ from largest width to
smallest width. We use the notation $w_{n}:=w(C_{n})$. By construction,
$w_{n}\leq w_{m}$ if and only if $n\geq m$. We further index
$n=0,1,\dots,M_{*}-1$ by $n_{\bm{i}}=n_{i_{1},i_{2},\dots,i_{d}}$, where

$n_{i_{1},i_{2},\dots,i_{d}}:=i_{1}+i_{2}M_{1}+i_{3}M_{1}M_{2}+\dots+i_{d}M_{1}M_{2}\dots
M_{d-1},$

for $i_{k}=0,\dots,M_{k}-1$ (with $k=1,2,\dots,d$). We use the notation
$C_{\bm{i}}:=C_{n_{\bm{i}}}$ and $w_{\bm{i}}:=w_{n_{\bm{i}}}$. Position each
$C_{\bm{i}}$ in $S$ as

$C_{\bm{i}}:=[x^{1}_{\bm{i}},x^{1}_{\bm{i}}+w_{\bm{i}}]\times[x^{2}_{\bm{i}},x^{2}_{\bm{i}}+w_{\bm{i}}]\times\dots\times[x^{d}_{\bm{i}},x^{d}_{\bm{i}}+w_{\bm{i}}],$

where for any $k\in\\{1,2,\dots,d\\}$, we define

$x^{k}_{\bm{i}}=\sum_{i_{k}^{\prime}=0}^{M_{k}-1}w_{i_{1},\dots,i_{k-1},i_{k}^{\prime},0,\dots,0}-\sum_{i_{k}^{\prime}=i_{k}}^{M_{k}-1}w_{i_{1},\dots,i_{k-1},i_{k}^{\prime},i_{k}+1,\dots
i_{d}}$

(see Figure 1).

Figure 1. The packing of the cubes $C_{\bm{i}}$ in $S$. Here, $d=2$,
$M_{1}=3$, $M_{2}=4$, $M_{*}=12$, and $i_{1}=i_{2}=1$. Note that the diagram
is not to scale.

We will verify that this is a legal packing shortly. Note that each
$(x_{\bm{i}}^{1},x_{\bm{i}}^{2},\dots,x_{\bm{i}}^{d})$ is asymptotically fixed
at a lattice point as $\mathrm{sd}(\mathcal{C})\ll\varepsilon$, namely

$(x_{\bm{i}}^{1},x_{\bm{i}}^{2},\dots,x_{\bm{i}}^{d})=(wi_{1},wi_{2},\dots,wi_{d})+O_{M}(\varepsilon
w).$ (4.1)

Collect the subset of cubes from $\mathcal{C}$ which form its exterior
“shell”:

$\tilde{\mathcal{C}}=\\{C_{\bm{i}}\in\mathcal{C}:i_{k}=0\text{ or
}i_{k}=M_{k}-1,\,\text{for some $k\in\\{1,2,\dots,d\\}$}\\}.$

Let $B$ be the smallest brick containing
$\mathcal{C}\setminus\tilde{\mathcal{C}}$. By (4.1), $B$ has dimensions

$[w,(M_{1}-1)w]\times[w,(M_{2}-1)w]\times\dots\times[w,(M_{d}-1)w]+o(1).$

Figure 2. The simple solid $K$ is constructed from $B$ and $\mathcal{C}$.

Define the simple solid $K=B\cup\mathcal{C}=B\cup\tilde{\mathcal{C}}$ (see
Figure 2). Observe that

$S\setminus K=(S\setminus B)\setminus\tilde{\mathcal{C}}.$

Observe that $S\setminus B$ can be partitioned into $O(M_{*}/M)$ bricks
$\mathcal{B}^{\prime}$ each with dimensions $O(w)$. Each brick
$B^{\prime}\in\mathcal{B}^{\prime}$ intersects at most $O(1)$ cubes in
$\tilde{\mathcal{C}}$. This means that we can partition
$B^{\prime}\setminus\tilde{\mathcal{C}}$ into $O(1)$ bricks, each with
weighted surface area less than that of $B^{\prime}$, which is $\ll
w^{d-1+\delta}$. Thus, $S\setminus K$ can be partitioned into bricks with
allowable weighted surface area $\ll\frac{M_{*}}{M}w^{d-1+\delta}$. It then
suffices to show that $K\setminus\mathcal{C}$ can be partitioned into bricks
with weighted surface area $\ll\frac{M_{*}}{M}w^{d-1+\delta}$.

By Lemma 2.3, it suffices to show that $\mathcal{C}$ is a packing that is
$O_{M}(\varepsilon)$-snug in $K$. First, observe that all of the cubes are
inside of $S$. This follows from the bound $w_{\bm{i}}\leq w$. Thus, we have
to check that none of the cubes’ interiors overlap, and that every cube in
$\mathcal{C}\setminus\tilde{\mathcal{C}}$ is touching the $2d$ adjacent cubes
(the cubes in $\tilde{\mathcal{C}}$ are already touching $K$ by construction).

Define $\pi_{k}$ to be the projection operator onto the $x_{k}$ axis for every
$k\in\\{1,2,\dots,d\\}$. Let $\tilde{E}$ be the collection of $2^{d}-1$
vectors $\bm{e}=(e_{1},e_{2},\dots,e_{d})$ such that every $e_{k}\in\\{0,1\\}$
but $\bm{e}\neq\bm{0}$. Let $\bm{e}_{k}$ be the $k$th unit vector, namely
$(0,\dots,0,1,0,\dots,0)$, with a $1$ in the $k$th component and let
$E\subset\tilde{E}$ be the collection of such $d$ unit vectors. Define $I$ be
the collection of $\bm{i}=(i_{1},i_{2},\dots,i_{d})$ such that
$i_{k}\in\\{0,1,\dots,M_{k}-2\\}$ for $k\in\\{1,2,\dots,d\\}$. By symmetry and
the asymptotic positioning of the cubes (4.1), we only have to worry about
checking overlap on “adjacent” cubes, reducing the proof to showing the
following two claims:

  1. (i)

Let $\bm{e}\in E$ and $\bm{i}\in I$. Then for at least one
$k\in\\{1,2,\dots,d\\}$ we have that
$\pi_{k}(C_{\bm{i}})\cap\pi_{k}(C_{\bm{i}+\bm{e}})$ is exactly one point.

  2. (ii)

Let $\bm{e}\in\tilde{E}\setminus E$ and $\bm{i}\in I$. Then for at least one
$k\in\\{1,2,\dots,d\\}$, we have that
$\pi_{k}(C_{\bm{i}})\cap\pi_{k}(C_{\bm{i}+\bm{e}})$ is at most one point.

To see why this is sufficient to complete the proof, observe that as long as
the boundary of the cubes are touching, the asymptotic positioning of the
cubes (4.1) ensures that the non-overlapping boundary will have area
$O_{M}(\varepsilon w)$, meaning that (i) will imply that the packing is
$O_{M}(\varepsilon)$-snug. Clearly, (ii) implies that none of the cubes’
interiors overlap, and thus our packing of $\mathcal{C}$ is valid.

To see (i), note that the construction of the $C_{i}$ immediately implies that
for every $k\in\\{1,2,\dots,d\\}$, we have

$\pi_{k}(C_{\bm{i}+\bm{e}_{k}})\cap\pi_{k}(C_{\bm{i}})=\\{x_{\bm{i}}^{k}+w_{\bm{i}}\\}.$

Now we show (ii). Fix some
$\bm{e}=(e_{1},e_{2},\dots,e_{d})\in\tilde{E}\setminus E$. Let $k$ be the
smallest index such that the component $e_{k}$ is nonzero. Clearly,
$k\in\\{1,2,\dots,d-1\\}$. Recall that

$x^{k}_{\bm{i}}=\sum_{i_{k}^{\prime}=0}^{M_{k}-1}w_{i_{1},\dots,i_{k-1},i_{k}^{\prime},0,\dots,0}-\sum_{i_{k}^{\prime}=i_{k}}^{M_{k}-1}w_{i_{1},\dots,i_{k-1},i_{k}^{\prime},i_{k}+1,\dots
i_{d}}.$

By the ordering of $w_{\bm{i}}$, we have that $w_{\bm{i}+\bm{e}}\leq
w_{\bm{i}+\bm{e}_{k}}$. Thus, $x_{\bm{i}+\bm{e}}\geq
x_{\bm{i}+\bm{e}_{k}}=x_{\bm{i}}+w_{\bm{i}}$. This shows that
$\pi_{k}(C_{\bm{i}+\bm{e}})\cap\pi_{k}(C_{\bm{i}})$ is at most a singleton, as
desired. This completes the proof. ∎

## References


# Deep Learning Methods for Retinal Blood Vessel Segmentation: Evaluation on
Images with Retinopathy of Prematurity

1st Gorana Gojić  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>2nd Veljko Petrović  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>3rd Radovan Turović  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>4th Dinu Dragan  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>5th Ana Oros  Dept. of Retinopathy of Prematurity and Retinal
Development
Institute of Neonatology
Belgrade, Serbia
<EMAIL_ADDRESS>6th Dušan Gajić  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>7th Nebojša Horvat  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>
###### Abstract

Automatic blood vessel segmentation from retinal images plays an important
role in the diagnosis of many systemic and eye diseases, including retinopathy
of prematurity. Current state-of-the-art research in blood vessel segmentation
from retinal images is based on convolutional neural networks. The solutions
proposed so far are trained and tested on images from few available retinal
blood vessel segmentation datasets, which might limit their performance when
given an image with retinopathy of prematurity signs. In this paper we
evaluate performance of three high-performing convolutional neural networks
for retinal blood vessel segmentation in context of blood vessel segmentation
on retinopathy of prematurity retinal images. The main motive behind the study
is to test if existing public datasets suffice to develop a high performing
predictor that could assist an ophthalmologist in retinopathy of prematurity
diagnosis. To do so, we create a dataset consisting solely of retinopathy of
prematurity images with retinal blood vessel annotations manually labeled by
two observers, where one is the ophthalmologist experienced in retinopathy of
prematurity treatment. Experimental results show that all three solutions have
difficulties in detecting the retinal blood vessels of infants due to a lower
contrast compared to images from public datasets as demonstrated by
significant drop in classification sensitivity. All three solutions segment
alongside retinal also choroidal blood vessels which are not used to diagnose
retinopathy of prematurity, but instead represent noise and are confused with
retinal blood vessels. By visual and numerical observations, we observe that
existing solutions for retinal blood vessel segmentation need improvement
toward more detailed datasets or deeper models in order to assist the
ophthalmologist in retinopathy of prematurity diagnosis.

###### Index Terms:

Retinal blood vessel, segmentation, retinopathy of prematurity, convolutional
neural networks, retinal image

##  I Introduction

Retinopathy of prematurity (ROP) [1] is an eye disease affecting prematurely
born infants. It manifests itself through anomalies in the retina which is the
back part of the eye. Abnormal chemical changes in the eye result in
interruption of normal blood vessel development causing excessive vessel
growth in affected areas and appearance of vessel loops as can be seen in Fig.
1. Since the normal vessel growth is suspended, parts of the retina stay
avascularized, and if not treated, detach over time and finally lead to
blindness.

To effectively treat ROP, and prevent the possibility of disease developing
into blindness, it is of great importance to diagnose it in early stages of
development. Prematurely born infants fulfilling certain criteria are screened
for ROP to check for the retinal blood vessel tortuosity and dilatation, which
are initial ROP symptoms. These changes can be observed in digital images
captured by specialized wide-field fundus cameras such as RetCam. In the first
column of Fig. 3 we show retinal images captured with a RetCam3 camera.
However, ROP diagnosis from an image can be challenging due to issues such as:
(1) low contrast between blood vessels and retina, (2) variable illumination
in image due to wide-field view, (3) relatively low image resolution and (4)
appearance of choroidal blood vessels which might be visible on image but do
not play part in ROP diagnosis [2].

Figure 1: Eye retinal image showing excessive blood vessel growth and
appearance of closed blood vessel paths.

To assist ophthalmologists in disease diagnosis from digital retinal images,
many general, fully-automated blood vessel segmentation solutions have been
developed. Since the subset of ROP symptoms manifests itself through blood
vessel changes, it can be expected that these solutions could assist an
ophthalmologist in disease diagnosis. Current state-of-the-art solutions for
retinal blood vessel segmentation are based on convolutional neural networks
(CNNs) [3, 4, 5, 6, 7, 8, 9]. It is reported that these solutions achieve
outstanding segmentation results confirmed both visually and numerically.
However, it is known that CNNs performance depends on the data used for
training. All high performing solutions in this field are trained using
publicly available retinal segmentation datasets counting up to 40 annotated
retinal images of adults or school children. We observe that images from
publicly available datasets are visually different and apparently easier to
segment compared to the retinal images of infants showing signs of ROP. Some
differences between images from those datasets and ROP images concern: (1) a
blood vessel network which is fully developed in images from publicly
available datasets occasionally showing signs of degradation, while ROP images
show just partially developed blood network, (2) blood vessels that are
thicker compared to infant blood vessels, (3) a contrast between blood vessels
and retina which is usually stronger, and (4) the illumination distribution
which is uneven in wide-field ROP images opposite to narrow-field images from
the public datasets.

Although the data class remains the same, existing solutions trained on public
datasets might degrade in performance when tested on ROP images. It is our
intention to measure the impact on the network performance when given a ROP
image and to determine if the existing solutions can be used as-is to assist
in ROP diagnosis. Since these solutions perform well when tested on publicly
available datasets, poor results in this study could indicate that existing
publicly available datasets do not suffice to develop high-performing CNN that
can be used for blood vessel segmentation on ROP images. To test the
hypothesis, we choose three recently developed, high-performing CNNs proposed
by Li et al. [9], Oliveira et al [4], and Zhuang [10]. We feed pretrained
models with 9 ROP images provided by the Institute of Neonatology in Belgrade.
The blood vessels on the images are manually labeled by two observers,
including one ROP expert with more than 20 years of experience in ROP
treatment, to create ground truth labels which are used to evaluate the
selected networks.

The rest of the paper is organized as follows: in Section II we give an
overview of the datasets commonly used to train CNNs for blood vessel
segmentation. Also, we introduce the reader to the recent progress in the
development of CNNs for blood vessel segmentation from color retinal images.
Section III discusses the details of the methodology used, including the
preparation of a ROP dataset, the experimental protocol, and used evaluation
metrics. Obtained results are presented and discussed in Section IV. The final
section offers the main conclusions as well as directions for the future work.

##  II Related Work

In this section we give a brief overview of publicly available datasets for
blood vessel segmentation from retinal images. We include DRIVE [11], STARE
[12] and CHASE_DB1 [13] datasets, which are mainstream datasets for automatic
vessel extraction, as well as less known HFR [14] image dataset. Example
images from these four datasets and ROP image are shown in Fig. 2. Next, we
introduce the reader to the recent work on automatic blood vessel segmentation
using convolutional neural networks. For comprehensive review of methods for
retinal vessel segmentation, including both machine learning approaches not
based on neural networks and non-machine learning approaches, we direct the
reader to the work of Jin et al. [3] and Oliveira et al. [4]. Here we discuss
solely CNN based solutions since by quantitative measurements they have proven
to be the most effective in blood vessel segmentation from a color retinal
image.

Figure 2: Example images from retinal vessel segmentation datasets. From left
to right are representatives of: (1) DRIVE, (2) STARE, (3) CHASE_DB1, (4) HFR,
and (5) ROP image.

###  II-A Retinal datasets for blood vessel segmentation

The DRIVE dataset [11] includes 40 color retinal images measuring $768\times
584$ pixels collected as a part of a Netherlands screening program for
diabetic retinopathy. Images are captured using a non-mydriatic 3CCD camera at
45° FoV. All images in the dataset belong to adults of age 25 to 90, with a
fully developed blood vessel network. Seven images show signs of mild vessel
degradation as a result of diabetic retinopathy. All images are accompanied
with manually labeled groundtruth vessel segmentation masks and FoV masks. One
groundtruth set is labeled by an expert in ophthalmology, while the other is
labeled by an observer trained by the ophthalmologist.

The STARE dataset [12] contains 20 color retinal images measuring $700\times
605$ pixels extracted from the larger database created as a part of STARE
project. The images are obtained using Topcon TRV-50 fundus camera at 35° FoV.
Ten images show signs of different diseases that distort or obscure the blood
vessel appearance, but none of the images show signs of retinopathy of
prematurity. For the full diagnosis list, the reader is referred to the STARE
project’s page. Remaining 10 images show no signs of disease. Similarly to
DRIVE, this dataset comes with two sets of manually annotated groundtruth
blood vessel masks.

CHASE_DB1 [13] dataset contains 28 color retinal images measuring $1280\times
960$ pixels. The images are captured by the Nidek NM-200-D handheld fundus
camera with a 30° FoV. Initially, the dataset was created to validate the
performance of a novel computer software for retinal vessel tortuosity
measurements in school-aged children. In the dataset, there are images of left
and right eyes for 14 different children, each image with two sets of
groundtruth blood vessel masks. Since tortuosity is one of the key factors in
diagnosis of retinopathy of prematurity, it is expected that networks trained
on this dataset could perform better on ROP images.

The HFR dataset [14] includes 45 color retinal images of size $1280\times 960$
pixels. The images are captured using a Canon CF-60UVi camera. The dataset
includes 15 images of healthy patients, 15 glaucoma images and 15 images
showing diabetic retinopathy. Images are accompanied with two manually labeled
groundtruth blood vessel masks, as well as FoV masks.

###  II-B Convolutional neural networks for blood vessel segmentation from
retinal images

The publication of UNet architecture [5] for medical image segmentation
boosted development of CNN based solutions in many medical applications
including retinal blood vessel segmentation [6, 3, 4, 7, 10, 9]. Confronting
the need for a large dataset to train a well performing CNN, Ronneberger et
al. have demonstrated that the UNet architecture can produce state-of-the-art
results for image segmentation task, while being trained on a small dataset.
Availability of small public datasets for retinal blood vessel segmentation
directed the efforts toward CNNs based on UNet architecture. There have been
many modifications of the original UNet architecture to address different
problems occurring in blood vessel segmentation, such as varying vessel width,
shape and orientation, low image quality, as well as vessel branching and
crossings which often result in misclassification of a single vessel as two
separate vessels. In [6] authors address the appearance of thin false vessels
by incorporating UNet as a generator of a Generative Adverserial Network
(GAN). Jan et al. [3] proposed the network inspired by UNet and Deformable
convolutional networks [15] to address variability in shapes and scales of
blood vessels. They replace the convolutional layer with a deformable
convolutional block resulting in adaptive adjustment of a receptive field to
the vessels’ scales and shapes during network training. Similar issues have
been addressed in [4] where authors cope with vessels’ scale and orientation
variation by calculating Stationary Wavelet Transform (SWT) and using it as
additional information about vessels in network training. In [7] authors
propose an addition to a CNN named Graph Convolutional Network [8] (GCN).
While the CNN filters learn local per-pixel features, the GCN learns the
global structure of the vessel network using the fact that the vessel network
can be represented as a graph. The final prediction is derived by combining
CNN and GCN learned features. Another approach to improve segmentation is to
chain multiple networks as demonstrated in [10, 9]. The author of the
LadderNet [10] propose UNet modification where multiple encoders and decoders
with skip connections are chained to improve data flow, resulting in a network
able to better learn features from the input image. Similar idea is employed
with IterNet [9] where multiple networks are chained to iteratively improve
the initial prediction. The first network in a chain is the UNet which
produces the initial prediction for the given input color image. The output of
the second last layer is then iteratively refined by lightweight UNets and the
final prediction is the output of the last network in a chain.

All mentioned solutions use DRIVE in test and training stages. Very often,
STARE [6] or CHASE_DB1 [10] datasets, or both [3, 4, 9], are used alongside
DRIVE to create and evaluate the proposed network. The exception is [3] where
WIDE [16] and SYNTHE [17] datasets are used to additionally test the proposed
network. These datasets are not often used since they do not have groundtruth
labels and the second dataset is synthetically generated. To measure the
network performance, authors often report network’s prediction accuracy [3, 4,
10, 9], sensitivity [3, 4, 10, 9], specificity [3, 4, 10, 9] and precision [6,
9], as well as derived measures. To be comparable to previously acquired
results, we also adopt these measures in our study.

##  III Methodology

In this section we report details of the dataset used for evaluation,
experimental protocol, and measures used for network performance evaluation.

Figure 3: Images used for the CNN evaluation. From left to right: (1) a color
image showing sings of ROP, (2) FoV mask, (3) the first groundtruth, (4) the
second groundtruth.

###  III-A The dataset

The dataset used in this experiment consists of 9 color retinal images shown
in the first column of Fig. 3. To capture the images, a wide-field RetCam3
fundus camera with 130° FoV was used. Resulting images are of size $480\times
640$ and are saved in JPEG file format. Abnormal vessel tortuosity and
dilatation are present on all images. To test the networks on complex vessel
patterns, we include images with visible excessive vessel growth and presence
of vessel loops (e.g. rows 3, 5, 6 and 9 in Fig. 3). For each color image, we
provide two groundtruth blood vessel masks (columns 3 and 4 in Fig. 3) and FoV
mask (the second column in Fig. 3). The groundtruth masks are labeled by two
different observers. The first groundtruth is labeled by a PhD student in
computing and control engineering which was previously trained by the
experienced ophthalmologist. This is not a new approach in blood vessel
labeling, and it has already been employed to create groundtruth labels for
widely accepted datasets such as DRIVE and STARE. The second groundtruth is
labeled by an expert ophthalmologist with more than 20 years of experience in
ROP screening. The expert refined the coarse groundtruth masks created by the
first observer to produce the second groundtruth. To achieve maximum accuracy,
FoV masks are also labeled manually. Since groundtruth labeling is time
consuming and requires expert involvement, we limit ourselves in this
experiment to 9 images which are used solely for testing purposes.

###  III-B Experimental protocol

We feed color retinal images from the ROP image dataset to the neural network
models proposed by Li et al. [9], Zhuang [10] and Oliveira et al. [4]. These
networks are considered to achieve state-of-the-art results in blood vessel
segmentation from retinal images in terms of accuracy, sensitivity and
specificity, as can be seen in Table I. For each network architecture we use
the latest pretrained model provided by authors. All pretrained models are
made open by authors and can be downloaded from GitHub repositories following
the work of Li et al., Oliveira et al. and Zhuang. Each network outputs a
probability map of the same size as an input image. Values in the probability
map range from 0 to 1 indicating probability that the corresponding pixel on
the input map is a part of the blood vessel. A probability map is then
thresholded with 0.5 threshold value to obtain a binary segmentation mask
where white pixels indicate that the corresponding input pixel is a part of
the blood vessel, whereas the black pixel indicates opposite. All experiments
are performed on a personal computer with 32GBs of RAM, Ryzen 3700X central
processing unit and NVIDIA RTX2060S graphics card.

###  III-C Performance Evaluation Metrics

To evaluate models performance, we adopt commonly used metrics in evaluation
of blood vessel segmentation solutions. These metrics include Accuracy (ACC),
True Positive Rate (TPR), True Negative Rate (TNR), F-measure ($F_{1}$), and
Area Under Curve (AUC) of a Receiver Operating Characteristics (ROC) graph.
The accuracy measures ratio of image pixels correctly classified as a blood
vessel or a background pixel. True Positive Rate and True Negative Rate, also
known as sensitivity and specificity, are more specific and measure ratio of
pixels correctly classified as a blood vessel or a background pixel,
respectively. F-measure is used as an integral measure of the model’s
precision and specificity.

##  IV Results and Discussion

Three models are compared in terms of accuracy, sensitivity, specificity,
F-measure and AUC on the images from our dataset as shown in Table I. All
values are calculated on predicted segmentation masks with applied FoV masks.
In Table I we also report numerical results from corresponding papers for
DRIVE, STARE, and CHASE_DB1 datasets to compare them to our results. We
evaluate models on two sets of groundtruth segmentation masks to evaluate the
models in predicting thicker, easier to see blood vessels annotated by non-
expert against small and hard to see blood vessels annotated by the expert.

TABLE I: Performance on DRIVE, STARE, CHASE_DB1 and our dataset (GT1 - first groundtruth, GT2 - second groundtruth) |  | Datasets
---|---|---
|  | ACC | TPR | TNR | F1 | AUC
Li et al. | DRIVE | 0.9574 | 0.7791 | 0.9831 | 0.8218 | 0.9816
STARE | 0.9782 | 0.7715 | 0.9919 | 0.8146 | 0.9915
CHASE | 0.9760 | 0.7969 | 0.9881 | 0.8073 | 0.9899
GT1 | 0.9678 | 0.4990 | 0.9881 | 0.5513 | 0.9378
GT2 | 0.9429 | 0.3457 | 0.9884 | 0.4307 | 0.8396
Oliveira et al. | DRIVE | 0.9576 | 0.8039 | 0.9804 | - | 0.9821
STARE | 0.9694 | 0.8315 | 0.9858 | - | 0.9905
CHASE | 0.9653 | 0.7779 | 0.9864 | - | 0.9855
GT1 | 0.9681 | 0.5271 | 0.9871 | 0.5685 | 0.9125
GT2 | 0.9436 | 0.3622 | 0.9875 | 0.4459 | 0.8250
Zhuang | DRIVE | 0.9561 | 0.7856 | 0.9810 | 0.8202 | 0.9793
STARE | - | - | - | - | -
CHASE | 0.9656 | 0.7978 | 0.9818 | 0.8031 | 0.9855
GT1 | 0.9639 | 0.5244 | 0.9830 | 0.5384 | 0.9307
GT2 | 0.9393 | 0.3633 | 0.9833 | 0.4291 | 0.8211

When compared to the first set of groundtruth values, all three models achieve
above 96% accuracy and above 98% specificity which is comparable to the
accuracy and specificity values reported on publicly available datasets.
However, models’ sensitivity significantly drops from approximately 77% for
public datasets to approximately 50% for our dataset, as the models
misclassified more blood vessel pixels as background pixels. Degradation in
model’s performance is also visible in AUC which is reduced up to 7.9%
indicating that models produce less reliable predictions on the ROP images.
The predicted segmentation masks are shown in Fig. 4. It can be seen that even
thick blood vessels are not segmented correctly if illumination issues are
present on the image which is often the case in infant retinal images due to
wide-field lenses used for screening.

Comparison with the second set of groundtruth values yields further
degradation in terms of segmentation accuracy, sensitivity and AUC, which drop
for approximately 2%, 15% and 15% compared to the performance calculated on
the first groundtruth dataset. The aggressive reduction in sensitivity and AUC
is caused by the inability of the existing models to segment tiny blood
vessels labeled by the expert as it can be seen in Fig. 4.

Existing solutions for retinal blood vessel segmentation also produce noise in
predicted segmentation maps, due to inability to differ between choroidal and
retinal blood vessels. Despite being blood vessels, choroidal blood vessels
are not retinal blood vessels and do not play part in ROP diagnosis.
Segmentation of these vessels alongside retinal vessels would aggravate ROP
diagnosis based on segmented vessels instead of alleviating it. Since the
chosen models perform outstandingly on public datasets, it is our educated
assumption that issues with choroidal blood vessels are caused by the lack of
images showing these vessels in public datasets used for network training.

Figure 4: Predicted segmentation maps. From top to bottom rows are shown: (1)
original images, (2) expert’s groundtruth and segmentation maps produced by
the methods of (3) Li et al., (4) Oliveira et al., and (5) Zhuang.

##  V Conclusions

In this paper we evaluated three state-of-the-art solutions for blood vessel
segmentation from color retinal images in context of retinopathy of
prematurity. The chosen solutions perform outstandingly when tested on
publicly available datasets. In this study we created a dataset consisting
solely of ROP images to test how well the solutions perform on ROP images and
can they be used as-is for blood vessel segmentation from ROP images.
Numerical measurements and visual results show that all three solutions
perform similarly and that all three have significant issues with uneven
illumination, infant’s thin blood vessels and choroidal blood vessels present
on images. There are two possible reasons for segmentation performance
degradation: (1) the images from the few public datasets that are used for
training do not have properties needed for effective retinal blood vessel
segmentation from ROP images and (2) trained models are not complex enough to
successfully address issues present on ROP images. Since the models tend to
under classify not only thin vessels, but in some cases also thick vessels
which is not the problem observed on the public datasets, we believe that not
all appearing issues can be solved by creating more complex CNN architecture
and that retinal blood vessel segmentation CNNs would benefit from a new, ROP-
aware retinal image dataset that would be used in CNN training. Additionally,
we observe that existing models have learned to classify blood vessels, but
not to distinguish between retinal and choroidal blood vessels which is
possibly due to the lack of images showing choroidal blood vessels in public
datasets.

As future work, more effort should be directed towards creating a more
adequate dataset that could be used to train efficient CNNs for blood vessel
segmentation from ROP retinal images. Alternatively, if ROP images are too
complex for proposed CNN architectures, more complex architecture could be
used to try to improve segmentation performance on existing datasets.

## References

  * [1] I. C. for the Classification of Retinopathy of Prematurity _et al._ , “The international classification of retinopathy of prematurity revisited.” _Archives of ophthalmology (Chicago, Ill.: 1960)_ , vol. 123, no. 7, p. 991, 2005.
  * [2] E. Poletti, D. Fiorin, E. Grisan, and A. Ruggeri, “Automatic vessel segmentation in wide-field retina images of infants with retinopathy of prematurity,” in _2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society_. IEEE, 2011, pp. 3954–3957.
  * [3] Q. Jin, Z. Meng, T. D. Pham, Q. Chen, L. Wei, and R. Su, “Dunet: A deformable network for retinal vessel segmentation,” _Knowledge-Based Systems_ , vol. 178, pp. 149–162, 2019.
  * [4] A. Oliveira, S. Pereira, and C. A. Silva, “Retinal vessel segmentation based on fully convolutional neural networks,” _Expert Systems with Applications_ , vol. 112, pp. 229–242, 2018.
  * [5] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [6] J. Son, S. J. Park, and K.-H. Jung, “Retinal vessel segmentation in fundoscopic images with generative adversarial networks,” _arXiv preprint arXiv:1706.09318_ , 2017.
  * [7] S. Y. Shin, S. Lee, I. D. Yun, and K. M. Lee, “Deep vessel segmentation by learning graphical connectivity,” _Medical image analysis_ , vol. 58, p. 101556, 2019.
  * [8] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.
  * [9] L. Li, M. Verma, Y. Nakashima, H. Nagahara, and R. Kawasaki, “Iternet: Retinal image segmentation utilizing structural redundancy in vessel networks,” in _The IEEE Winter Conference on Applications of Computer Vision_ , 2020, pp. 3656–3665.
  * [10] J. Zhuang, “Laddernet: Multi-path networks based on u-net for medical image segmentation,” _arXiv preprint arXiv:1810.07810_ , 2018.
  * [11] J. Staal, M. D. Abràmoff, M. Niemeijer, M. A. Viergever, and B. Van Ginneken, “Ridge-based vessel segmentation in color images of the retina,” _IEEE transactions on medical imaging_ , vol. 23, no. 4, pp. 501–509, 2004.
  * [12] A. Hoover, V. Kouznetsova, and M. Goldbaum, “Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response,” _IEEE Transactions on Medical imaging_ , vol. 19, no. 3, pp. 203–210, 2000\.
  * [13] C. G. Owen, A. R. Rudnicka, R. Mullen, S. A. Barman, D. Monekosso, P. H. Whincup, J. Ng, and C. Paterson, “Measuring retinal vessel tortuosity in 10-year-old children: validation of the computer-assisted image analysis of the retina (caiar) program,” _Investigative ophthalmology & visual science_, vol. 50, no. 5, pp. 2004–2010, 2009.
  * [14] A. Budai, R. Bock, A. Maier, J. Hornegger, and G. Michelson, “Robust vessel segmentation in fundus images,” _International journal of biomedical imaging_ , vol. 2013, 2013.
  * [15] J. Dai, H. Qi, Y. Xiong, Y. Li, G. Zhang, H. Hu, and Y. Wei, “Deformable convolutional networks,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 764–773.
  * [16] R. Estrada, C. Tomasi, S. C. Schmidler, and S. Farsiu, “Tree topology estimation,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 37, no. 8, pp. 1688–1701, 2014.
  * [17] H. Zhao, H. Li, S. Maurer-Stroh, and L. Cheng, “Synthesizing retinal and neuronal images with generative adversarial nets,” _Medical image analysis_ , vol. 49, pp. 14–26, 2018.