for all $p>0$ (cf. [H1], Chap. I, Theorem 4.5). Moreover, for $p=1$, it is
equivalent to the condition (b) of Proposition 12.4.

(b) It is likely that Theorem 12.7 (and hence Corollary 12.8) is valid for any
$\mathfrak{g}$ and any $\lambda,\mu\in{\mathscr{P}^{+}_{S}}^{o}.$

## References

  * [BJK] Besson, M., Jeralds, S., Kiers, J., Multiplicity in root components via geometric Satake, arXiv: 1909.05103.
  * [Bou] Bourbaki, N., Lie groups and Lie algebras, Chapters 4-6, Springer (2005).
  * [BrKu] Brown, M., Kumar, S., A study of saturated tensor cone for symmetrizable Kac–Moody algebras, Math. Ann. 360 (2014).
  * [H1] Hartshorne, R., On the De Rham cohomology of algebraic varieties, Publ. Math. IHES 45 (1976).
  * [H2] Hartshorne, R., Algebraic Geometry, Springer-Verlag (1977).
  * [Kac] Kac, V., Infinite dimensional Lie algebras, Cambridge University Press (1990).
  * [KRR] Kac, V., Raina, A., Rozhkovskaya, N., Bombay lectures on highest weight representations of infinite dimensional Lie algebras, Adv. Series in Mathematical Physics, vol. 29, World Scientific (2013).
  * [Kas1] Kashiwara, M., On crystal bases of the Q-analogue of universal enveloping algebras, Duke Math Journal 63 (1991).
  * [Kas2] Kashiwara, M., Global bases of quantum groups, Duke Math Journal 69 (1993).
  * [Kos] Kostant, B., A formula for the multiplicity of a weight, Trans. Am. Math. Soc. 93 (1959).
  * [Ku1] Kumar, S., Proof of Wahl’s conjecture on surjectivity of the Gaussian map for flag varieities, American Journal of Mathematics 114 (1992).
  * [Ku2] Kumar, S., Existence of certain components in the tensor product of two integrable highest weight modules for Kac–Moody algebras, in: “Infinite dimensional Lie algebras and groups” (ed. by V.G. Kac). Advanced series in Mathematical Physics Vol. 7, World Scientific (1989).
  * [Ku3] Kumar, S., Kac–Moody groups, their flag varieties, and representation theory, Birkhäuser (2002).
  * [Ku4] Kumar, S., Tensor product decomposition, Proceedings of the International Congress of Mathematician, Hyderabad, India (2010).
  * [M1] Mathieu, O., Formules de caractères pour les algèbres de Kac-Moody générales, Asterisque vol. 159-160 (1988).
  * [M2] Mathieu, O., Construction d’un groupe de Kac-Moody et applications, Compositio Math. 69 (1989).
  * [Wah] Wahl, J., Gaussian maps and tensor products of irreducible representations, Manuscripta Math. 73 (1991).


11institutetext:  Department of Computer Science, Bar-Ilan University, Ramat-
Gan 5290002, Israel
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Gator: Customizable Channel Pruning
of Neural Networks with Gating

Eli Passov  Eli O. David  Nathan S. Netanyahu Nathan Netanyahu is also
affiliated with the Department of Computer Science at the College of Law and
Business, Ramat-Gan 5257346, Israel

###### Abstract

The rise of neural network (NN) applications has prompted an increased
interest in compression, with a particular focus on channel pruning, which
does not require any additional hardware. Most pruning methods employ either
single-layer operations or global schemes to determine which channels to
remove followed by fine-tuning of the network. In this paper we present Gator,
a channel-pruning method which temporarily adds learned gating mechanisms for
pruning of individual channels, and which is trained with an additional
auxiliary loss, aimed at reducing the computational cost due to memory,
(theoretical) speedup (in terms of FLOPs), and practical, hardware-specific
speedup. Gator introduces a new formulation of dependencies between NN layers
which, in contrast to most previous methods, enables pruning of non-sequential
parts, such as layers on ResNet’s highway, and even removing entire ResNet
blocks. Gator’s pruning for ResNet-50 trained on ImageNet produces state-of-
the-art (SOTA) results, such as $50\%$ FLOPs reduction with only $0.4\%$-drop
in top-5 accuracy. Also, Gator outperforms previous pruning models, in terms
of GPU latency by running 1.4 times faster. Furthermore, Gator achieves
improved top-5 accuracy results, compared to MobileNetV2 and SqueezeNet, for
similar runtimes. The source code of this work is available at:
https://github.com/EliPassov/gator.

5in(0mm, -65mm) Ref: An earlier version of this arXiv appeared at the
International Conference on Artificial Neural Networks (ICANN), Springer LNCS,
Vol. 12894, pp. 46–58, Bratislava, Slovakia, September 2021.

##  1 Introduction

The use of NNs, in particular, convolutional neural networks (CNNs), in a
variety of Computer Vision applications has surged considerably in recent
years. Most NNs consist of tens to hundreds of millions of weight parameters
and require typically billions of floating point operations (FLOPs) for
processing a single image. This limits their use on available hardware,
especially on small devices. Consequently, many works in recent years have
focused on reducing the computational cost of network architectures without
severely impairing their performance, either by creating more efficient
network architectures or by applying network compression. Pruning is a popular
method among compression methods, where network weights and activations are
removed. Pruning methods can be categorized as unstructured, i.e., referring
to the removal of individual weights, or structured, i.e., referring mostly to
the removal of convolutional channels. Unstructured pruning is considered much
less practical, as it requires special hardware for efficient application.

More recent techniques also utilize dynamic computation in conjunction with
gating mechanisms to control the network operations that are invoked for each
input instance; thus, the computational cost is determined dynamically for
each individual input. Our method utilizes dynamic gating mechanisms only
during training, to learn which channels should be removed, thereby producing
a pruned network without any gating mechanisms.

The effectiveness of pruning methods is usually assessed by these parameters:

  * •

Reduction of memory required for storing weights;

  * •

Theoretical speedup, measured by the reduction of the number of FLOPs;

  * •

Latency/runtime on a specific hardware.

Note that latency is not commonly used as an optimization objective, and that
theoretical speedup is never translated into the exact same practical speedup.
For unstructured methods, this is due to the fact that most hardware does not
support sparse operations, leading to no optimization in practice; for
structured pruning, the sole focus on the FLOPs metric ignores hardware-
specific behavior.

This paper makes the following contributions:

  * •

Presentation of Gator, an innovative pruning technique that utilizes channel
hard-gating to remove channels during training while also adjusting the
remaining weights;

  * •

Introduction of a novel formulation for dependencies between NN layers
allowing Gator to prune non-sequential parts of NN architectures (such as
ResNet’s highway), which are not handled by most pruning methods;

  * •

Presentation of a customizable cost function, which reflects various practical
compression criteria, such as memory reduction, FLOPs reduction, and hardware-
specific speedup;

  * •

Obtaining new SOTA results for ResNet-50 pruning on ImageNet, in particular,
50% FLOPs reduction with $0.4\%$ top-5 accuracy drop and runtime speedup
(latency reduction) of 1.4 for the same accuracy as other methods;

  * •

Presentation of additional results for higher pruning rates, outperforming
MobileNet V2 and SqueezeNet for the same runtime.

##  2 Related Work

###  2.1 Network Compression

Most methods for NN compression can be categorized into three categories.
First, tensor decomposition ([20], [15], [17]), where the convolutional
operation is approximated by decomposing it to smaller operations.

Second, quantization methods ([33], [34], [2]), in particular, binary methods
([24], [3]), where the network weights are approximated with lower precision
variables allowing for memory reduction.

And third, pruning, where parts of layer outputs or weights are removed. In
early works [21] and more recent ones [8], pruning methods are characterized
by weight elimination and a ”fine tuning” training stage. A thorough review of
pruning methods can be found in [1], which also highlights the distinction
between structured and unstructured pruning. Since unstructured pruning cannot
reduce computational costs on most hardware types, we will focus and compare
Gator only to other structured pruning methods.

###  2.2 Structured Channel Pruning

Most channel pruning works differ in the criteria by which channels are
pruned. In [22] the smallest sums of weights of channels from all layers are
pruned, while non sequential layers are being considered. The pruning
criterion in [23] is the first order Taylor expansion of the loss
($w\frac{\partial L}{\partial w}$). Weights closest to the geometric median
are pruned in [9], in an iterative process for each layer. And in [32],
individually trained squeeze and excitation [12] blocks are used to prune
layers with the lower attention values.

Other methods, such as [11], remove the convolutional input channels with the
least impact on the output. While [35] first train auxiliary classification
losses on intermediate layers and define reconstruction losses on potential
channels to be remove. The channels chosen to be removed are those which have
the smallest impact on the sum total of all classification and reconstruction
losses.

###  2.3 Gating

Dynamic gating modules determine, for each individual input, which parts are
executed and which are skipped within a computational NN graph. In [30],
squeeze and excitation blocks [12] are used with a gumbel softmax [16] gate on
individual residual blocks, enabling skipping some blocks dynamically. In [4],
the same gating mechanism is applied on individual channels.

Methods which do not rely on gumbel softmax, such as [7], use regression on
the global pooling vector of the convolutional input to determine dynamically
which channels are passed to the convolution. In [13], convolutions are
divided into two parts, such that the output of the first gates the execution
of the second.

