Efficient Column Generation for Cell Detection and
Segmentation

arXiv:1709.07337v1 [] 21 Sep 2017

Chong Zhanga,∗, Shaofei Wangb , Miguel A. Gonzalez-Ballestera,c , Julian
Yarkonyd,∗
a SimBioSys,

DTIC, Universitat Pompeu Fabra, Barcelona, Spain
b A&E Technologies, Beijing, China
c ICREA, Spain
d Experian Data Lab, USA

Abstract
We study the problem of instance segmentation in biological images with
crowded and compact cells. We formulate this task as an integer program where
variables correspond to cells and constraints enforce that cells do not overlap.
To solve this integer program, we propose a column generation formulation
where the pricing program is solved via exact optimization of very small scale
integer programs. Column generation is tightened using odd set inequalities
which fit elegantly into pricing problem optimization. Our column generation
approach achieves fast stable anytime inference for our instance segmentation
problems. We demonstrate on three distinct light microscopy datasets, with
several hundred cells each, that our proposed algorithm rapidly achieves or
exceeds state of the art accuracy.
Keywords: Combinatorial optimization, Column generation, Integer
programming, Large scale optimization, Linear programming

1. Introduction
Cell detection and instance segmentation are fundamental tasks for the study
of bioimages (Meijering, 2012) in the era of big data. Detection corresponds
∗ Corresponding

author
Email addresses: chong.zhang@upf.edu (Chong Zhang), julian.e.yarkony@gmail.com
(Julian Yarkony)

Preprint submitted to European Journal of Operational Research

September 22, 2017

to the problem of identifying individual cells and instance segmentation corresponds to the problem of determining the pixels corresponding to each of the
cells. Cells are often in close proximity and/or occlude each other. Traditionally bioimages were manually analyzed, however recent advances in microscope
techniques, automation, long-term high-throughput imaging, etc, result in vast
amounts of data from biological experiments, making manual analysis, and even
many computer aided methods with hand-tuned parameters, infeasible (Meijering et al., 2016; Hilsenbeck et al., 2017). The large diversity of cell lines and
microscopy imaging techniques require the development of algorithms for these
tasks to perform robustly and well across data sets.
In this paper we introduce a novel approach for instance segmentation specialized to bioimage analysis designed to rapidly produce high quality results
with little human intervention. The technique described in this paper is applicable to images that have crowded and compact cell regions acquired from
different modalities and cell shapes, as long as they produce intensity changes
at cell boundaries. Such patterns result from several microscopy imaging techniques, such as trans-illumination (e.g. bright field, dark field, phase contrast)
and fluorescence (e.g. through membrane or cytoplasmic staining) images. Thus
it is specifically suitable for images from which cells are almost transparent.
We formulate instance segmentation as the problem of selecting a set of
visually meaningful cells under the hard constraint that no two cells overlap
(share a common pixel). This problem corresponds to the classic integer linear
programming (ILP) formulation of the set packing problem (Karp, 1972) where
sets correspond to cells and elements correspond to pixels. The number of
possible cells is very large and can not be easily enumerated. We employ a
column generation approach (Barnhart et al., 1996) for solving the combinatorial
problem where the pricing problem is solved via exact optimization of very
small scale integer programs (IPs). Inference is made tractable by relying on
the assumption that cells are small and compact. When needed we tighten
the linear programming (LP) relaxation using odd set inequalities (Heismann
and Borndörfer, 2014). The use of odd set inequalities in our context does
2

not destroy the structure of the pricing problem so branch and price (Barnhart
et al., 1996) is not needed.
For the purpose of dimensionality reduction we employ the common technique of aggregating pixels into superpixels. Superpixels (Levinshtein et al.,
2009; Achanta et al., 2012) are the output of a dimensionality reduction technique that groups pixels in a close proximity with similar visual characteristics
and is commonly used as a preprocessing for image segmentation(Arbelaez et al.,
2011). Superpixels provide a gross over-segmentation of the image meaning
that they capture many boundaries not in the ground truth but miss very few
boundaries that are part of the ground truth. Hence we apply our set packing
formulation on the superpixels meaning that each set corresponds to a subset
of the superpixels and each element is a superpixel.
Our contributions consist of the following:
• Novel formulation of cell instance segmentation amenable to the tools and
methodology of the operations research community
• Structuring our formulation to admit tightening the corresponding LP
relaxation outside of branch/price methods
• Achieve benchmark level results on real microscopy datasets
We structure this document as follows. In Section 2 we consider the related
work in the fields of bioimage analysis and operations research. Next in Section 3
we introduce our set packing formulation of instance segmentation and our
column generation formulation with its corresponding pricing problem. Next
in Sections 4 and 5 we consider the production of anytime integral solutions
and lower bounds respectively. In Section 6 we demonstrate the applicability
of our approach to real bioimage data sets. Finally we conclude and consider
extensions in Section 7.

3