A soft sigmoid gate is used in [6] to prune connections in fully-connected
layers. Our method also uses a sigmoid; however, it gates entire channels,
using a hard gating mechanism with a different learning process which leads to
practical results on ImageNet.

###  2.4 Network Architecture Search

Network architecture search (NAS) first introduced in [36]. It focuses on
finding an optimal network architecture within a wide search space using
methods such as reinforcement learning. Many NAS works, such as [29],
typically regard the number of layer channels as part of the search space
while minimizing the FLOPs metric. The search space of [10] consists of layer
pruning ratios minimizing either FLOPs or memory metrics for a given target
accuracy. Other methods, such as [31], optimize the latency metric. Although
many NAS methods also share the idea of optimizing the number FLOPs, memory
footprint or latency, and even yield SOTA performance in many cases, they are
extremely computationally intensive compared to pruning methods, and thus
exhibit sub-optimal latency.

##  3 Approach

In this section we describe the components of Gator. We start with the gating
mechanism in Subsection 3.1, and then explain in Subsection 3.2 how modeling
the network architecture as a hypergraph encapsulates the channel dependency
between different layers in the network architecture, thereby allowing for
non-sequential pruning. Finally, in Subsection 3.3 we present the loss
function for computational cost and demonstrate how to customize it for FLOPs,
memory, and latency reduction.

###  3.1 Channel Gating

Figure 1: Logistic sigmoid gating model for an output of a layer.

Gator utilizes a gating mechanism which determines for each individual channel
whether it should be passed or gated for each individual sample within a
batch. The gating mechanism (see Eqn. 1) uses a binary version of the Gumbel
softmax distribution (introduced in [16]), referred to as the logistic sigmoid
distribution.

$\begin{split}&\sigma(\theta)=\frac{1}{1+e^{-\frac{\theta+x}{\tau}}}\quad{\rm
where}\quad x\sim Logistic(0,1)\\\
&g(\theta)=\begin{cases}1&\sigma(\frac{\theta+x}{\tau})\geq 0.5\\\
0&\sigma(\frac{\theta+x}{\tau})<0.5\\\ \end{cases}\simeq
sign(\theta+x)\;.\end{split}$ (1)

The gating binary value $g(\theta)$ is determined by adding the gating weight
$\theta$ to a random variable $x$ sampled from a logistic distribution. The
purpose of the sigmoid function is to enable a continuous gradient while the
threshold test provides a discrete output equivalent to a sign function in the
forward pass. The temperature parameter $\tau$, which controls the gradient
magnitude, is set to the default value of 1 and is, henceforth, ignored. The
addition of the random variable $x$ is according to the reparameterization
trick in [18], which enables a computational flow of the gradients for the
gating weights as part of the backpropagation process, while providing a
discrete value in the forward pass.

A gating module can be placed before or after a layer to pass or block
individual channels, thereby simulating the effect of pruning input or output
channels. Fig. 1 illustrates gating on an output of a model.

###  3.2 Hypergraph of Neural Network Channels

Early pruning methods on sequential network architectures, such as VGG [27],
focus naturally on individual convolutions, so they are applicable to all
layers. However, most recent architectures consist of many non-sequential
parts. Considering ResNet architectures, for example, note that inputs of some
convolutions depend on outputs of multiple convolutions residing in different
ResNet blocks, as a result of the “highway” path configuration. Thus, using a
channel pruning scheme, in an attempt to remove a channel from the output of a
single convolution may not remove the same channel from the input of a
dependent convolution or the output of another convolution, as implied by the
network’s non-sequential architecture. In other words, one needs to consider
the impact of a channel removal on the input and output channels of all
convolutions dependent on that channel. This aspect is illustrated for the
ResNet-18 architecture in Fig.2. We shall refer to a collection of channel-
dependent convolutions, as far as their input and output channels are
concerned, as a channel dependency group.

Figure 2: Layer 2 blocks 1,2 and layer 3 block 1 from ResNet-18 architecture,
demonstrating channel dependency between 6 different convolutions illustrated
by red lines; Thicker lines indicate the highway path which is the main cause
of this dependency.

We define an undirectional hypergraph $H=(V,E)$, where the set of vertices $V$
corresponds to the union of the input and output channel sets over all
convolutions. Let $C^{in}$ and $C^{out}$ denote these input and output sets,
respectively. Specifically, the vertices $c_{i}^{in}\in C^{in}$ and
$c_{i}^{out}\in C^{out}$ represent the (numbers of the) input and output
channels of convolution $i$, respectively. The set of edges $E$ is composed of
subsets of $V$, where an edge $e_{j}\in E$ corresponds to channel dependency
group $j$. (Each edge in $E$ is referred to as a dependency edge.) We
represent the dependency group as $e_{j}=C_{j}^{in}\cup C_{j}^{out}\subseteq
V$, where $C_{j}^{in}\subseteq C^{in}$ and $C_{j}^{out}\subseteq C^{out}$
denote, respectively, the input and output channel subsets associated with
dependency group $j$. We can now define pruning as a removal of channels from
each dependency edge.

For sequential networks, each dependency edge is composed of a pair of a
convolutional layer’s output and the following convolutional layer’s input.
However, we can now also map channel dependency for non-sequential networks,
such as ResNet-50, as detailed in Table 1. Specifically, 53 convolutions and a
fully-connected layer are mapped to 37 dependency edges, 32 of which are
trivial edges representing sequential connections inside bottleneck blocks,
and to 5 additional edges, which include more than two convolutions (some of
which are located along the “highway” path).

Table 1: Dependency mapping for ResNet-50; Each row represents one edge composed of output and input layers (vertices); encoding used, L1B1-(example of) layer and block, C0-first convolution, C1,C2,C3-bottleneck convolutions, D-downsampling convolution, and FC-fully connected; for brevity, trivial edges inside bottlenecks are listed once output layers | input layers | # edges
---|---|---
c0 | l1b1d, l1b1c1 | 1
c1 | c2 | 16
c2 | c3 | 16
l1b1d, l1b1c3, l1b2c3, l1b3c3 | l1b2c1, l1b3c1, l2b1d, l2b1c1 | 1
| l2b1d, l2b1c3, l2b2c3, l2b3c3,
---
l2b4c3
| l2b2c1, l2b3c1, l2b4c1, l3b1d,
---
l3b1c1
1
| l3b1d, l3b1c3, l3b2c3 l3b3c3,
---
l3b4c3, l3b5c3, l3b6c3
| l3b2c1, l3b3c1, l3b4c1, l3b5c1,
---
l3b6c1, l4b1d, l4b1c1
1
l4b1d, l4b1c3, l4b2c3, l4b3c3 | l4b2c1, l4b3c1, fc | 1

###  3.3 Loss Function for Computational Cost

For each edge $e_{j}$ we have $c_{j}(t)$ non-gated channels at iteration $t$
from a total of $c_{j}$ channels, computed from the binary gating values
defined in Eqn. (1):

$c_{j}(t)=\sum\limits_{i=1}^{c_{j}}g(\theta_{j}^{i})\;.$ (2)

Let $\lambda_{j}(t)$ denote the per-channel cost factor for iteration $t$. The
change in cost for different iterations stems from the fact that the number of
non-gated input and output channels of a convolution changes at each
iteration.

We now define the auxiliary computational loss

$L_{computation}(t)=\sum\limits_{j=1}^{m}c_{j}(t)\widehat{\lambda_{j}}(t)=\frac{\sum\limits_{j=1}^{m}c_{j}(t)\lambda_{j}(t)}{\sum\limits_{j=1}^{m}c_{j}\lambda_{j}}\;,$
(3)

where $m$ is the number of dependency edges, $\widehat{\lambda_{j}}(t)$ is the
normalized, per-channel cost factor, and $\lambda_{j}$ is the initial cost
factor. The normalization sets the initial loss to $L_{computation}(0)=1$, so
that it is easier to calibrate the total loss

$L_{total}(t)=L_{original}+\alpha L_{computation}(t)\;,$ (4)

where $L_{original}$ is the original loss and $\alpha$ is a global weight
being calibrated.

The memory footprint cost factor is the number of weights in all layers of
each individual edge, i.e.,

$\lambda_{j}(t)=\sum\limits_{q_{i}\in\bf{C}_{j}^{in}}k_{q_{i}}^{w}k_{q_{i}}^{h}c_{q_{i}}^{out}(t)+\sum\limits_{q_{i}\in\bf{C}_{j}^{out}}k_{q_{i}}^{w}k_{q_{i}}^{h}c_{q_{i}}^{in}(t)\;,$
(5)

where $\bf{C}_{j}^{out}$ and $\bf{C}_{j}^{in}$ are the convolution sets
corresponding to $C_{j}^{out}$ and $C_{j}^{in}$, respectively. For each
convolution $q_{i}$, $k_{q_{i}}^{w}$ and $k_{q_{i}}^{h}$ denote the kernel
width and kernel height, respectively. Also, $c_{q_{i}}^{in}(t)$ and
$c_{q_{i}}^{out}(t)$ denote, respectively, the number of non-gated input and
output channels at iteration $t$.

The theoretical speed cost factor is computed for a single input pixel, by
multiplying the number of weights of each layer by the downsample factor
$d_{i}$, i.e., the ratio between the number of pixels in the current layer’s
output and the original input image. Typically $d_{i}=4^{-s_{i}}$, where
$s_{i}$ is the number of stride-2 operations (e.g., convolutions, maxpooling,
etc.) encountered on a computational path from the input image to the current
layer $i$. Thus the FLOPs cost factor is given by

$\lambda_{j}(t)=\sum\limits_{q_{i}\in\bf{C}_{j}^{in}}d_{i}k_{q_{i}}^{w}k_{q_{i}}^{h}q_{i}^{out}(t)+\sum\limits_{q_{i}\in\bf{C}_{j}^{out}}d_{i}k_{q_{i}}^{w}k_{q_{i}}^{h}q_{i}^{in}(t)\;.$
(6)

For latency, the cost factor cannot be generally defined, as it is hardware-
specific. Instead, we approximate, in this case, the cost factor
$\lambda_{j}(t)$, by measuring the latency difference between the original
network and a network with $50\%$ of the channels in $e_{j}$ pruned.

###  3.4 Learning Rate Gradient Update Adjustment

During training, the loss function described in Eqn. 3, produces a gradient
update equivalent to the cost factor. The cost factor can be quite different
in scale for different layer. For FLOPs cost, for instance, the cost factor is
smaller for deeper layers (since they have more channels). As a result,
pruning of deeper layers will be slower (take more iterations) compared to
shallow layers. In turn, this could also lead to over pruning of shallow
layers, causing extensive adjustment of the rest of the network weights to
compensate, which eventually produces a sub optimal pruning result.

To combat both issues, we define a different time-dependent learning rate for
the weights of each individual dependency edge

$\eta_{j}(t)=\frac{\gamma\eta(t)}{\widehat{\lambda_{j}}(t)}\;,$ (7)

where $\eta_{j}(t)$ is the learning rate of dependency edge $j$, $\eta(t)$ is
the global learning rate, and $\gamma$ is a constant parameter.

Breaking down the gradient update for an individual gating weight $i$ on
dependency edge $j$, $\theta_{j}^{i}$, with a global pruning weight $\alpha$
we get

$\begin{split}\eta_{j}(t)\nabla_{\theta_{j}^{i}}&=\eta_{j}(t)\frac{\partial
L}{\partial\theta_{j}^{i}}=\eta_{j}(t)\frac{\partial L}{\partial
g(\theta_{j}^{i})}\frac{\partial g(\theta_{j}^{i})}{\partial\theta_{j}^{i}}\\\
&=\eta_{j}(t)\frac{\partial
g(\theta_{j}^{i})}{\partial\theta_{j}^{i}}\Big{(}\alpha\frac{\partial
L_{computational}}{\partial g(\theta_{j}^{i})}+\frac{\partial
L_{original}}{\partial g(\theta_{j}^{i})}\Big{)}\\\
&=\frac{\gamma\eta(t)}{\widehat{\lambda_{j}}(t)}\frac{\partial
g(\theta_{j}^{i})}{\partial\theta_{j}^{i}}\Big{(}\alpha\widehat{\lambda_{j}}(t)\frac{\partial
c_{j}(t)}{\partial g(\theta_{j}^{i})}+\frac{\partial L_{original}}{\partial
g(\theta_{j}^{i})}\Big{)}\\\ &=\gamma\eta(t)\frac{\partial
g(\theta_{j}^{i})}{\partial\theta_{j}^{i}}\big{(}\alpha+\frac{1}{\widehat{\lambda_{j}}(t)}\frac{\partial
L_{original}}{\partial g(\theta_{j}^{i})}\big{)}\;,\end{split}$ (8)

where we substituted, respectively, Eqn. 4 in the second line, Eqn. 3 and Eqn.
7 in the third line, and Eqn. 2 in the last line. The role of parameter
$\gamma$ is to balance out the large values of the term
$\frac{1}{\widehat{\lambda}(t)}$, since $\widehat{\lambda}(t)$ is generally
small.

In Eqn. 8 we observe the two components of the gradient produced by the
computational loss and the original loss. The first component attempts to gate
the channel, while the other retains as much relevant information as possible.
The gradient adjustment from Eqn. 7 is applied on both parts, which is crucial
for a balanced optimization. The last line shows how the gradient adjustment
cancels out $\widehat{\lambda_{j}}(t)$ for the pruning cost component, making
the impact of the pruning loss function on the gradient update similar for all
edges.

Another insight from Eqn. 8 is the self-pruning phenomenon we have
occasionally observed, where the original loss component contributes to the
removal of certain channels on its own (even with a zero pruning loss). This
could occur if a channel has a negative impact on the original loss, which in
turn, produces a larger update on the gating weight compared to the update on
the channel weight; thus, the channel is pruned faster then it can adjust. To
alleviate this issue, we calibrate $\gamma$ by performing training with no
pruning loss, and setting it to balance reasonable conversion time with
minimal self-pruning.

##  4 Experiments

###  4.1 Methodology and Setup

The lack of a comparison standard for pruning results is discussed in [1],
listing different types of pruning, different benchmark types, and a plethora
of different network architectures on which various pruning methods are
benchmarked. Hence, comparing between pruning methods is only relevant for
compatible criteria, i.e., same pruning type (e.g., structured pruning),
optimized parameter (e.g., FLOPs reduction), and architecture. (From a
practical consideration, we only compared with structured channel pruning
methods.) Similarly to most other works, we performed pruning on the ImageNet
ILSVRC2012 [25] and CIFAR-10 [19] classification challenges. However, we
believe that results for CIFAR-10 are not very indicative, due to the small
dataset, the large variation in accuracy [5], and the fact that the dataset
does not adequately represent “real-world” data. We chose the ResNet
architecture, which enabled demonstrating Gator’s pruning capability for a
non-sequential architecture and comparing it to other works. Note that newer
and smaller architectures would require a much greater effort in recreating
proper training and pruned versions, not to mention that there are very few
results to compare to.

In both experiments, we trained our model with the standard SGD optimizer,
setting the momentum to $0.9$ and the weight decay to $0.0001$. We initialized
the gating probability to $0.005$, which was found empirically to have a low
impact on the network’s initial accuracy, while still enabling the pruning in
a relatively low number of epochs. Note that a lower initial gating
probability would have required more steps to reach the pruning threshold
probability set to $0.5$. Crossing the threshold probability marks a channel
as permanently pruned. We performed iterative pruning with an increasing
global weight $\alpha$. Each pruning iteration consisted of training with
gating, and of additional fine-tuning without gating. The next iteration
loaded the weights produced by the gating phase.

###  4.2 Results for ResNet-50 on ImageNet

Figure 3: ImageNet (ILSVRC-12) top-5 accuracy (%) as a function of latency
speedup for ResNet-50 pruning by Gator compared with (a) other pruning
networks and (b) smaller networks for larger scale pruning. Table 2:
Comparison of pruned ResNet-50 on ImageNet for Gator optimized on FLOPs and
GPU latency; number of Gator variant indicates pruning weight $\alpha$;
results are ordered and divided to categories by FLOPs reduction ($\%$)

| Pruning
---
Algorithm
| Top-1
---
acc.%
| Top-5
---
acc.%
| FLOPs %
---
reduction
| Memory %
---
reduction
| Practical
---
Speed-up 1
channel pruning [11] | 72.30 | 90.80 | 24.14 | 15.52 | 1.069
Gator latency 0.25 | 76.00 | 92.94 | 26.75 | 5.22 | 1.310
Gator flops 0.25 | 75.89 | 92.93 | 31.48 | 8.88 | 1.200
DCP 30 [35] | 76.352 |  93.09 2 | 36.04 | 33.45 | 1.145
Geometric 30% [9] | 75.59 | 92.63 | 41.89 | 37.64 | 1.162
Gator latency 0.5 | 75.28 | 92.50 | 48.27 | 18.05 | 1.647
Gator flops 0.5 | 75.19 | 92.61 | 49.91 | 18.96 | 1.442
Geometric 40% | 74.83 | 92.32 | 53.16 | 47.72 | 1.251
DCP 50 | 74.88 | 92.19 | 55.77 | 51.63 | 1.477
PCAS [32] | 72.68 | 91.09 | 56.70 | 51.20 | –
Gator latency 1 | 74.24 | 91.95 | 61.17 | 31.18 | 1.857
Gator flops 1 | 74.14 | 91.99 | 62.64 | 31.40 | 1.612
DCP 70 | 72.65 | 91.03 | 71.10 | 66.00 | 1.569
Gator flops 2 | 72.36 | 90.97 | 76.57 | 55.15 | 2.106

  * 1

Speedup multiplier for batch 32 on single GTX 1080 Ti GPU

  * 2

Better accuracy then baseline ResNet-50.

Table 3: Comparison of small networks to high rate pruned ResNet-50 on ImageNet for Gator optimized on GPU latency; number of Gator variant indicates pruning weight $\alpha$; results are ordered and divided to categories by speedup | Pruning
---
Algorithm
| Top-1
---
acc.%
| Top-5
---
acc.%
FLOPs |  | Speed
---
up
Gator latency 2 | 72.69 | 90.98 | 973M | 2.273
Mobilenet V2 | 71.87 | 90.29 | 301M | 2.399
squeezenet 1.0 | 58.10 | 80.42 | 888M | 2.843
Gator latency 4 | 69.48 | 89.3 | 551M | 2.959
Gator latency 16 | 61.07 | 83.72 | 191M | 4.777
squeezenet 1.1 | 58.19 | 80.62 | 387M | 5.070
Gator latency 20 | 56.89 | 80.67 | 134M | 5.404