2. Related Work
2.1. Optimization in computer vision and bioimage analysis
Our work should be considered in the context of methods that are based on
cell boundary information and clustering of super-pixels. Relevant methodologies include contour profile pattern (Kvarnström et al., 2008; Mayer et al., 2013;
Dimopoulos et al., 2014), constrained label cost model (Zhang et al., 2014a),
correlation clustering (Zhang et al., 2014b; Yarkony et al., 2015), structured
learning (Arteta et al., 2012; Liu et al., 2014; Funke et al., 2015), and deep
learning (Ronneberger et al., 2015), etc. A comprehensive review can be found
in (Xing and Yang, 2016). Here we discuss the most relevant work.
The method of (Zhang et al., 2014b) frames instance segmentation as correlation clustering on a planar or nearly planar graph and relies heavily on
the planarity of their clustering problem’s structure in order to achieve efficient inference. Our work differs from (Zhang et al., 2014b) primarily from the
perspective of optimization. Notably, our model is not bound by planarity restrictions and instead relies on the assumption that cells are typically small and
compact. Therefore, our model is also applicable to 3D segmentation.
In (Zhang et al., 2015) the authors use depth to transform instance segmentation into a labeling problem and thus break the difficult symmetries found in
instance segmentation. They formulate the optimization as an ILP and solve it
using greedy network flow methods (Boykov et al., 2001), notably the Quadratic
Pseudo-Boolean optimization (QPBO) (Boros and Hammer, 2002; Rother et al.,
2007). The approach in (Zhang et al., 2015) also requires prior knowledge of the
number of labels present in the image, which is not realistic for images crowded
with hundreds or thousands of cells. In contrast our proposed ILP framework
does not require knowledge of the number of objects in the image.
Our inference approach is inspired by (Wang et al., 2017a) which tackles
multi-object tracking using column generation where the corresponding pricing
problem is solved using dynamic programming. This echoes the much earlier
operations research work in diverse areas such as vehicle routing (Ropke and

4

Cordeau, 2009), and cutting stock (Gilmore and Gomory, 1961) which use dynamic programming for pricing. In contrast our pricing problem optimization
is solved by many small ILPs, which can be run in parallel which echoes the
work in the operations research community of (Barahona and Jensen, 1998).
2.2. Column generation in operations research
Column generation (Gilmore and Gomory, 1961; Desaulniers et al., 2006;
Barnhart et al., 1996) is a popular approach for solving ILPs in which compact
formulations result in loose LP relaxations where the fractional solution tends
to be uninformative. Here uninformative means that the fractional solution can
not easily be rounded to a low cost integer solution. Column generation replaces
the LP with a new LP over a much larger space of variables which corresponds
to a tighter LP relaxation of the ILP (Geoffrion, 2010; Armacost et al., 2002).
The new LP retains the property from the original LP that it has a finite number
of constraints.
To solve the new LP, the dual of the new LP is considered which has a finite
number of variables and a huge number of constraints. Optimization considers
only a limited subset of the primal variables, which is initialized as empty, or
set heuristically. Optimization alternates between solving the LP relaxation
over the limited subset of the primal variables (called the master problem) and
identifying variables that correspond to violated dual constraints (which is called
pricing). Pricing often corresponds to combinatorial optimization which is often
an elegant dynamic program which has the powerful feature that many primal
variables (violated dual constraints) are generated at once. Approaches with
dynamic programming based pricing include (but are not limited to) the diverse
fields of cutting stock (Gilmore and Gomory, 1961, 1965), routing crews (Lavoie
et al., 1988; Vance et al., 1997), and routing vehicles (Ropke and Cordeau, 2009),
Column generation formulations can be tightened using branch-price methods (Barnhart et al., 2000, 1996; Vance, 1998) which is a variant of branch and
bound (Land and Doig, 1960) that is structured as to not disrupt the structure
of the pricing problem.
5

Term

Form

Index

Meaning

D

set

d

set of super-pixels

Q

set

q

set of cells

Q

{0, 1}|D|×|Q|

d, q

Qdq = 1 indicates that d in cell q

Γ

R|Q|

q

Γq is the cost of cell q

C

set

c

set of triples

γ

{0, 1}|Q|

q

γq = 1 indicates that cell q is selected.

θ

R|D|

d

θd is the cost of including d in a cell

ω

R

none

ω is the cost of instancing a cell

φ

R|D|×|D|

d1 , d2

φd1 d2 is the cost of including d1 , d2 in the same cell

V

R0+

S

|D|

d

Vd is the volume of super-pixel d

R0+

d1 , d2

Sd1 d2 is the distance between the centers of super-pixels d1 , d2

mV

R+

none

maximum volume of a cell

mR

R+

none

maximum radius of a cell

Q̂

set

q

set of cells generated during column generation

Cˆ

set

c

set of triples generated during column generation

Q̇

set

q

set of cells generated during a given iteration of column generation

C˙

set

c

set of triples generated during a given iteration of column generation

λ

R0+

d

Lagrange multipliers corresponding to super-pixels

κ

R0+

c

Lagrange multipliers corresponding to triples

x

{0, 1}|D|

d

xd = 1 indicates that super-pixel d is included in the column being generated

|D|×|D|

|D|
|C|

Table 1: Summary of Notation