For pruning ResNet-50 trained on ImageNet, we used the pretrained weights. The
batch size was 256, which was split between two GPUs. However, due to the lack
of RAM, we used gradient accumulation and split each half batch into 4 parts,
training equivalenty on 8 GPUs with sub-batches of size 32 each. (This only
affects the batch normalization.) We trained with gates for $30$ epochs, using
learning rates of 0.01 and 0.001 for the first 20 epochs and the latter 10
epochs, respectively. Table 2 contains fine-tuning results for lower pruning
rates ($\alpha\leq 2$). We trained, in this case, for $40$ epochs with a
learning rate step policy of 0.001, 0.0001, and 0.00001 for 20,10,and 10
epochs, respectively. For high pruning rates, we trained for 20 epochs (given
the available resources under the time constraints), with a learning rate step
policy of 0.001, 0.0001, and 0.00001 for 10,5,and 5 epochs, respectively (see
Table 3). For image preprocessing, we applied mean subtraction and standard
deviation normalization for each color channel. Data augmentation consisted of
selecting a random $224\times 224$ crop from the original image and applying
to it a random horizontal flip. Practical speedup was measured by timing an
inference run in Pytorch on GTX 1080 Ti GPU with a batch of size 32 for
$224\times 224$ images.

Table 2 contains comparative results of optimized Gator pruning on FLOPs and
latency, with other methods (also illustrated in Fig. 3(a)). Note the very
similar accuracy obtained by the two optimization methods and how each type is
better suited for its designated reduction metric. (Although there is a single
instance that outperforms Gator (optimized on FLOPs reduction) for lower
pruning rates, note that the accuracy reported in this case [35] exceeds the
original ResNet-50 accuracy; this could be explained by the use of auxiliary
classification losses, known to enhance the performance [28]. Thus, we believe
this isolated instance is not a fair comparison.)

Regarding the latency speedup, Gator outperforms all other methods, achieving
a speedup of 1.4 for some pruning rates.

Results for higher pruning rates are available in Table 3. We compared Gator’s
latency optimization to that of much smaller architectures, such as SqueezeNet
[14] and MobileNet V2 [26], a newer architecture containing speed-efficient
elements, such as inverse bottlenecks and grouped convolutions. According to
Fig. 3(b), our method slightly outperforms both of these small architectures.

###  4.3 Results for ResNet-56 on CIFAR-10

We trained ResNet-56 on CIFAR-10 from scratch for $240$ epochs with a batch
size of $128$. The learning rate was $0.1$ for the first $120$ epochs, and was
changed to 0.01 and 0.001 at epochs 120 and 180, respectively. For pruning, we
trained with gating for $90$ epochs, starting with at a learning rate of
$0.01$ and decreasing it to $0.001$ at epoch $60$. This was followed by fine-
tuning for another $30$ epochs at the same learning rate. Image preprocessing
consisted of mean subtraction and standard deviation division for each
channel. Data augmentation consisted of random cropping and random horizontal
flipping. The results are available in Table 4, where we see a significant
reduction of FLOPs with little accuracy loss for Gator which also and best
performance on memory reduction compared to other models. In addition Gator
was able to prune entire ResNet blocks, thereby reducing the number of layers.

Table 4: Pruning results of ResNet-56 on CIFAR-10 using Gator optimization on FLOPs and memory vs. other methods; baseline network accuracy is listed for comparison as it differs across methods | Pruning
---
Algorithm
| Baseline
---
acc. %
| Acc.
---
%
| FLOPs %
---
reduction
| Memory %
---
reduction
Gator flops 1 | 93.63 | 93.58 | 43.39 | 33.50
Gator flops 2 | 93.63 | 93.14 | 54.62 | 44.39
Gator memory 0.25 | 93.63 | 93.66 | 33.21 | 42.70
Gator memory 0.5 | 93.63 | 92.59 | 57.98 | 75.08
DCP | 93.80 | 93.49 | 49.75 | 49.24
DCP Adapt | 93.80 | 93.81 | 47.09 | 70.33
Geometric 40 | 93.59 | 93.49 | 52.30 | –
Channel Pruning | 92.80 | 91.80 | 50.00 | –
PCAS | 93.04 | 93.58 | 54.80 | 53.70
Gator flops 16 | 93.63 | 91.50 | 79.33 | 73.40
Gator memory 2 | 93.63 | 90.99 | 76.05 | 90.87

##  5 Conclusions

In this paper we presented Gator, a novel structural channel-pruning method,
which utilizes dynamic channel gating during training to obtain smaller
effective networks. We showed how to prune non-sequential architectures by
using a hyper-graph formulation to encapsulate the mapping of dependencies
between input and output channels of different convolutions. We also
demonstrated how Gator can accommodate a variety of auxiliary cost functions
to optimize FLOPs, memory, and latency reduction due to pruning. Latency
reduction, in particular, is important for allowing to produce very light
networks designed for small devices. Finally, we showcased SOTA results of
pruning ResNet-50, trained on ImageNet, with respect to FLOPs reduction and
practical speedup on actual hardware.

## References

  * [1] Blalock, D., Ortiz, J.J.G., Frankle, J., Guttag, J.: What is the state of neural network pruning? arXiv preprint arXiv:2003.03033 (2020)
  * [2] Chen, W., Wilson, J., Tyree, S., Weinberger, K., Chen, Y.: Compressing neural networks with the hashing trick. In: International Conference on Machine Learning. pp. 2285–2294 (2015)
  * [3] Courbariaux, M., Hubara, I., Soudry, D., El-Yaniv, R., Bengio, Y.: Binarized neural networks: Training deep neural networks with weights and activations constrained to + 1 or -1. arXiv preprint arXiv:1602.02830 (2016)
  * [4] Ehteshami Bejnordi, B., Blankevoort, T., Welling, M.: Batch-shaping for learning conditional channel gated networks. arXiv pp. arXiv–1907 (2019)
  * [5] Frankle, J., Carbin, M.: The lottery ticket hypothesis: Finding sparse, trainable neural networks. arXiv preprint arXiv:1803.03635 (2018)
  * [6] Gain, A., Kaushik, P., Siegelmann, H.: Adaptive neural connections for sparsity learning. In: The IEEE Winter Conference on Applications of Computer Vision. pp. 3188–3193 (2020)
  * [7] Gao, X., Zhao, Y., Dudziak, Ł., Mullins, R., Xu, C.z.: Dynamic channel pruning: Feature boosting and suppression. arXiv preprint arXiv:1810.05331 (2018)
  * [8] Han, S., Pool, J., Tran, J., Dally, W.: Learning both weights and connections for efficient neural network. Advances in Neural Information Processing Systems 28, 1135–1143 (2015)
  * [9] He, Y., Liu, P., Wang, Z., Hu, Z., Yang, Y.: Filter pruning via geometric median for deep convolutional neural networks acceleration. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 4340–4349 (2019)
  * [10] He, Y., Lin, J., Liu, Z., Wang, H., Li, L.J., Han, S.: Amc: Automl for model compression and acceleration on mobile devices. In: Proceedings of the European Conference on Computer Vision. pp. 784–800 (2018)
  * [11] He, Y., Zhang, X., Sun, J.: Channel pruning for accelerating very deep neural networks. In: Proceedings of the IEEE International Conference on Computer Vision. pp. 1389–1397 (2017)
  * [12] Hu, J., Shen, L., Sun, G.: Squeeze-and-excitation networks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 7132–7141 (2018)
  * [13] Hua, W., Zhou, Y., De Sa, C.M., Zhang, Z., Suh, G.E.: Channel gating neural networks. In: Advances in Neural Information Processing Systems. pp. 1886–1896 (2019)
  * [14] Iandola, F.N., Han, S., Moskewicz, M.W., Ashraf, K., Dally, W.J., Keutzer, K.: Squeezenet: Alexnet-level accuracy with 50x fewer parameters and < 0.5 MB model size. arXiv preprint arXiv:1602.07360 (2016)
  * [15] Jaderberg, M., Vedaldi, A., Zisserman, A.: Speeding up convolutional neural networks with low rank expansions. arXiv preprint arXiv:1405.3866 (2014)
  * [16] Jang, E., Gu, S., Poole, B.: Categorical reparameterization with gumbel-softmax. arXiv preprint arXiv:1611.01144 (2016)
  * [17] Kim, Y.D., Park, E., Yoo, S., Choi, T., Yang, L., Shin, D.: Compression of deep convolutional neural networks for fast and low power mobile applications. arXiv preprint arXiv:1511.06530 (2015)
  * [18] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
  * [19] Krizhevsky, A., Hinton, G.: Learning multiple layers of features from tiny images (2009)
  * [20] Lebedev, V., Ganin, Y., Rakhuba, M., Oseledets, I., Lempitsky, V.: Speeding-up convolutional neural networks using fine-tuned CP-decomposition. arXiv preprint arXiv:1412.6553 (2014)
  * [21] LeCun, Y., Denker, J., Solla, S.: Optimal brain damage. Advances in Neural Information Processing Systems 2, 598–605 (1989)
  * [22] Li, H., Kadav, A., Durdanovic, I., Samet, H., Graf, H.P.: Pruning filters for efficient convnets. arXiv preprint arXiv:1608.08710 (2016)
  * [23] Molchanov, P., Tyree, S., Karras, T., Aila, T., Kautz, J.: Pruning convolutional neural networks for resource efficient inference. arXiv preprint arXiv:1611.06440 (2016)
  * [24] Rastegari, M., Ordonez, V., Redmon, J., Farhadi, A.: Xnor-net: Imagenet classification using binary convolutional neural networks. In: European Conference on Computer Vision. pp. 525–542. Springer (2016)
  * [25] Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z., Karpathy, A., Khosla, A., Bernstein, M., et al.: Imagenet large scale visual recognition challenge. International Journal of Computer Vision 115(3), 211–252 (2015)
  * [26] Sandler, M., Howard, A., Zhu, M., Zhmoginov, A., Chen, L.C.: Mobilenetv2: Inverted residuals and linear bottlenecks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 4510–4520 (2018)
  * [27] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556 (2014)
  * [28] Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., Wojna, Z.: Rethinking the inception architecture for computer vision. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 2818–2826 (2016)
  * [29] Tan, M., Le, Q.V.: Efficientnet: Rethinking model scaling for convolutional neural networks. arXiv preprint arXiv:1905.11946 (2019)
  * [30] Veit, A., Belongie, S.: Convolutional networks with adaptive inference graphs. In: Proceedings of the European Conference on Computer Vision. pp. 3–18 (2018)
  * [31] Wu, B., Dai, X., Zhang, P., Wang, Y., Sun, F., Wu, Y., Tian, Y., Vajda, P., Jia, Y., Keutzer, K.: Fbnet: Hardware-aware efficient convnet design via differentiable neural architecture search. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 10734–10742 (2019)
  * [32] Yamamoto, K., Maeno, K.: Pcas: Pruning channels with attention statistics for deep network compression. arXiv preprint arXiv:1806.05382 (2018)
  * [33] Zhou, A., Yao, A., Guo, Y., Xu, L., Chen, Y.: Incremental network quantization: Towards lossless cnns with low-precision weights. arXiv preprint arXiv:1702.03044 (2017)
  * [34] Zhu, C., Han, S., Mao, H., Dally, W.J.: Trained ternary quantization. arXiv preprint arXiv:1612.01064 (2016)
  * [35] Zhuang, Z., Tan, M., Zhuang, B., Liu, J., Guo, Y., Wu, Q., Huang, J., Zhu, J.: Discrimination-aware channel pruning for deep neural networks. In: Advances in Neural Information Processing Systems. pp. 875–886 (2018)
  * [36] Zoph, B., Le, Q.V.: Neural architecture search with reinforcement learning. arXiv preprint arXiv:1611.01578 (2016)