Column generation has had few applications in computer vision until recently but has included diverse variants of correlation clustering (Yarkony and
Fowlkes, 2015; Yarkony et al., 2012; Yarkony, 2015) with applications to image
partitioning, multi-object tracking (Wang et al., 2017a), and multi-human pose
estimation (Wang et al., 2017b). Column generation in (Yarkony and Fowlkes,
2015; Yarkony et al., 2012; Yarkony, 2015) is notable in that the pricing problem is solved using the max cut on a planar graph (Shih et al., 1990; Barahona,
1982, 1991; Barahona and Mahjoub., 1986) which is known to be polynomial
time solvable via a reduction to perfect matching (Fisher, 1966).

3. Problem formulation
We now discuss our approach in detail. Given an image we start with computing a set of super-pixels (generally named super-voxels in 3D), which provides

6

an over-segmentation of cells. These super-pixels are then clustered into “perceptually meaningful” regions by constructing an optimization problem that
either groups the super-pixels into small coherent cells or labels them as background. The solution to this optimization problem corresponds to the globally
optimal selection of cells according to our model, which we formulate/solve as an
ILP. We consider our model below and summarize the corresponding notation
in Table 1.
Definitions. Let D be the set of super-pixels in an image, Q be the set of
all possible cells, and G ∈ {0, 1}|D|×|Q| be the super-pixel/cell incidence matrix
|D|×|D|

where Gdq =1 if and only if super-pixel d is part of the cell q. We use S ∈ R0+

to describe the Euclidean distance between super-pixels; where Sd1 d2 indicates
the distance between the centers of the super-pixel pair d1 and d2 . We use
|D|

V ∈ R+ to describe the area of super-pixels, with Vd being the area of superpixel d. The indicator vector γ ∈ {0, 1}|Q| gives a feasible segmentation solution,
where γq =1 indicates that cell q is included in the solution and γq =0 otherwise.
A collection of cells specified by γ is a valid solution if and only if each superpixel is associated with at most one active cell.
We use Γ ∈ R|Q| to define a cost vector, where Γq is the cost associated with
including cell q in the segmentation. Here we model such a cost with terms
θ ∈ R|D| and φ ∈ R|D|×|D| which are indexed by d and d1 , d2 respectively. We
use θd to denote the cost for including d in a cell and φd1 d2 to denote the cost
for including d1 and d2 in the same cell. We use ω ∈ R to denote the cost of
instancing a cell. We now define Γq in terms of θ, φ and ω.
Γq = ω +

X
d∈D

X

θd Gdq +

φd1 d2 Gd1 q Gd2 q ,

(1)

d1 ,d2 ∈D

Constraints. For most biological problems, it is valid to model cells of a given
type as having a maximum radius mR and a maximum area (volume if in 3D)
mV . Clearly, mV and mR are model defined parameters that vary from one
application to another, but they are also often known a-priori. The radius

7

Term

Effect of Positive Offset

θ

Decrease total volume of cells

φ

Fewer pairs of super-pixels in a common cell

ω

Decrease number of cells detected

mV

Increase maximum volume of cells

mR

Increase maximum radius of cells

Table 2: Summary of effect of offsetting values by positive offset

constraint can be written as follows:
∃[d∗ ; Gd∗ q = 1]

X

s.t. 0 =

Gd2 q [Sd∗ ,d2 > mR ] ∀q ∈ Q.

(2)

d2 ∈D

For any given q ∈ Q, any argument d∗ ∈ D satisfying Eq 2 is called an anchor
of q. Similarly, we write the area constraint as follows.
X
mV ≥
Gdq Vd ∀q ∈ Q.

(3)

d∈D

ILP formulation. Given the above variable definitions, we frame instance
segmentation as an ILP that minimizes the total cost of the selected cells:
X
min
Γq γq = min Γ> γ
(4)
P

γq ∈{0,1}∀q∈Q
Gdq γq ≤1 ∀d∈D q∈Q

q∈Q

γ∈{0,1}|Q|
Gγ≤1

The effect of our modeling parameters is summarized in Table 2.
Primal and Dual formulation. The LP relaxation of Eq 4 only contains
constraints for cells that share a common super-pixel. This generally results in
a tight relaxation, although not always. We tighten the relaxation using odd
set inequalities (Heismann and Borndörfer, 2014). Specifically we use odd set
inequalities of size three, (called triples), as similarly imposed in (Wang et al.,
2017a).
Triples are defined as follows: for any set of three unique super-pixels (called
a triple) the number of selected cells of Q that include two or more of superpixels in {d1 , d2 , d3 } can be no larger than one, i.e.
X
[Gd1 q + Gd2 q + Gd3 q ≥ 2]γq ≤ 1.
q∈Q

8

(5)

We denote the set of triples as C and describe it by a constraint matrix C∈{0, 1}|C|×|Q| ,
where Ccq = 1 if and only if cell q contains two or more members of set c. The
P
constraint matrix has a row for each triple: Ccq =[ d∈c Gdq ≥2], ∀c∈C, q∈Q.
The primal and dual LP relaxations of instance segmentation with constraints
on inequalities corresponding to triples are written below. The dual is expressed
|D|

|C|

using Lagrange multipliers λ ∈ R0+ and κ ∈ R0+ .
min Γ> γ =

γ≥0
Qγ≤1
Cγ≤1

max

λ≥0
κ≥0
>
Γ+Q λ+C > κ≥0

1> λ + 1> κ.

(6)