# A Random Point Initialization Approach to Image Segmentation with
Variational Level-sets

J.N. Mueller, J.N. Corcoran
University of Colorado Boulder Postal Address: Department of Applied
Mathematics, University of Colorado, Box 526 Boulder CO 80309-0526, USA;
email<EMAIL_ADDRESS>phone: 303-492-0685

###### Abstract

Image segmentation is an essential component in many image processing and
computer vision tasks. The primary goal of image segmentation is to simplify
an image for easier analysis, and there are two broad approaches for achieving
this: edge based methods, which extract the boundaries of specific known
objects, and region based methods, which partition the image into regions that
are statistically homogeneous. One of the more prominent edge finding methods,
known as the _level set method_ , evolves a zero-level contour in the image
plane with gradient descent until the contour has converged to the object
boundaries. While the classical level set method and its variants have proved
successful in segmenting real images, they are susceptible to becoming stuck
in noisy regions of the image plane without _a priori_ knowledge of the image
and they are unable to provide details beyond object outer boundary locations.
We propose a modification to the variational level set image segmentation
method that can quickly detect object boundaries by making use of _random
point initialization_. We demonstrate the efficacy of our approach by
comparing the performance of our method on real images to that of the
prominent Canny Method.

00footnotetext: Keywords: image segmentation, level set methods, computer
vision, Canny method
AMS Subject classification: 68U10, 62M40

##  1 Introduction

Image segmentation is one of the most challenging tasks in image processing,
which is essential for myriad computer vision tasks including facial
recognition, traffic control systems management, autonomous vehicle
development, robotics, and medical imaging and tomography. The primary goal of
image segmentation is to simplify the image for easier analysis. This can be
done in several ways, for example, by breaking up an image into homogeneous
regions made up of pixels with some measure of similarity or detecting edges
based on continuity changes in pixel intensity.

The development of the first techniques for image segmentation can be traced
back more than 50 years. In 1965, an operator for detecting edges between
different parts of an image, known as the the Roberts edge detector (also
referred to as the Roberts operator), was introduced [1] as a first step
toward decomposing an image into its constituent components. Since the
introduction of this operator, many approaches to image segmentation that
range from clustering algorithms to curve evolution have been discovered.
Broadly speaking, they can be split into two categories. The first is to find
the boundaries of specific known objects, for example, in MRI scans.
Approaches on this front have been very successful as in [2], [3], and [4].
The second problem is to partition an image into disjoint regions that are
statistically homogeneous. Based on an image alone, this problem is generally
ill-posed with potentially multiple reasonable partitions. This is especially
true for low quality images where a lack of assumptions about smoothness will
result in a noisy segmentation.

##  2 Active Contours

The active contour approach to edge detection is a partial differential
equation (PDE)-based segmentation method that treats image segmentation as an
energy optimization task. Broadly speaking, an initial curve (or set of
curves) is defined and driven into the interior of the image plane along the
negative gradient of an energy functional using gradient descent methods [5,
6, 7, 8, 9]. In particular, the energy functional is composed of terms
representing the _internal energy_ of the evolving contour and _external
energy_ of the image [5]. The external energy is computed from underlying
image properties such as image gradients and acts as a “shrinkage” term,
drawing the contour into the image plane toward object boundaries, while the
internal energy acts as a regularization term which controls the smoothness
and continuity of the deforming curve [5, 9, 10]. The energy functional must
be carefully chosen so that the propagating curve stops once it encounters
object boundaries within the image. Often the chosen active contour models
include a gradient dependent edge-function, $g$, since object boundaries in
the image occur at the minimum (or critical) value(s) of the energy functional
[5].

Commonly, active contour methods either construct the initial curve explicitly
as a parameterized set of Lagrangian curves or implicitly represent the
initial curve as a certain contour of a higher dimensional function in an
Eulerian framework [6]. Explicit active contour methods were first introduced
in [11] and are robust with respect to image noise and boundary gaps [8].
However, explicit parameterization schemes are often sensitive to the initial
conditions used and struggle to handle significant changes in curve geometry
such as merging or splitting that may naturally arise in an image segmentation
task. Thus, Lagrangian methods rely on a continual reparameterization, which
usually becomes more complex as the segmentation task progresses, and they can
break down as the evolving curve’s shape changes. Methods which handle such
topological changes naturally and avoid the complexities of reparameterization
alternately define the curve implicitly as a level set.

###  2.1 Level Set Methods

The _level set method_ was introduced in [12] as a means of following
complicated dynamics in front propagation and was separately extended to the
task of edge detection in [13] and [14]. The basic idea is to represent the
initial curve implicitly by embedding it as the _zero-level contour_ (or zero-
level set) of a higher-dimensional hypersurface known as the _level set
function_ and denoted $\varphi$. That is, we define the curve of interest to
be the contour

$\mathcal{C}:=\big{\\{}\vec{x}\ |\ \varphi(\vec{x},t)=0\big{\\}}.$

The level set function (and, thus, the zero-level set curve) is optimized
according to a level set evolution equation such as

$\displaystyle\frac{\partial\varphi}{\partial t}$
$\displaystyle=-\frac{\partial E}{\partial t}=-F|\nabla\varphi|\,.$ (1)

A detraction of traditional level set methods is the shape of the curve can
become distorted as it evolves in the image plane so that portions of the
curve form “shocks,” becoming excessively steep or flat over time [6, 15]. The
formation of shocks can result in instability in the algorithm if unchecked
and the computed solution where shocks form is inaccurate. It follows the
evolving curve must be periodically reshaped (known as “re-initialization”)
according to some metric describing how deformed the curve has become [7].
However, there are no systematic criteria for determining when or how often a
curve should be reshaped [16] and in some cases reshaping can introduce
additional error by spatially shifting the reformed curve away from the true
hypersurface [6].

While early level set methods are “pure PDE methods” which construct the
evolution equation (1) by first defining an evolution PDE in the Lagrangian
framework and then converting it to an evolution PDE for the level set
function [6], a method known as _variational level sets_ obtains the desired
evolution equation by minimizing an energy functional defined on the level set
function itself [6, 9, 10]. First proposed in [17], this approach was
popularized by [5]. Not only does the variational level set approach allow the
evolution equation to be obtained in a more straightforward manner, but it
also allows additional details such as region-based information [5, 18] or
shape information [19] to be included in the energy functional.

More recently, variational level set methods such as [6], [7], [16], and [20]
have been used to address the problem of shock formation. In particular, these
methods which have emerged avoid the creation of shocks by including _distance
regularization terms_ in the energy functional. The purpose of the
regularization terms is to locally monitor the evolving curve and to penalize
deviations away from a _signed distance function_ which may eventually lead to
shock formation.111The choice of penalty term reflects the usual method of
initializing the level set curve $\varphi_{0}$ as a signed distance function.

For example, the variational form introduced in [6] and expanded upon in [7]
exploits the property that any signed distance function satisfies
$|\nabla\varphi|=1$ to define a metric

$\displaystyle\mathcal{P}$
$\displaystyle=\int_{\Omega}\frac{1}{2}\big{(}|\varphi|-1)^{2}d\vec{x}$ (2)

which measures deviations of the level set curve $\varphi$ from a signed
distance function in a neighborhood $\Omega\subset\mathbb{R}^{2}$ around the
desired zero-level set. This term depends only on $\varphi$ and represents the
internal energy of the curve.

The energy functional which must be minimized is then given by a _variational
formula_

$\displaystyle E$ $\displaystyle=\mu\mathcal{P}+\mathcal{E}_{\text{\tiny
ext}}$ (3)

where $\mu>0$ controls the effect of the penalty on the evolving curve and
$\mathcal{E}_{\text{\tiny ext}}$ is comprised of external and boundary
energies that depend only on the image data.

The gradient flow which propagates the level set curve $\mathcal{C}$ into the
image plane is computed using a first G$\hat{\text{a}}$teaux derivative [6, 7]
of the energy functional (3)

$\displaystyle\frac{\partial\varphi}{\partial t}$
$\displaystyle=-\frac{\partial
E}{\partial\varphi}=-\mu\frac{\partial\mathcal{P}}{\partial\varphi}-\frac{\partial\mathcal{E}_{\text{\tiny
ext}}}{\partial\varphi}\,.$ (4)

Specific expressions for $\mathcal{E}_{\text{\tiny ext}}$ and the completely
differentiated right-hand size of (4) are given in [6] while finite difference
schemes for implementing this method can be found in [21, 22].

##  3 The Random Point Initialization Method

###  3.1 Motivation

We motivate our modification to the variational level set method described in
[6] with the following example. Consider the image of a spider [23] given in
the far-left frame of Figure 1. To segment the image and obtain the boundary
of the spider, define an initial level set curve $\varphi_{0}$ as a step
function

$\displaystyle\varphi_{0}$ $\displaystyle=\begin{cases}\ \
c_{0},\quad\text{outside contour}\\\ -c_{0},\quad\text{inside
contour}\end{cases}$ (5)

where $c_{0}>0$ is constant and drive the curve into the image plane along the
energy gradient given by (4) as in [21, 22].

Although the level set appears to converge quickly to the object boundary
where the foreground (spider) and background are reasonably distinct and the
image gradient is smooth, the evolving curve stagnates around the spider’s
feet; even after tens of thousands of iterations the position of the level set
is relatively unchanged.

Figure 1: _Left_ Gray-scale image of spider with initial level set curve
(pink). _Center-left_ Level set curve (pink) after 500 steps. _Center-right_
Level set curve (pink) after 1,000 steps. _Right_ Level set curve (pink)
function after 50,000 steps.

Recognizing that the curve becomes stuck in regions where the image gradient
is not smooth, we attempted to prompt convergence of the curve to the object
boundary in these regions by periodically “kicking” the points of the level
set by a small, random amount
$\bm{\varphi_{k}}\mapsto\bm{\varphi}_{k}+\varepsilon,\ \ \varepsilon\sim
N(0,\sigma^{2})$. However, perturbing the level set in this way did not result
in a significant improvement to the convergence speed. Perturbations large
enough to noticeably displace the level set destroyed the curve and, while
future iterations collected the scattered points into a salient boundary once
again, experimental evidence suggests that the reformed boundary remained
stuck in troublesome regions of the image. On the other hand, perturbations
which locally displaced points were usually too small to nudge the evolving
curve to regions in the image where the gradient were smooth.

Figure 2: _Left_ Gray-scale image of spider with initial level set curves
(pink solid: iterated with variational level set method; green dashed:
iterated with perturbed level set method). _Center-left_ Level set curves
after 500 steps. _Center-right_ Level set curves after 1,000 steps. _Right_
Level set curves after 50,000 steps.

Additionally, the segmentation obtained by periodically perturbing the level
set often suffered from boundary leakage since it is possible for the
perturbation to push points on the evolving curve across the object boundary.
This is illustrated in the far-right frame of Figure 2, wherein the white tops
of the spider’s feet are excluded from the segmentation obtained with a
perturbed level set (green dashed) curve.

###  3.2 The Random Point Initialization

Our efforts to perturb the level set in this way provided the following key
insight: _points scattered by random perturbations are collected together in
future iterations of the algorithm._ Edge detection algorithms depend on the
image gradients to capture features of the underlying image topography and use
this information to drive the defined hypersurface and the embedded level set
curve to object boundaries (i.e. global minima) within the image plane.
Importantly, the algorithm uses this gradient information to propagate any
collection of points, regardless of whether these points define a coherent
curve or are randomly scattered. If the algorithm iterates a curve, then the
net affect is to contract the curve around the object(s) of interest (i.e. a
lassoing effect). If instead the algorithm iterates a set of free points, the
net affect is to push the points “downhill” along the image gradients until
they settle into minima of the image plane.

Figure 3: A nano-system of chemical reactions.

We take advantage of this and make the following modification to the
algorithm: instead of embedding a level set in a hypersurface as in (4), embed
an initial set of randomly scattered points $\bm{\phi}\in\mathbb{R}^{I\times
J}$ where $\phi_{ij}\sim N(0,\sigma^{2})$ and $\sigma$ is the size of the
perturbation. The desired object boundary (edge-set) resulting from iterating
$\bm{\phi}$ with gradient descent is still given by the zero-level set contour
after $k$ steps. We call this the _random point initialization_ (RPI) method.

To illustrate the RPI approach clearly, consider the image of a system of
chemical reactions depicted in Figure 3 [24]. Scattering points randomly in
the image plane (viewed on top of the grayscale image for reference) and
completing several gradient descent steps using [21, 22] results in the noisy
segmentation shown in the far-right frame of Figure 4.

Figure 4: _Left:_ Random points (red) scattered on top of the grayscale image.
_Center:_ The result of $k=4$ steps of the variational level set method
starting from randomly placed points, plotted with the original grayscale
image for reference. _Right:_ The result of $k=4$ steps of the variational
level set method starting from randomly placed points. The grayscale image has
been removed from the background to make the coalescence of boundary points
more apparent.

Clearly, the points coalesce about the desired object boundaries. However, the
segmented image is excessively noisy and the obtained boundary is diffuse. The
problem of filling boundary gaps is easily solved by constructing the desired
edge set from multiple runs of the RPI method. We demonstrate this idea first
by plotting multiple edge sets obtained with the RPI method on top of one
another. (The background image is omitted for clarity.)

Figure 5: Plotting the edge sets resulting from multiple RPI runs on top of
one another. _Left:_ A single RPI zero-level set. _Center:_ Four RPI zero-
level sets plotted on top of one another. _Right:_ Eight RPI zero-level sets
plotted on top of one another.

Notice that as the zero-level sets are layered together the boundaries around
objects within the image become more distinct.

To replicate this outcome in a systematic way, we periodically (after $k$
iterations) re-initialize the evolving set as a random matrix $\bm{\phi}_{0}$.
Each time we re-initialized the set to $\bm{\phi}_{0}$ and perform $k$
iterations, we complete a single run of the RPI method. The desired edge set
from each run is the zero-level set, as before, which can be flattened into an
associated vector of values.

$\displaystyle\bm{\phi}_{k}=\begin{bmatrix}\phi_{11}&\phi_{12}&\ldots&\phi_{1J}\\\
\phi_{21}&\phi_{22}&\ldots&\phi_{2J}\\\ \vdots&\vdots&\ddots&\vdots\\\
\phi_{I1}&\phi_{I2}&\ldots&\phi_{IJ}\end{bmatrix}\longrightarrow\begin{bmatrix}\phi_{11}\ldots\phi_{1J}&\phi_{21}\ldots\phi_{2J}&\ldots&\phi_{I1}\ldots\phi_{IJ}\end{bmatrix}$

We collect the flatten the edge sets from each RPI run by stacking the vectors
vertically in a segmentation matrix. Here the super-script denotes the RPI
run.

$\displaystyle\begin{bmatrix}\bm{\phi}_{k}^{(1)}\\\ \bm{\phi}_{k}^{(2)}\\\
\vdots\\\ \bm{\phi}_{k}^{(m)}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\phi_{11}^{(1)}\ldots\phi_{1J}^{(1)}&\ldots&\phi_{I1}^{(1)}\ldots\phi_{IJ}^{(1)}\\\
\phi_{11}^{(2)}\ldots\phi_{1J}^{(2)}&\ldots&\phi_{I1}^{(2)}\ldots\phi_{IJ}^{(2)}\\\
\vdots&\ddots&\vdots\\\
\phi_{11}^{(m)}\ldots\phi_{1J}^{(m)}&\ldots&\phi_{I1}^{(m)}\ldots\phi_{IJ}^{(m)}\\\
\end{bmatrix}$ $\displaystyle\hskip
25.29494pt\begin{matrix}\downarrow&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\ldots}&\downarrow&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\ldots}&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\ldots}&\downarrow&{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\ldots}&\downarrow\end{matrix}$
$\displaystyle\overline{\bm{\phi}}$ $\displaystyle=\hskip
7.22743pt\begin{bmatrix}\overline{\phi}_{11}\hskip 1.8063pt\ldots\hskip
1.8063pt\overline{\phi}_{1J}\hskip 1.8063pt&\ldots&\hskip
7.22743pt\overline{\phi}_{I1}\hskip 3.61371pt\ldots\hskip
3.61371pt\overline{\phi}_{IJ}\end{bmatrix}$