3.1. Algorithm
Since Q, C are intractably large, we use cutting plane method in the primal
and dual to build a sufficient subsets of Q,C.
We denote the nascent subsets of Q, C as Q̂, Cˆ respectively. In Alg 1 we write
column generation algorithm. We define the cutting plane/column generation
in Sections 3.2 and 3.3 respectively and display optimization in Fig 1. We use
Q̇, C˙ to refer to the columns and rows generated during a given iteration of our
algorithm.

9

Algorithm 1 Dual Optimization
Q̂ ← {}
Cˆ ← {}
repeat
Q̇ ← {}
C˙ ← {}
[λ, κ, γ] ← Solve Primal and Dual of Eq 6 over Q̂, Cˆ
for d∗ ∈ D do
q∗ ← arg min
0=

P

d∈c

P

d2

q∈Q
Qd∗ q =1
∈D Gd2 q [Sd∗ ,d2 >mR ]

Γq +

P

d∈D

Qdq λd +

P

c∈C

κc [2 ≤

Qdq ]
if Γq∗ +

P

d∈D

Qdq∗ λd +

P

c∈C

κc [2 ≤

P

d∈c

Qdq ] < 0 then

Q̇ ← [Q̇ ∪ q∗ ]
end if
end for
P
c∗ ← maxc∈C q∈Q Ccq γq
P
if q∈Q Ccq γq > 1 then
C˙ ← c∗
end if
Q̂ ← [Q̂, Q̇]
ˆ C]
˙
Cˆ ← [C,
until Q̇ = [] and C˙ = []

3.2. Row generation
Finding the most violated row consists of the following optimization.
max
c∈C

X

Ccq γq

(7)

q∈Q

Enumerating C is unnecessary and we generate its rows as needed by considering
only c = {dc1 dc2 dc3 } such that for each of pair dci , dcj there exists an index q
such that γq > 0 and Qdi q = Qdj q = 1. Generating rows is done only when
no (significantly) violated columns exist. Triples are only added to Cˆ if the
10

corresponding constraint is violated. We can add one or more than one per
iteration depending on a schedule chosen.
3.3. Generating columns
Violated constraints in the dual correspond to primal variables (cells) that
may improve the primal objective. To identify such primal variables we compute for each d∗ ∈ D the most violated dual constraint corresponding to a cell
such that d∗ is an anchor of that cell. The corresponding cell is described using indicator vector x∈{0, 1}|D| , where the corresponding column is defined as
Qdq ← xd , ∀d ∈ D. We write the pricing problem as an IP below.
min

x∈{0,1}|Q|

X

(θd + λd )xd +

d∈D

X

φd1 d2 xd1 xd2 +

d1 ,d2 ∈D

X
c∈C

κc ([2 ≤

X

xd ])

d∈c

s.t. xd∗ = 1
xd = 0 ∀d ∈ D
X
V d xd ≤ m V

s.t. Sd,d∗ > mR
(8)

d∈D

For our data sets of images crowded with several hundreds of cells, the maximum radius of a cell is relatively small and the number of super-pixels within the
radius of a given super-pixel is of the order of tens and often around ten. Therefore, solving Eq 8 is efficient and can be done in parallel for each d∗ ∈ D. We
tackle this by converting Eq 8 to an ILP and then solving it with an off-the-shelf
ILP solver.

4. Anytime Integral Solutions
We now consider the anytime production of integral solutions in the master
problem. While set packing NP hard in general, in practice the LP relaxations
are integral at termination and generally integral after each step of optimization. For cases where the LP is loose, we find that solving the ILP given the
primal variables generated takes little additional time beyond solving the LP.
However we can use rounding procedures (Wang et al., 2017a) when difficult

11

Algorithm 2 Upper Bound Rounding
while ∃q ∈ Q s.t. γq ∈
/ {0, 1} do
P
∗
q ← arg min q∈Q Γq γq − q̂∈Q⊥q Γq̂ γq̂
1>γq >0

γq̂ ← 0

∀q̂ ∈ Q⊥q∗

γq ∗ ← 1
end while
RETURN γ

ILPs occur. Specifically, we tackle the rounding of a fractional γ with a greedy
iterative approach. At each iteration, it selects the cell q with non-binary γq
that minimizes Γq γq discounted by the fractional cost of any cells that share
a super-pixel with q; hence can no longer be added to the segmentation if q is
already added. We write the rounding procedure in Alg 2 using the notation
Q⊥q to indicate the set of cells in Q that intersect cell q (excluding q itself).
5. Lower bounds
We now consider the production of anytime lower bounds on the optimal
integral solution. We first write the ILP for cell instance segmentation and then
introduce Lagrange multipliers.
min

γ∈{0,1}|Q|
Qγ≤1
Cγ≤1

Γ> γ =

min

max Γ> γ + (−λ> 1 + λ> Qγ) + (−κ> 1 + κ> Cγ) (9)

γ∈{0,1}|Q| λ≥0
κ≥0
Qγ≤1

We now relax the constraint in Eq 9 that the dual variables are optimal producing the following lower bound.
Eq 9 ≥

min

γ∈{0,1}|Q|
Qγ≤1

Γ> γ + (−λ> 1 + λ> Qγ) − (κ> 1 + κ> Cγ)

= −κ> 1 − λ> 1 +

min

γ∈{0,1}|Q|
Qγ≤1

(Γ + Q> λ + C > κ)> γ

(10)

Recall that every cell is associated with at least one anchor. We denote the set
of anchors associated with a given cell q as Nq . We use Q:,q , C:,q to refer to the
12

column q of the matrices Q, C respectively. Given any fixed γ ∈ {0, 1}|Q| such
that Qγ ≤ 1 observe the following.
(Γ + Q> λ + C > κ)> γ ≥

X

>
min[0, min γq (Γq + Q>
:,q λ + C:,q κ)]
q∈Q
d∈Nq

d∈D

(11)

We now use Eq 11 to produce following lower bound on Eq 10.
Eq 10 ≥

−κ> 1 − λ> 1 +

X

min

γ∈{0,1}|Q|
d∈D
Qγ≤1

>
min[0, min γq (Γq +Q>
:,q λ+C:,q κ)] (12)
q∈Q
d∈Nq

We now relax the constraint in Eq 12 that Qγ ≤ 1 producing the following lower
bound.
Eq 12 ≥ −κ> 1−λ> 1+
= −κ> 1 − λ> 1 +

X

min

X

γ∈{0,1}|Q|

d∈D

>
min[0, min γq (Γq +Q>
:,q λ+C:,q κ)]
q∈Q
d∈Nq

>
min[0, min (Γq + Q>
:,q λ + C:,q κ)]

d∈D

q∈Q
d∈Nq

(13)

>
Observe that the term min q∈Q (Γq + Q>
:,q λ + C:,q κ) is identical to the optimizad∈Nq

tion computed at every stage of column generation.

6. Results
The technique described in this paper is applicable to images crowded with
cells which are mainly discernible by boundary cues. Such images can be acquired from different modalities and cell types. Here we evaluate our algorithm
on three datasets. Challenges of these datasets include: densely packed and
touching cells, out-of-focus artifacts, variations on shape and size, changing
boundaries even on the same cell, as well as other structures showing similar
boundaries.
6.1. Experiment settings
To ensure detecting cell boundaries with varying patterns, a trainable classifier seems to be the right choice. For each dataset, we choose to train a Random
Forest (RF) classifier from the open source software, ilastik (Sommer et al.,
13

2011), to discriminate: (1) boundaries of in-focus cells; (2) in-focus cells; (3)
out-of-focus cells; and (4) background. The posterior probabilities for class (1)
is used as the pairwise potentials. For training, we used < 1% pixels per dataset
with generic features e.g. Gaussian, Laplacian, Structured tensor. Subsequent
steps use the posterior probabilities to calculate parameters and require no more
training. The prediction from the class boundaries of in-focus cells is also used
to generate super-pixels. And those for classes (3) and (4) are combined and
inverted to create a foreground prediction. Here foreground corresponds to the
superpixels that are part of cells which are background otherwise. For each
super-pixel, the proportion of its foreground part defines the unary potential θ
which we then offset by a constant fixed for each dataset. A summary about
the parameters used in our experiments are shown in Table 3.
6.2. Evaluation
A visualization of the results can be seen in Fig 2. Quantitatively, we compare the performance of our algorithm with those reported in the state-of-the-art
methods (Arteta et al., 2012, 2016; Funke et al., 2015; Hilsenbeck et al., 2017;
Dimopoulos et al., 2014; Ronneberger et al., 2015; Zhang et al., 2014b), in terms
of detection (precision, recall and F-score) and segmentation (Dice coefficient
and Jaccard index). For detection, we establish possible matches between found
regions and ground truth (GT) regions based on overlap, and find a Hungarian matching using the centroid distance as minimizer. Unmatched GT regions
are FN, unmatched segmentation regions are FP. Jaccard index is computed
between the area of true positive (TP) detection regions Rtpd and the area of
GT region Rgt : (Rtpd ∩Rgt )/(Rtpd ∪Rgt ). They are summarized in Tables 4
and 5. In general, our method achieves or exceed state of the art performance.
Additionally, our method requires very little training for the RF classifiers, as
opposed to methods like (Arteta et al., 2012, 2016; Funke et al., 2015), which
require fully labeled data for training. This is an advantage of relieving human
annotations when several hundreds of cells need to be labeled per image. Also,
our method can handle very well large variations of cell shape/size even in the
14

Cell image
Compute
θ, φ

Merge
λ, κ, θ

Column

Column

Column

Column

(λ, κ, 1)

(λ, κ, 2)

(λ, κ, d)

(λ, κ, D)

Cˆ

Q̂

Q̂
C˙
Cˆ

λ, κ

LP Solver

γ

Row(γ)

γ
Round LP
γ
Cell
Pipeline

Figure 1: Overview of our system. We use colors to distinguish between the different parts
of the system, which are defined as follows: user input (green), sub-problem solution (blue),

15 (orange). We use Column(λ, κ, d) to refer
triples (red), rounding the output of the LP solver
to generating the column where d is an anchor.

Figure 2: Example cell segmentation results of Datasets 1-3 (left to right). Rows are (top
to bottom): original image, cell of interest boundary classifier prediction image, super-pixels,
color map of segmentation, and enlarged views of the inset (black square)

16

Table 3:

Summary of experimental datasets on the number of cells, cell radius, image size, the number of super-pixels and region adjacent graph

(RAG) edges.

Dataset

# cells

avg. cell radius

image size

# super-pixels

# RAG edges

1 (Zhang et al., 2014a)

1768

30

1024×1024

1225±242

3456±701

2 (Peng et al., 2013)

2340

50

1024×1024

3727±2450

10530±7010

3 (Arteta et al., 2012)

1073

20

400×400

1081±364

3035±1038

Table 4: Evaluation and comparison of detection for Datasets 1-3 (Fig. 2) on precision (P), recall (R), F-score (F), dice coefficient (D) and Jaccard
index (J) are reported for the proposed method, as well as those reported in the state-of-the-art methods. Here (Zhang et al., 2014b) uses the
algorithms planar correlation clustering (PCC) and non-planar correlation clustering (NPCC).

Dataset 1

Dataset 2

Dataset 3

18

P

R

F

P

R

F

P

R

F

(Arteta et al., 2012)

-

-

-

-

-

-

0.89

0.86

0.87

(Arteta et al., 2016)

-

-

-

-

-

-

0.99

0.96

0.97

(Funke et al., 2015)

0.93

0.89

0.91

0.99

0.90

0.94

0.95

0.98

0.97

(Hilsenbeck et al., 2017)

-

-

-

-

-

-

-

-

0.97

(Ronneberger et al., 2015)

-

-

-

-

-

-

-

-

0.97

PCC (Zhang et al., 2014b)

0.95

0.86

0.90

0.80

0.75

0.76

0.92

0.92

0.92

NPCC (Zhang et al., 2014b)

0.71

0.96

0.82

0.75

0.83

0.78

0.85

0.97

0.90

Proposed

0.99

0.97

0.98

1.00

0.94

0.97

1.00

0.97

0.99

Table 5:

Evaluation and comparison of segmentation for Datasets 1-3 (Fig. 2) on dice coefficient (D) and Jaccard index (J) are reported for

the proposed method, as well as those reported in the state-of-the-art methods. Here (Zhang et al., 2014b) uses the algorithms planar correlation
clustering (PCC) and non-planar correlation clustering (NPCC).

19

Dataset 1

Dataset 2

Dataset 3

D

J

D

J

D

J

(Funke et al., 2015)

0.90

0.82

0.90

0.83

0.84

0.73

(Hilsenbeck et al., 2017)

-

-

-

-

-

0.75

(Dimopoulos et al., 2014)

-

0.87

-

-

-

-

(Ronneberger et al., 2015)

-

-

-

-

-

0.74

PCC (Zhang et al., 2014b)

0.87

0.84

0.91

0.85

0.79

0.72

NPCC (Zhang et al., 2014b)

0.86

0.89

0.91

0.84

0.80

0.70

Proposed

0.91

0.90

0.90

0.83

0.82

0.71

Histogram of Inference Time: Cells

150

total time
time without upper bound

counts

100

50

0
0

100

200

300

400

500

600

700

time

Figure 3: Histogram of inference time for Dataset 1.

same image, as shown in Fig 2 for Dataset 2.
6.3. Timing and bounds
We now consider the performance of our approach with regard to the gap
between the upper and lower bounds produced by our algorithm. We normalize these gaps by dividing by the absolute value of the lower bound. For our
three data set the proportion of problem instances that achieve normalized gaps
under 0.1 are 99.28 %, 80 % and 100 %, on Datasets 1,2,3 respectively. The
peak histogram of inference time are around 150, 500 and 100 seconds without
parallelization. As an example, the histogram of inference time for Dataset 1
is shown in Fig 3. Our approach is approximately an order of magnitude faster
than that of (Zhang et al., 2014b).

7. Conclusion
In this article we introduce a novel column generation strategy that efficiently
optimizes an ILP formulation of instance segmentation through clustering superpixels. We use our approach to detect and segment crowded clusters of cells in

20

distinct microscopy image datasets and achieves state of the art or near state
of the art performance.
We now consider some extensions of our approach. The use of odd set
inequalities may prove useful for traditional set cover formulations of vehicle
routing problems. In this context for triples the corresponding inequality is
defined as follows: For any set of three unique depots the number of routes that
pass through one or two or those depots plus two times the number of routes
that pass through all three depots is no less than two. Dynamic programming
formulations for pricing can be adapted to include the corresponding Lagrange
multipliers (Irnich and Desaulniers, 2005) ( in either the elementary or nonelementary (Kallehauge et al., 2005) setting). The approach in (Wang et al.,
2017a) can also be adapted which employs dynamic programming in a branch
and bound context in the pricing problem (never the master problem).
Another extension considers multiple types of cells with a unique model for
each cell type and its own pricing problem. Such types can include rotations,
scalings, or other transformations of a common model which may be useful for
cells that highly non-circular in shape.
In future work one should apply dual feasible inequalities (Ben Amor et al.,
2006; Yarkony and Fowlkes, 2015). In this case one would create separate variable for each pair of cell, feasible anchor for that cell (where the anchor is called
the the main anchor). Then the ILP would be framed as selecting a set of cells
such that (1) no super-pixel is included more than once and (2) no main anchor
is included in more than one cell. However the Lagrange multipliers for (1) can
be bounded from above by the increase in cost corresponding to removing the
super-pixel d from a cell. For a super-pixel d one trivial such bound is minus
one times the sum of the non-positive cost terms involving d.

Acknowledgment
This work was partly supported by the Spanish Ministry of Economy and
Competitiveness under the Maria de Maeztu Units of Excellence Programme
(MDM-2015-0502).
21

References
Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., Ssstrunk, S., 2012.
SLIC superpixels compared to state-of-the-art superpixel methods. IEEE
Transactions on Pattern Analysis and Machine Intelligence 34, 2274–2282.
Arbelaez, P., Maire, M., Fowlkes, C., Malik, J., 2011. Contour detection and
hierarchical image segmentation. IEEE Transactions on Pattern Analysis and
Machine Intelligence 33, 898–916.
Armacost, A.P., Barnhart, C., Ware, K.A., 2002. Composite variable formulations for express shipment service network design. Transportation Science 36,
1–20.
Arteta, C., Lempitsky, V., Noble, J., Zisserman, A., 2012. Learning to detect
cells using non-overlapping extremal regions, in: International Conference on
Medical Image Computing and Computer Assisted Intervention (MICCAI).
volume 7510 of Lecture Notes in Computer Science, pp. 348–356.
Arteta, C., Lempitsky, V., Noble, J., Zisserman, A., 2016. Detecting overlapping
instances in microscopy images using extremal region trees. Medical Image
Analysis 27, 3–16.
Barahona, F., 1982. On the computational complexity of ising spin glass models.
Journal of Physics A: Mathematical, Nuclear and General 15, 3241–3253.
Barahona, F., 1991. On cuts and matchings in planar graphs. Mathematical
Programming 36, 53–68.
Barahona, F., Jensen, D., 1998. Plant location with minimum inventory. Mathematical Programming 83, 101–111.
Barahona, F., Mahjoub., A., 1986. On the cut polytope. Mathematical Programming 60, 157–173.

22

Barnhart, C., Hane, C.A., Vance, P.H., 2000. Using branch-and-price-and-cut
to solve origin-destination integer multicommodity flow problems. Operations
Research 48, 318–326.
Barnhart, C., Johnson, E.L., Nemhauser, G.L., Savelsbergh, M.W.P., Vance,
P.H., 1996. Branch-and-price: Column generation for solving huge integer
programs. Operations Research 46, 316–329.
Ben Amor, H., Desrosiers, J., Valério de Carvalho, J.M., 2006. Dual-optimal
inequalities for stabilized column generation. Operations Research 54, 454–
463.
Boros, E., Hammer, P., 2002. Pseudo-boolean optimization. Discrete Applied
Mathematics 123, 155–225.
Boykov, Y., Veksler, O., Zabih, R., 2001. Fast approximate energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and Machine
Intelligence 23, 1222–1239.
Desaulniers, G., Desrosiers, J., Solomon, M.M., 2006. Column Generation.
volume 5. Springer Science & Business Media.
Dimopoulos, S., Mayer, C., Rudolf, F., Stelling, J., 2014. Accurate cell segmentation in microscopy images using membrane patterns. Bioinformatics 30,
2644–2651.
Fisher, M.E., 1966. On the dimer solution of planar ising models. Journal of
Mathematical Physics 7, 1776–1781.
Funke, J., Hamprecht, F., Zhang, C., 2015. Learning to segment: Training
hierarchical segmentation under a topological loss, in: International Conference on Medical Image Computing and Computer Assisted Intervention
(MICCAI). volume 9351 of Lecture Notes in Computer Science, pp. 268–275.
Geoffrion, A.M., 2010. Lagrangian relaxation for integer programming, in: et al.,
M.J. (Ed.), 50 Years of Integer Programming 1958-2008. Springer. chapter 9,
pp. 243–281.
23

Gilmore, P., Gomory, R.E., 1965. Multistage cutting stock problems of two and
more dimensions. Operations research 13, 94–120.
Gilmore, P.C., Gomory, R.E., 1961. A linear programming approach to the
cutting-stock problem. Operations research 9, 849–859.
Heismann, O., Borndörfer, R., 2014. A generalization of odd set inequalities for
the set packing problem, in: Operations Research Proceedings 2013. Springer,
pp. 193–199.
Hilsenbeck, O., Schwarzfischer, M., Loeffler, D., Dimopoulos, S., Hastreiter, S.,
Marr, C., Theis, F., Schroeder, T., 2017. fastER : a user-friendly tool for ultrafast and robust cell segmentation in large-scale microscopy. Bioinformatics
33, 2020–2028.
Irnich, S., Desaulniers, G., 2005. Shortest path problems with resource constraints. Column Generation , 33–65.
Kallehauge, B., Larsen, J., Madsen, O.B., Solomon, M.M., 2005. Vehicle routing
problem with time windows. Column Generation , 67–98.
Karp, R.M., 1972. Reducibility among combinatorial problems, in: Complexity
of computer computations. Springer, pp. 85–103.
Kvarnström, M., Logg, K., Diez, A., Bodvard, K., Kall, M., 2008. Image analysis
algorithms for cell contour recognition in budding yeast. Optics Express 16,
1035–1042.
Land, A.H., Doig, A.G., 1960. An automatic method of solving discrete programming problems. Econometrica: Journal of the Econometric Society ,
497–520.
Lavoie, S., Minoux, M., Odier, E., 1988. A new approach for crew pairing
problems by column generation with an application to air transportation.
European Journal of Operational Research 35, 45–58.

24

Levinshtein, A., Stere, A., Kutulakos, K.N., Fleet, D.J., Dickinson, S.J., Siddiqi, K., 2009. Turbopixels: Fast superpixels using geometric flows. IEEE
Transactions on Pattern Analysis and Machine Intelligence 31, 2290–2297.
Liu, F., Xing, F., Yang, L., 2014. Robust muscle cell segmentation using region
selection with dynamic programming, in: IEEE International Symposium on
Biomedical Imaging (ISBI), pp. 1381–1384.
Mayer, C., Dimopoulos, S., Rudolf, F., Stelling, J., 2013. Using CellX to Quantify Intracellular Events. Current Protocols in Molecular Biology , 14.22.1–
14.22.20.
Meijering, E., 2012. Cell segmentation: 50 years down the road. IEEE Signal
Processing Magazine , 140–145.
Meijering, E., Carpenter, A., Peng, H., Hamprecht, F., Olivo-Marin, J.C., 2016.
Imagining the future of bioimage analysis. Nature Biotechnology 34, 1250–
1255.
Peng, J.Y., Chen, Y.J., Green, M.D., Sabatinos, S.A., Forsburg, S.L., Hsu, C.N.,
2013. PombeX: Robust cell segmentation for fission yeast transillumination
images. PLoS One 8, e81434.
Ronneberger, O., Fischer, P., Brox, T., 2015. U-Net: convolutional networks
for biomedical image segmentation, in: Frangi, A., et al. (Eds.), International
Conference on Medical Image Computing and Computer Assisted Intervention
(MICCAI). Springer. volume 9351 of Lecture Notes in Computer Science, pp.
234–241.
Ropke, S., Cordeau, J.F., 2009. Branch and cut and price for the pickup and
delivery problem with time windows. Transportation Science 43, 267–286.
Rother, C., Kolmogorov, V., Lempitsky, V., Szummer, M., 2007. Optimizing
binary mrfs via extended roof duality, in: 2007 IEEE Conference on Computer
Vision and Pattern Recognition (CVPR).

25

Shih, W.K., Wu, S., Kuo, Y., 1990. Unifying maximum cut and minimum cut
of a planar graph. IEEE Transactions on Computers 39, 694–697.
Sommer, C., Straehle, C., Koethe, U., Hamprecht, F.A., 2011. ilastik: Interactive Learning and Segmentation Toolkit, in: IEEE International Symposium
on Biomedical Imaging (ISBI).
Vance, P.H., 1998. Branch-and-price algorithms for the one-dimensional cutting
stock problem. Computational Optimization and Applications 9, 211–228.
Vance, P.H., Barnhart, C., Johnson, E.L., Nemhauser, G.L., 1997. Airline crew
scheduling: A new formulation and decomposition algorithm. Operations
Research 45, 188–200.
Wang, S., Wolf, S., Fowlkes, C., Yarkony, J., 2017a. Tracking objects with
higher order interactions using delayed column generation, in: International
Conference on Artificial Intelligence and Statistics (AISTATS).
Wang, S., Zhang, C., Gonzalez-Ballester, M.A., Ihler, A., Yarkony, J.,
2017b. Multi-person pose estimation via column generation. arXiv preprint
arXiv:1709.05982 .
Xing, F., Yang, L., 2016. Robust nucleus/cell detection and segmentation in
digital pathology and microscopy images: A comprehensive review. IEEE
Reviews in Biomedical Engineering 9, 234–263.
Yarkony, J., 2015. Next generation multicuts for semi-planar graphs. arXiv
preprint arXiv:1511.01994 .
Yarkony, J., Fowlkes, C., 2015. Planar ultrametrics for image segmentation, in:
Neural Information Processing Systems.
Yarkony, J., Ihler, A., Fowlkes, C., 2012. Fast planar correlation clustering for
image segmentation, in: Proceedings of the 12th European Conference on
Computer Vision (ECCV).

26

Yarkony, J., Zhang, C., Fowlkes, C., 2015. Hierarchical planar correlation clustering for cell segmentation, in: Energy Minimization Methods in Computer
Vision and Pattern Recognition (EMMCVPR), pp. 492–504.
Zhang, C., Huber, F., Knop, M., Hamprecht, F., 2014a. Yeast cell detection and
segmentation in bright field microscopy, in: IEEE International Symposium
on Biomedical Imaging (ISBI).
Zhang, C., Yarkony, J., Hamprecht, F., 2014b. Cell detection and segmentation
using correlation clustering, in: International Conference on Medical Image
Computing and Computer Assisted Intervention (MICCAI). volume 8673 of
Lecture Notes in Computer Science, pp. 9–16.
Zhang, Z., Schwing, A., Fidler, S., Urtasun, R., 2015. Monocular object instance
segmentation and depth ordering with CNNs, in: International Conference on
Computer Vision (ICCV).

27