Once the iterations have been collected we perform element-wise averaging of
the sets (that is, the average of each column) to find the mean “image pixel”
value in the segmentation matrix. The vector that results from this averaging
process can then be reshaped into an $I\times J$ segmentation matrix that is
the zero-level set of interest.

$\displaystyle\overline{\bm{\phi}}=\begin{bmatrix}\overline{\phi}_{11}\hskip
1.8063pt\ldots\hskip 1.8063pt\overline{\phi}_{1J}\hskip 1.8063pt&\ldots&\hskip
7.22743pt\overline{\phi}_{I1}\hskip 3.61371pt\ldots\hskip
3.61371pt\overline{\phi}_{IJ}\end{bmatrix}\longrightarrow\begin{bmatrix}\overline{\phi}_{11}&\overline{\phi}_{12}&\ldots&\overline{\phi}_{1J}\\\
\overline{\phi}_{21}&\overline{\phi}_{22}&\ldots&\overline{\phi}_{2J}\\\
\vdots&\vdots&\ddots&\vdots\\\
\overline{\phi}_{I1}&\overline{\phi}_{m2}&\ldots&\overline{\phi}_{IJ}\end{bmatrix}$

The segmentation of the chemical system resulting from combining multiple runs
of the RPI method (including several post-processing steps discussed in Sec.
3.3) is given in Figure 6.

Figure 6: Segmentation of the chemical system with $m$-random point
initialization steps.

We emphasize that the random re-initialization is _not_ the same as the re-
initialization methods found in the literature and discusses in Section 2, but
is instead a numerically expedient means of repeatedly segmenting the same
image and combining the resulting edge sets.

###  3.3 Edge set Post Processing

Once the pixel-averaged segmentation matrix $\bm{\overline{\phi}}$ has been
constructed from multiple RPI runs, we found that the edge set usually
contains some amount of noise from random points settling in local minima of
the image and becoming stuck. Furthermore, we noticed that some points
generally settled very close to the true boundary (true zero-level set) but
were not pushed all the way into the minimum by the $k$th iteration. Further
de-noising and smoothing of the segmented image is done with a series of post-
processing steps.

####  3.3.1 Collecting Near-Boundary Points

We collect near-boundary points and include them in the edge set by
thresholding, i.e. defining a small range of values around the zero-level
contour that we wish to include in the edge set. The set of points found in
this range, which we denote $\bm{\overline{\phi}_{\text{\tiny Th}}}$, are then
used as an approximation to the true zero-level set.

As an aside, we found it was most effective to first normalize the values in
$\bm{\overline{\phi}}$ before extracting the (approximate) zero-level set
$\bm{\overline{\phi}_{\text{\tiny Th}}}$. We used the normalization

$\displaystyle\dfrac{\bm{\overline{\phi}}}{||\bm{\overline{\phi}}||}$
$\displaystyle=\dfrac{\bm{\overline{\phi}}}{\max_{i,j}|\bm{\overline{\phi}}|+\min_{i,j}|\bm{\overline{\phi}}|}$
(6)

but we note that other choices for normalization, such as dividing by the
absolute maximum pixel value, appeared to yield similarly good results.

Note that choosing

$\displaystyle\bm{\overline{\phi}_{\text{\tiny
Th}}}=\bigg{\\{}\frac{\bm{\overline{\phi}}}{||\bm{\overline{\phi}}||}\
\bigg{|}\ p_{\text{\tiny
low}}\leq\frac{\bm{\overline{\phi}}}{||\bm{\overline{\phi}}||}\leq
p_{\text{\tiny up}}\bigg{\\}}$ (7)

eliminates much of the noise in the edge set, as points which have settled in
local minimum values will be automatically excluded from this set.

Obtaining the edge set in (7) using conditional relationships results in an
edge set of binary entries of the form

$(\overline{\phi}_{i,j})_{\text{\tiny
Th}}=\begin{cases}1,\qquad\overline{\phi}_{i,j}\in[p_{\text{\tiny
low}},p_{\text{\tiny up}}]\\\ 0,\qquad\text{otherwise}\end{cases}$

where pixels of value $1$ are black and pixels of value $0$ are white. The
final post-processing steps used to thin and smooth the edge set require
binary entries such as these.

####  3.3.2 Further Curve Smoothing and Image De-noising

Segmentation algorithms are evaluated by both the accuracy and quality of the
obtained boundary, usually defined by the thickness and continuity of the
curve in the edge set. Since the multiple-RPI method is based on collecting
free points along edges within the image plane, the obtained segmented image
often has an undesirably thick boundary. Once the zero-level edge set
$\displaystyle\bm{\overline{\phi}_{\text{\tiny Th}}}$ is constructed we
further refine the segmented image by smoothing and thinning the edges.

The edges are smoothed using morphological operations which examine small non-
overlapping neighborhoods of pixels. In particular, the operator compares the
central pixel $p$ in a $3\times 3$ neighborhood to its eight-neighbor pixels
and assigns to $p$ the value ($0$ or $1$) shared by a majority of its
neighbors. In the event of a tie, $p$ is treated as a pixel of noise and is
assigned $p=0$ [25]. This removes extraneous pixels from the segmented image
while also filling in small spurious gaps in the boundary.

The edges can be thinned with a “skeletonizing” algorithm such as [26] (first
introduced in [27] but later modified in [26] to reduce the likelihood of
introducing boundary gaps). The thinning algorithm uses contour-following
methods to count the number of times a pixel is traced. Pixels only traced
once are deleted (set to $0$), while pixels traced multiple times are
considered essential to maintaining connectivity in the boundary (keep the
value of $1$). In a basic sense, the algorithm thins the boundary while
retaining two fundamental types of geometry: diagonal lines and $2\times 2$
squares [26, 27]. While the algorithm may be run until the image stops
changing, in practice it is often only necessary to run two or three
iterations to sufficiently reduce the thickness of the boundary found in these
post-processing steps.

###  3.4 Other modifications

At times scatting $IJ$ random points in the image plane can lead to points
piling up in certain locations of the image plane, resulting in object
boundaries that are excessively thick or that merge with the boundaries of
nearby objects. When this happens, it may obscure subtle details in the image
that we may wish to capture.

One solution which often improves the result of segmenting with the $m$-RPI
approach is to re-size the image using bicubic interpolation to spatially
separate objects within the image. By pushing the boundaries of objects apart
from one another, points scattered in the image plane collect around
individual objects since the objects are distinct and reduces the likelihood
of merged boundaries.

Alternately, the number of random points scattered in the image plane may be
reduced by some fraction $\alpha$ to $\alpha IJ,\ \ 0<\alpha<1$, of the total
number of pixels in the image.

$\longrightarrow$ $\longrightarrow$

Figure 7: Left A dense $I\times J$ matrix $\bm{\phi}$ of random values. Center
A block-matrix sparse $\bm{\phi}$ with $\frac{1}{4}IJ$ random entries. Right
The sparse version of $\bm{\phi}$ with the rows and columns randomly permuted
to distribute the non-zero entries throughout the full image matrix.

However, using a sparse-$\bm{\phi}_{0}$ set for noisy or complicated images
can result in unwanted boundary gaps in the final edge set. For many images, a
useful approach is to scatter $IJ$ random points for the first RPI run, which
in theory should find the general locations of object boundaries. The
subsequent RPI runs can then be re-initialized with sparse random sets, which
can be though of as filling in gaps along the boundary without adding too many
extra points to regions where objects are close together.

###  3.5 Post-processing example

We demonstrate the efficacy of post-processing by segmenting the synthetic
image shown in Figure 8 using the $m$-RPI method.

Figure 8: The original synthetic image

We initialize the first run of the RPI method using a dense random point set
and re-initialize subsequent RPI runs with sparse random sets where
$\alpha=0.25$ is the proportion of non-zero elements. Both the dense and
sparse random sets have entries distributed normally as $\phi_{ij}\sim
N(0,\sigma^{2})$ where $\sigma=0.01$. We restart the method every $k=8$
iterations and iterate until a total of 15 RPI runs are completed. Note that
each post-processing step is applied sequentially to the segmentation matrix
in the order presented.

We start by averaging $15$ $m$-RPI runs to find the average segmentation
matrix $\overline{\bm{\phi}}$ as described in Section 3. The resulting
segmentation in Figure 9 clearly outlines the boundary of object in the image,
but this segmentation matrix also contains excessive noise.

Figure 9: Segmentation for Figure 8 after averaging $15$ runs of the $m$-RPI
method.

Normalizing the entries of $\overline{\bm{\phi}}$ as in (6) pushes all of the
values to the range
$0\leq\big{|}\frac{\bm{\overline{\phi}}}{||\bm{\overline{\phi}}||}\big{|}\leq
1$. This does not de-noise the image but it accentuates the boundary of the
object we wish to find.

Figure 10: Normalization of $\overline{\bm{\phi}}$.

We then threshold to approximate the zero-level set
$\overline{\bm{\phi}}_{\text{\tiny log}}$ as in (7) by taking a small range of
values around zero. In this example we use $[p_{\ell},p_{u}]=[-0.175,0.075]$.
Experimental evidence suggests that using a larger proportion of negative
entries (values within the object in question) gives a cleaner segmentation.

Figure 11: Narrow band approximation to the zero-level set.

The thresholding step effectively de-noises the edge set but also results in a
_doubled_ boundary curve around the object. We apply both curve smoothing and
thinning to fill and thin the boundary, producing the final image segmentation
for the sample image. (Note that the built-in Matlab® function bwmorph was
used to implement the curve smoothing and thinning routines. The morphological
operation majority with the additional option Inf was used to smooth the curve
while thin, which implements the thinning algorithm described in Section
3.3.2, was used to skeletonize the boundary.)

Figure 12: _Left:_ The ground truth segmentation of the image in Figure 8.
_Right:_ The final segmentation obtained with the $m$-RPI method after the
boundary curve has been smoothed and thinned.

The expected segmentation of the image (the “ground truth”) is included in
Figure 12 for comparison.

##  4 Experimental Results: Comparison to the Canny Method

Returning to the image of the spider from Section 3 and now using the $m$-RPI
method, we successfully segmented the image using $15$-RPI runs. In addition
to finding the full boundary of the spider, including the troublesome region
around his legs, the $m$-RPI method located additional details within the
image – such as the pattern of spots on the spider’s back and texture of the
leaf he is standing on – that were previously missed by the classical level
set method.

Figure 13: Comparison of the $m$-RPI method (center-left) to the Canny method
when the Canny parameters are both default (center-right) and optimized (far-
right).

As image segmentation is an ill posed problem, each segmentation task depends
heavily on the image that will be segmented (“synthetic” or “real”) and the
type of information sought from the image. In general, it is difficult to
establish any consistent criteria for evaluating segmentation algorithms.
While synthetic images often allow for pixel-to-pixel comparisons of a ground
truth to the segmented image, these are not representative of most real images
of interest and are not useful for evaluating an algorithm on real datasets.

To evaluate the performance of the $m$-RPI method on the spider image, we make
a visual comparison to the segmentation results found using the Canny Edge
Detector, which is considered state-of-the-art. Notice that two Canny results
are given. In an effort to compare the methods as honestly as possible,
significant time and effort was dedicated to optimizing the parameters used in
the Canny Edge Detector provided by Matlab.222We chose an upper thresholding
value of $0.2$ and a lower thresholding value of $0.1$ as optimal values. It
should be noted that, prior to optimizing the Canny thresholding values in
Matlab, the $m$-RPI method equalled or surpassed the Canny segmentation of
most images that we examined. When the $m$-RPI results for the spider are
compared to the results of the optimized Canny method, we see similar levels
of detail in the segmentation. However, the $m$-RPI method clearly outperforms
the Canny segmentation obtained with the default parameter settings.

We conclude by demonstrating the results of segmenting real images with the
$m$-RPI method and we visually compare our results to the Canny Edge
Detector.333Test images belong to the authors and may not be re-used without
permission.

##  5 Conclusion

The Canny Method is considered the state-of-the-art image segmentation
algorithm, but the performance of this method depends heavily on the chosen
upper and lower thresholding values that are used. Careful selection of the
parameters used in the Canny algorithm is required for the method to
consistently perform well, and therefore significant experimentation is
required to optimize the thresholding values. Because of this, our algorithm
has the advantage of successfully segmenting many real image without extensive
parameter tuning and may be easily implemented for many image segmentation
tasks.

While the Canny Method in Matlab and Python segments images more quickly
(clock time) than our current $m$-RPI code does – for example, the spider
image was segmented with Matlab’s built-in Canny in approximately 6 seconds
while the $m$-RPI method with post-processing took approximately 29 second –
the $m$-RPI method promises substantial speed-up under a parallelization of
the code as no RPI run depends on any of the others.

###  4.1 Comparison to Canny Segmentation

## References

  * [1] Roberts LG.  Machine perception of three-dimensional solids.  In: et al JT, editor. Optical and electro-optical information processing; 2018. p. 159-97.
  * [2] Bouman CA, Shapiro M.  A multiscale random field model for Bayesian image segmentation.  IEEE Transactions on Image Processing. 1994;3(2):162-77.
  * [3] Cheng W, Ma L, Yang T, Liang J, Zhang Y.  Joint Lung CT Image Segmentation: A Hierarchical Bayesian Approach.  PLoS ONE. 2016;11(9).
  * [4] Puonti O, Iglesias JE, , Leemput KV.  Fast and Sequence-Adaptive Whole-Brain Segmentation Using Parametric Bayesian Modeling.  NeuroImage. 2016;143:235-49.
  * [5] Chan T, Vese LA.  Active Contours without Edges.  IEEE Transactions on Image Processing. 2001;10.
  * [6] Li C, Xu C, Gui C, Fox MD.  Level set Evolution without Re-initialization: A New Variational Formulation.  2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05). 2005;1:430-6.
  * [7] Li C, Xu C, Gui C, Fox MD.  Distance Regularized Level Set Evolution and Its Application to Image Segmentation.  IEEE Transactions on Image Processing. 2010:3243-54.
  * [8] Weeratunga SK, Kamath C.  An Investigation of Implicit Active Contours for Scientific Image Segmentation.  Proc SPIE 5308, Visual Communications and Image Processing. 2004.
  * [9] Yuan Y, He C.  Variational level set methods for image segmentation based on both L2 and Sobolev gradients.  Nonlinear Analysis: Real World Applications. 2012;13(2):959-66.
  * [10] Cremers D, Rousson M, Deriche R.  Review of Statistical Approaches to Level Set Segmentation: Integrating Color, Texture, Motion, and Shape.  International Journal of Computer Vision. 2007;72(2):195-215.
  * [11] Kass M, Witkin A, Terzopoulos D.  Snakes: Active Contour Models.  International Journal of Computer Vision. 1987;1:321-31.
  * [12] Osher S, Sethian J.  Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations.  Journal of Computational Physics. 1977;79:12-49.
  * [13] Caselles V, Cattè F, Coll T, Dibos F.  A geometric model for active contours in image processing.  Numer Math. 1993;66:1-31.
  * [14] Malladi R, Sethian JA, Vemuri BC.  Shape Modeling with Front Propagation: A Level Set Approach.  IEEE Transactions on Pattern Analysis and Machine Intelligence. 1995;17(2):158-75.
  * [15] Suri JS, Liu K, Singh S, Laxminarayan SN, Zeng X, Reden L.  Shape recovery algorithms using level sets in 2-D/3-D medical imagery: a state-of-the-art review.  IEEE Trans Inf Technol Biomed. 2002;6(1):8-28.
  * [16] Gomes J, Faugeras O.  Reconciling distance functions and Level Sets.  J Visual Communic and Imag Representation. 2000;11:209-23.
  * [17] Mumford D, Shaw J.  Optimal approximations by piecewise smooth functions and associated variational problems.  Communications on Pure and Applied Mathematics. 1989;42:577-685.
  * [18] Zhao HK, Chan T, Merriman B, Osher S.  A Variational Level Set Approach to Multiphase Motion.  Journal of Computational Physics. 1996;127(1):179-95.
  * [19] Vemuri B, Chen Y.  Joint image registration and segmentation.  In: Geometric Level Set Methods in Imaging, Vision and Graphics. Springer; 2003. p. 251-69.
  * [20] Liu C, Dong F, Zhu S.  New Variational Formulations for Level Set Evolution without Reinitialization with Applications to Image Segmentation.  J Math Imaging Vis. 2011;41:194-209.
  * [21] Li C. ISICDM; 2010 (accessed October 2020).  http://www.imagecomputing.org/~cmli.
  * [22] Pikachu C. Classical Image Processing Algorithms: Level Set. GitHub; 2019 (accessed October 2020).  https://github.com/CharlesPikachu/Classical-Image-Processing-Algorithms/tree/master/LevelSet.
  * [23] Tsankashvili N. Comparing Edge Detection Methods. A Medium Corporation; 2018.  Available from: https://medium.com/@nikatsanka/comparing-edge-detection-methods-638a2919476e.
  * [24] Wang LP, Titov A, McGibbon R, Liu F, Pande V, Martinez TJ.  Discovering chemistry with ab initio nanoreactor.  Nature Chemistry. 2014;6 (12):1044-8.
  * [25] Thompson CM, Shure L. Image Processing Toolbox: For Use with MATLAB;[user’s Guide]. MathWorks; 1995.
  * [26] Lam L, Lee SW.  Thinning Methodologies – A Comprehensive Survey.  IEEE Transactions of Pattern Analysis and Machine Intelligence. 1992;14(9):869-83.
  * [27] Zhang TY, Suen CY.  A Fast Parallel Algorithms for Thinning Digital Patterns.  Communications of the ACM. 1984;27(3):236-9.