%% sample.tex VERSION 2 (February 9th, 2007 - revised February 12th, 2007)
%% sampledocument.tex VERSION 3 (April 27th, 2007)

%% A sample full paper template you modify and save as fullpaperdocument.tex
%% to produce your camera-ready full paper

%% +-------------+
%% | PAPER TITLE |
%% +-------------+

\title[Partitioned IFT]{A partitioned algorithm for the image \\ foresting transform}


%% +---------+
%% | AUTHORS |
%% +---------+


\begin{Authors}
\Author{Felipe P.G. Bergo}
  \Affil[IC]{LIV, Institute of Computing, UNICAMP, Campinas, SP, Brazil \\*
  \email{bergo@liv.ic.unicamp.br, afalcao@ic.unicamp.br}}
\Author{Alexandre X. Falc{\~{a}}o}
  \Affilref[IC]
\end{Authors}


%% +-----------------------+
%% | ABSTRACT and KEYWORDS |
%% +-----------------------+
%% Abstract and keywords are written inside two environments,
%% the "Abstract" amd the "Keywords" environments, respectively.
%% Note the uppercase "A" and "K".
%% Keywords are separated with comma. The last keyword ends
%% with a period.

\begin{Abstract}
The Image Foresting Transform (IFT) is a powerful graph-based
framework for the design and implementation of image processing
operators. In this work we present the Partitioned IFT (PIFT), an
algorithm that computes any IFT operator as a series of independent
IFT-like computations. The PIFT makes parallelization of existing IFT
operators easy, and allows the computation of IFTs in systems with
scarce memory. We evaluate the PIFT for two image processing
applications: watershed segmentation and Euclidean distance
transforms.
\end{Abstract}

\begin{Keywords}
  \Index{graph algorithms},
  \Index{parallel algorithms},
  \Index{mathematical morphology},
  \Index{image foresting transform},
  \Index{distance transforms}.
\end{Keywords}

%% +----------+
%% | SECTIONS |
%% +----------+

\newcommand{\fns}{\footnotesize}
\newcommand{\ssz}{\scriptsize}
\newcommand{\M}{$\times 10^6$}
\newcommand{\CA}{{\cal A}}
\newcommand{\CE}{{\cal E}}
\newcommand{\CF}{{\cal F}}
\newcommand{\CM}{{\cal M}}
\newcommand{\CS}{{\cal S}}
\newcommand{\CR}{{\CS''}}
\newcommand{\CLS}{L_0}
\newcommand{\BI}{{\bf I}}
\newcommand{\oo}{\infty}
\newcommand{\nil}{{\it nil}}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\seq}[1]{\langle #1 \rangle}
\newcommand{\Ppath}{P^\ast}

\section{Introduction}

The Image Foresting Transform (IFT)~\cite{Bergo:Falcao04a} is a graph-based
framework for the design and implementation of image processing
operators. It reduces image processing operations, such as
watersheds~\cites{Bergo:Vincent91,Bergo:Beucher93}, morphological
reconstructions~\cite{Bergo:Falcao01}, skeletonization~\cite{Bergo:Falcao02} and
distance transforms~\cite{Bergo:Danielsson80}, to the computation of a
minimum-cost path forest over an implicit graph representation of the
image. The IFT runs in linear time, but it does not take advantage of
parallel and distributed computer systems. Its data structures also
require considerable memory space~\cite{Bergo:Felkel01}, and this can be a
limitation to the processing of large 3D images.

In this work we present the Partitioned IFT, an algorithm that
computes any IFT as a set of independent IFTs over partitions of the
input image. Both time and memory required to compute the IFT of each
partition are proportional to the size of that partition.  The
minimum-cost path forests of the partitions are merged by fast
differential IFTs~\cite{Bergo:Falcao04b}. This scheme provides the means to
take advantage of parallel and distributed computer systems (by
assigning each partition's IFT to a different central processing unit
(CPU)) and to allow the computation of IFTs with a reduced memory
footprint (by computing partition forests sequentially).

\section{Related works}

\subsection{Related algorithms}

Moga et al.~\cite{Bergo:Moga1998a} presented two parallel watershed
algorithms that treat the image as a graph and perform independent
flooding simulations in image partitions. Parallel flooding
simulations are repeated while plateaus overflow to adjacent
partitions. The same group~\cite{Bergo:Moga1998b} presented a similar
parallel algorithm for the computation of the watershed-from-markers
transform.  Both works achieve scalable speedups in parallel
architectures, but the speedup factor does not scale linearly with the
number of processors. Moga et al.~\cite{Bergo:Moga1998a} achieve speedup
factors\footnote{The speedup factor of a parallel algorithm on an
$n$-CPU parallel system is calculated as $\frac{t_1}{t_N}$, where
$t_1$ is the time required to perform the computation on a single-CPU
system, and $t_N$ is the time required to perform the computation on
an $n$-way system.} around $2$ for 4-CPU systems, and $3.5$ for 8-CPU
systems.  Bruno and Costa~\cite{Bergo:Bruno04} present a distributed
algorithm for the computation of Euclidean distance transforms (EDT)
based on morphological dilations. Their algorithm achieves a speedup
factor of 3.5 on a 4-CPU system.

\subsection{The image foresting transform}

The IFT algorithm is essentially Dijkstra's
algorithm~\cite{Bergo:Ahuja93}, modified for multiple sources and general
path cost functions~\cite{Bergo:Falcao04a}. The image is interpreted as a
directed graph whose nodes are the pixels. The edges are defined
implicitly by an \emph{adjacency relation} $\CA$. Tree roots are drawn
from a set $\CS$ of \emph{seed nodes} and path costs are given by a
\emph{path cost function} $f$. We use $\Ppath(s)$ to denote the
current path reaching pixel $s$, $\seq{s}$ to denote a \emph{trivial
path} containing a single node, and $\seq{s,t}$ to denote the edge
from pixel $s$ to pixel $t$. $\Ppath(s)\cdot\seq{s,t}$ is the path
that results from the concatenation of $\Ppath(s)$ and an edge
$\seq{s,t}$.

The choice of $\CA$, $\CS$ and $f$ define an IFT operator. The IFT
algorithm can compute by ordered propagation any forest property that
uses the seed set as reference. Usually, the IFT computes 4 maps: the
cost map $C$ stores the cost of the optimal path that reaches each
pixel, the predecessor map $P$ stores the predecessor of each pixel in
the forest, the root map $R$ stores the root of each pixel's optimal
path, and the label map $L$ stores object labels for each pixel.
Algorithm~\ref{\fullpaperid:a.ift} below computes the IFT.

\begin{nicealgo}{\fullpaperid:a.ift}
\naTITLE{\textbf{IFT}}
\naPREAMBLE
\naINPUT{Image $\BI$, Path-cost function $f$, Adjacency relation $\CA$, Seed set $\CS$ and 
Seed label map $\CLS$.}
\naOUTPUT{Cost map $C$, Predecessor map $P$, Root map $R$ and Label map $L$.}
\naAUX{Priority queue $Q$.}
\naBODY
\na{Set $Q\mget\emptyset$.}
\naBEGIN{\naFOREACH pixel $s \in \BI \setminus \CS$, \naDO}
\naEND{Set $C(s)\mget\oo$, $P(s)\mget \nil$, $R(s)\mget s$ and $L(s)\mget \nil$.}
\naBEGIN{\naFOREACH pixel $s \in \CS$, \naDO}
\na{Set $C(s)\mget f(\seq{s})$ and $L(s)\mget\CLS(s)$.}
\naEND{Insert $s$ in $Q$.}
\naBEGIN{\naWHILE $Q \neq \emptyset$, \naDO}
\na{Remove a pixel $s$ from $Q$ such that $C(s)$ is minimum.}
\naBEGIN{\naFOREACH $t$ such that $(s,t) \in \CA$, \naDO}
\na{Compute $cost\mget f(\Ppath(s)\cdot\seq{s,t})$.}
\naBEGIN{\naIF $cost < C(t)$ \naTHEN}
\na{\naIF $t \in Q$ \naTHEN remove $t$ from $Q$.}
\na{Set $P(t) \mget s$, $C(t) \mget cost$, $L(t) \mget L(s)$, $R(t) \mget R(s)$.}
\naENDN{3}{Insert $t$ in $Q$.}
\end{nicealgo}

Lines 1--3 set the forest to an initial state where every node's
optimum path is a trivial path with infinite cost. Lines 4--6 insert
the seed pixels in the priority queue with a trivial path cost
computed by $f$, and initialize seed labels for ordered
propagation. The loop of lines 7--14 uses the priority queue to
propagate the optimum paths and conquer the entire image. As long as
$f$ is finite and smooth~\cite{Bergo:Falcao04a}, an optimum path with finite
cost will be assigned to all pixels connected to $\CS$. Once a pixel
is removed from the queue (line 8), it is never inserted
again. Therefore, the main loop is repeated $|\BI|$ times. For integer
path costs with limited increments, $Q$ can be efficiently implemented
such that insertions and removals take $O(1)$
time~\cite{Bergo:Ahuja93}. With small adjacency relations ($|\CA| \ll
|\BI|$) and $O(1)$ queue operations, the IFT algorithm runs in
$O(|\BI|)$ time~\cite{Bergo:Falcao04a}.

%The memory required for the maps $C$, $P$, $R$, and $L$ and the
%priority queue $Q$ is considerable. Pixel addresses require a 32-bit
%representation: 16-bit addresses would restrict images to $256\times
%256$. Labels must be able to enumerate all interest objects. The
%object concept depends on the application, and usual labels
%representations need from 8 to 32 bits. The cost domain depends on the
%path cost function $f$, and also range from 8 to 32 bits for most
%applications. Dial's priority queue requires the storage of
%$2|\BI|+2K$ pixel addresses~\cite{Bergo:Falcao98}, where $K$ is the largest
%path cost increment -- the maximum value of
%$f(\Ppath(s)\cdot\seq{s,t}) - f(\Ppath(s))$. For most applications, $K
%\ll |\BI|$. In the worst scenario, the IFT requires about $24|\BI|$
%bytes of memory ($24 = 4 \mbox{ bytes} \times (|C|+|P|+|R|+|L|) + 8
%\mbox{ bytes} \times |I|$). The path cost function $f$ may require
%additional data, such as the image intensities. The memory requirement
%of the IFT is an order of magnitude above the space needed to store
%the input image. This can be a serious limitation to the processing of
%large 3D images and to the use of the IFT on portable and embedded
%devices with reduced memory.

Two common path cost functions for IFT operators are $f_{max}$ and
$f_{euc}$, shown in Eqs.~\ref{\fullpaperid:eq.fmax}--\ref{\fullpaperid:eq.feuc} below. Both
$f_{max}$ and $f_{euc}$ are \emph{smooth}, as required to ensure the
correctness of the IFT~\cite{Bergo:Falcao04a}.

\begin{equation}
f_{max}(\seq{s_1,\ldots,s_n}) = \left\{ \begin{array}{ll} max_{i=1}^{n}\left(I\left(s_i\right)\right) & \mbox{ if } n>1 \mbox{,} \\ h(s_1) & \mbox{otherwise.} \\ \end{array} \right.
\label{\fullpaperid:eq.fmax}
\end{equation}

\begin{equation}
f_{euc}(\seq{s_1,\ldots,s_n}) = \mbox{Euclidean distance between } s_1 \mbox{ and } s_n
\label{\fullpaperid:eq.feuc}
\end{equation}

where $I(s)$ is some value associated to pixel $s$ (such as intensity
or gradient intensity) and $h$ is a handicap function for trivial
paths. A watershed-from-markers transform can be implemented as an IFT
where $f$ is $f_{max}$ (Eq.~\ref{\fullpaperid:eq.fmax}), $h=0$ (for marker
imposition), $\CA$ is an adjacency with radius between $1$ and
$\sqrt{2}$ and $\CS$ contains the watershed
markers~\cites{Bergo:Falcao04a,Bergo:Falcao04b}. A classical watershed can be
implemented using $f=f_{max}$, $h(s)=I(s)+1$ and
$\CS=\BI$~\cite{Bergo:Lotufo00}. Function $f_{euc}$ (Eq.~\ref{\fullpaperid:eq.feuc})
allows the computation of distance transforms~\cite{Bergo:Danielsson80},
discrete Voronoi diagrams, skeletonizations and shape
saliences~\cites{Bergo:Falcao02,Bergo:Torres04,Bergo:Falcao04a,Bergo:Bergo06}.

\subsection{The differential image foresting transform}

The differential IFT~\cite{Bergo:Falcao04b} (DIFT) was motivated by
interactive 3D image segmentation applications where the user
interactively selects the seed pixels. It is quite common for the user
to add new seeds and remove previous ones based on the visualization
of the segmentation result. The first IFT is computed by
Algorithm~\ref{\fullpaperid:a.ift} as usual, from a seed set $\CS_0$.  The maps
$C$, $P$, $R$ and $L$ must be initialized to a forest of trivial paths
with infinite costs before the first DIFT is computed. Given a set
$\CS'$ of seeds to be added and a set $\CS''$ of tree roots to be
removed, the DIFT computes the optimum path forest for the effective
seed seet $\CS_1 = (\CS_0 \setminus \CS'') \cup \CS'$. The DIFT
processes only pixels affected by the seed set editing, and runs in
sublinear time. Instead of providing $\CS''$ directly, the DIFT takes
a set $\CM$ of removal markers, and $\CS''$ is computed as the set of
roots of the pixels in $\CM$. Algorithm~\ref{\fullpaperid:a.dift} below is the main
DIFT algorithm. The DIFT-TreeRemoval subroutine referenced in line 2
visits all pixels that belong to removed trees, sets their optimum
paths to trivial paths with infinite costs (forcing their
recalculation by Algorithm~\ref{\fullpaperid:a.dift}), and builds the set $\CF$ of
frontier pixels.

\begin{nicealgo}{\fullpaperid:a.dift}
\naTITLE{\textbf{DIFT}}
\naPREAMBLE
\naINPUT{Image $\BI$, Cost map $C$, Predecessor map $P$, Root map $R$, 
Label map $L$, Path-cost function $f$, Adjacency relation $\CA$, 
Set $\CS'$ of new seed pixels, Set $\CM$ of marking pixels, Seed label map $\CLS$.}
\naOUTPUT{$C$, $P$, $R$ and $L$.}
\naAUX{Priority queue $Q$, Frontier set $\CF$.}
\naBODY
\na{Set $Q\mget\emptyset$.}
\na{$(C, P, \CF) \mget $\naname{DIFT-TreeRemoval}$(C, P, R, L, \CA, \CM)$.}
\na{$\CF \mget \CF \setminus \CS'$.}
\naBEGIN{\naWHILE $\CS' \neq \emptyset$, \naDO}
\na{Remove any $t$ from $\CS'$.}
\naBEGIN{\naIF $f(\seq{t}) < C(t)$ \naTHEN}
\na{Set $C(t) \mget f(\seq{t})$, $R(t) \mget t$, $L(t)\mget \CLS(t)$, $P(t) \mget \nil$.}
\naENDN{2}{Set $\CF \mget \CF \cup \{ t \}$.}
\naBEGIN{\naWHILE $\CF \neq \emptyset$, \naDO}
\naEND{Remove any $t$ from $\CF$ and insert $t$ in $Q$.}
\naBEGIN{\naWHILE $Q \neq \emptyset$, \naDO}
\na{Remove a pixel $s$ from $Q$, such that $C(s)$ is minimum.}
\naBEGIN{\naFOREACH $t$ such that $(s,t) \in \CA$, \naDO}
\na{Compute $cost \mget f(\Ppath(s)\cdot\seq{s,t})$.}
\naBEGIN{\naIF $cost < C(t)$ or $P(t)=s$ \naTHEN}
\na{\naIF $t \in Q$ \naTHEN remove $t$ from $Q$.}
\na{Set $P(t) \mget s$, $C(t) \mget cost$, $R(t) \mget R(s)$, $L(t) \mget L(s)$.}
\naENDN{3}{Insert $t$ in $Q$.}
\end{nicealgo}

Lines 2--3 compute a set $\CF$ of frontier pixels that belong to
non-removed trees but share edges with pixels in removed trees. Lines
4--10 insert the new seeds and the frontier pixels in the queue. Lines
11--18 are very much like the main loop of the IFT Algorithm
(Algorithm~\ref{\fullpaperid:a.ift}), except for the condition $P(t)=s$ in line 14,
which forces the update of all pixels that had their optimum paths
modified. The result of the DIFT is an optimum path forest for the
``effective seed set'' $\CS_1 = (\CS_0 \setminus \CS'') \cup \CS'$.

\section{The partitioned image foresting transform}

In the Partitioned IFT (PIFT), we split the input image and seed set
in $N_P$ \emph{partitions}. The number of partitions can be chosen to
match the number of available processing nodes, or so that the
computer system has enough memory to run the IFT algorithm on each
image partition. Partitions do not need to be equally sized. We
compute independent IFTs on each partition. At this point, we have an
optimum forest that ignores the inter-partition edges of the
graph. Figure~\ref{\fullpaperid:f.pift1}a shows an example of this
partial result for the EDT using a set of random pixels as seeds and 3
partitions. To allow propagation through the inter-partition graph
edges, we consider the paths obtained by the concatenation of each
edge $\seq{s,t}$ to $\Ppath(s)$
(Figure~\ref{\fullpaperid:f.pift1}c). When $f(\Ppath(s)\cdot\seq{s,t})$
is less than the current cost of $t$, or the edge was part of the
destination pixel's previous optimal path, the endpoint is added as
seed in a differential IFT so that it can be propagated. If more than
one inter-partition edge share a same endpoint $t$, the one that
provides the lower path cost $\Ppath(t)$ is propagated.  A new
iteration of differential IFTs is computed for each partition.  The
PIFT halts when no inter-partition edge satisfies the criteria for
addition. Figure~\ref{\fullpaperid:f.pift1}b shows the complete EDT,
obtained after 2 iterations over the 3 partitions.

\begin{figure}[!htb]\begin{center}\begin{tabular}{ccc}
\includegraphics[width=3.5cm]{\fullpaperpdirectory/images/fig1a.png} &
\includegraphics[width=3.5cm]{\fullpaperpdirectory/images/fig1b.png} &
\includegraphics[width=3.5cm]{\fullpaperpdirectory/images/fig0.png} \\ 
(a) & (b) & (c) \\
\end{tabular}\end{center}
\caption{Labels of an EDT with the Partitioned IFT: (a) Partial
result after the first iteration and (b) final result after the second
iteration. (c) PIFT notation: $\seq{s,t}$ is an inter-partition edge,
$\Ppath(s)$ is the optimum path assigned to $s$, and $R(s)$ the root
of $\Ppath(s)$.}\vspace{-5pt}
\label{\fullpaperid:f.pift1}
\end{figure}

The differential IFTs used in the Partitioned IFT always have an empty
set of removal markers. The Partition-IFT algorithm below
(Algorithm~\ref{\fullpaperid:a.pift}) computes the IFT within a partition. It is
essentially the differential IFT algorithm without tree removal, and
with special treatment of inter-partition edges.

\begin{nicealgo}{\fullpaperid:a.pift}
\naTITLE{\textbf{Partition-IFT}}
\naPREAMBLE
\naINPUT{Image partition $\BI'$, Cost map $C$, Predecessor map $P$, 
Root map $R$, Label map $L$, Path-cost function $f$, Adjacency relation $\CA$, 
Set $\CS$ of seed pixels, Seed label map $\CLS$, Set $\CE_I$ of incoming
inter-partition edges.}
\naOUTPUT{Maps $C$, $P$, $R$, $L$ and Set $\CE_O$ of outgoing inter-partition
edges.}
\naAUX{Priority queue $Q$.}
\naBODY
\na{Set $Q\mget\emptyset$, $\CE_O\mget\emptyset$.}
%% process \CS
\naBEGIN{\naIF $\CS \neq \emptyset$ \naTHEN}
\naBEGIN{\naFOREACH pixel $s \in \BI'$, \naDO}
\naEND{Set $C(s)\mget\oo$, $P(s)\mget \nil$, $R(s)\mget s$ and $L(s)\mget \nil$.}
\naBEGIN{\naFOREACH pixel $s \in \CS$, \naDO}
\na{Set $C(s)\mget f(\seq{s})$ and $L(s)\mget\CLS(s)$.}
\naENDN{2}{Insert $s$ in $Q$.}
%% process \CE_I
\naBEGIN{\naFOREACH edge $\seq{s,t} \in \CE_I$, \naDO}
\na{Compute $cost \mget f(\Ppath(s)\cdot\seq{s,t})$.} 
\naBEGIN{\naIF $cost < C(t)$ or $P(t)=s$ \naTHEN}
\na{Set $C(t)\mget cost$, $P(t)\mget s$, $R(t)\mget R(s)$ and $L(t)\mget L(s)$.}
\naENDN{2}{Insert $t$ in $Q$.}
%% main loop
\naBEGIN{\naWHILE $Q \neq \emptyset$, \naDO}
\na{Remove a pixel $s$ from $Q$, such that $C(s)$ is minimum.}
\naBEGIN{\naFOREACH $t$ such that $(s,t) \in \CA$, \naDO}
\naBEGIN{\naIF $t \in \BI'$ \naTHEN}
\na{Compute $cost \mget f(\Ppath(s)\cdot\seq{s,t})$.}
\naBEGIN{\naIF $cost < C(t)$ or $P(t)=s$ \naTHEN}
\na{\naIF $t \in Q$ \naTHEN remove $t$ from $Q$.}
\na{Set $P(t) \mget s$, $C(t) \mget cost$, $R(t) \mget R(s)$, $L(t) \mget L(s)$.}
\naENDN{2}{Insert $t$ in $Q$.}
\naENDN{2}{\naELSE Insert $\seq{s,t}$ in $\CE_O$.}
\end{nicealgo}

The DIFT is unable to tell whether the algorithm is on the first
iteration, therefore the initial state of the forest must be set before
the first iteration. In the PIFT, the seed set $\CS$ will only be
non-empty in the first iteration. We use this property to initialize
the partition's forest to trivial paths with infinite costs in lines
2--4. Lines 5--7 queue and initialize the seed pixels in the same way
the IFT does. Lines 8--12 process the incoming inter-partition edges
$\CE_I$. Edges that offer lower costs to their endpoints or belonged
to the previous forest are queued for propagation. If multiple edges
in $\CE_I$ reach the same endpoint, the edge that provides the lower
cost for the endpoint takes precedence. %
%%Duplicated pixel queueing
%%from the loops of lines 5--7 and 8--12 never occurs: $\CE_I =
%%\emptyset$ in the first PIFT iteration, and $\CS = \emptyset$ in the
%%subsequent iterations.  Therefore, only one of these loops are
%%executed during a PIFT iteration. 
The main loop in lines 13--22 is very similar to the main loop of the
DIFT, with the addition of the partition test $t \in \BI'$ in line
16. Edges within the current partition are processed normally.
Inter-partition edges are added to the outgoing edge set $\CE_O$ (line
22).  Note that the cost computation in line 9 may require additional
information about $\Ppath(s)$, which can contain pixels of several
partitions. All path information required to compute
$f(\Ppath(s)\cdot\seq{s,t})$ must be passed along with the set
$\CE_I$. For $f_{max}$, only $C(s)$ is required. For $f_{euc}$, only
$R(s)$ is required. Since $L(s)$ may be propagated in line 11, it must
also be part of the input.  Passing each element of $\CE_I$ as
$\{s,t,C(s),R(s),L(s)\}$ is enough to compute the PIFT with either
$f_{max}$ or $f_{euc}$. The PIFT algorithm (Algorithm~\ref{\fullpaperid:a.pift2}) that
computes the IFT of an image $\BI$ from its partitions is shown below.

\begin{nicealgo}{\fullpaperid:a.pift2}
\naTITLE{\textbf{Partitioned IFT}}
\naPREAMBLE
\naINPUT{Image $\BI$, Path-cost function $f$, Adjacency relation $\CA$, 
Set $\CS$ of seed pixels, Seed label map $\CLS$, Number of partitions $N_P$.}
\naOUTPUT{Cost map $C$, Predecessor map $P$, Root map $R$, Label map $L$.}
\naAUX{Edge sets $\CE$, $\CE'$, $\CE''$ and $\CE'''$, Seed set $\CS'$.}
\naBODY
\na{Set $\CE\mget\emptyset$.}
\na{Split $\BI$ in $N_P$ partitions $\BI[1]\ldots\BI[N_P]$.}
\naBEGIN{\naFOR $i=1$ \naTO $N_P$, \naDO}
\na{Set $\CS'=\{ s \;|\; s\in \CS \land s \in \BI[i] \}$.}
\na{Set $(C[i],P[i],R[i],L[i],\CE') \mget$}
\nac{\naname{Partition-IFT}$(\BI[i],C[i],P[i],R[i],L[i],f,\CA,\CS',\CLS,\emptyset)$.}
\naEND{Set $\CE \mget \CE \cup \CE'$.}
\naBEGIN{\naREPEAT}
\na{Set $\CE'''\mget\emptyset$.}
\naBEGIN{\naFOR $i=1$ \naTO $N_P$, \naDO}
\na{Set $\CE''=\{ \seq{s,t} \;|\; \seq{s,t} \in \CE \land t \in \BI[i] \}$.}
\na{Set $(C[i],P[i],R[i],L[i],\CE') \mget$}
\nac{\naname{Partition-IFT}$(\BI[i],C[i],P[i],R[i],L[i],f,\CA,\emptyset,\nil,\CE'')$.}
\naEND{Set $\CE''' \mget \CE''' \cup \CE'$.}
\na{Set $\CE\mget\CE'''$.}
\naEND{\naUNTIL $\CE = \emptyset$.}
\na{Set $C\mget\cup_{i=1}^{N_P}C[i]$, $P\mget\cup_{i=1}^{N_P}P[i]$, $R\mget\cup_{i=1}^{N_P}R[i]$ and $L\mget\cup_{i=1}^{N_P}L[i]$.} 
\end{nicealgo}

Lines 1--2 initialize the inter-partition edge set $\CE$ and split the
input image in $N_P$ partitions. The loop in lines 3--6 run the first
IFT iteration on each partition. All inter-partition edges are
accumulated in the set $\CE$. The loop in lines 7--14 run the
remaining IFT iterations on the partitions, until no propagation
occurs and the set $\CE$ of inter-partition edges is empty (line 14).

For parallel architectures, both loops (lines 3--6 and 7--14) can be
done in parallel. For distributed systems, the executions of
Partition-IFT (Algorithm~\ref{\fullpaperid:a.pift}) can be performed as remote procedure
calls. Note that the partitioned maps ($C[i]$, $P[i]$, $R[i]$ and
$L[i]$) are only needed at the end of the algorithm, to compose the
final IFT maps. In a distributed implementation, these maps can be
kept on the remote processing nodes and do not need to be transferred
at each call to Partition-IFT, as they are not modified by the
caller.

\paragraph{Performance Considerations.} The overall number of pixels
processed by the PIFT is larger than $|\BI|$. After the loop of lines
3--6 of Algorithm~\ref{\fullpaperid:a.pift2}, the PIFT has already
processed $|\BI|$ nodes. However, the number of pixels processed by
the loop of lines 7--14 decreases at each iteration, and the algorithm
converges rapidly to the optimum path forest. The number of PIFT
iterations --- i.e., one iteration of the loop of lines 3--6 plus the
number of iterations of the loop of lines 7--14 --- is bounded by the
maximum number of inter-partition edges contained by an optimum path,
plus one. Each inter-partition edge postpones the resolution of the
optimum path to the next PIFT
iteration. Figure~\ref{\fullpaperid:f.bound} illustrates some
examples. In an Euclidean distance transform
(Figure~\ref{\fullpaperid:f.bound}a), all paths flow away from the
roots, and a path may cross at most $N_P - 1$ partition boundaries,
requiring at most $N_P$ PIFT iterations. For path cost functions like
$f_{max}$, there is no restriction to the shape of optimum paths, and
cases like the one in Figure~\ref{\fullpaperid:f.bound}b can
occur. However, as the number of iterations increases, the number of
pixels processed by each iteration tends to decrease, and the PIFT
converges more rapidly to the optimum forest.

%%The experiments in the next section show that the
%%total number of pixels processed does not exceed $|\BI|\times N_P$
%%even in cases with a large number of iterations, and the PIFT offers
%%real performance advantage for parallel architectures, as long as
%% inter-iteration communication .

\begin{figure}[!htb]\begin{center}\begin{tabular}{cc}
\includegraphics[width=3.5cm]{\fullpaperpdirectory/images/fig2a.png} &
\includegraphics[width=3.5cm]{\fullpaperpdirectory/images/fig2b.png} \\
(a) & (b) \\
\end{tabular}\end{center}
\caption{Partition crossings and PIFT iterations: In the PIFT-EDT,
paths cross at most $N_P - 1$ partition boundaries. In (a),
$\Ppath(p)$ crosses 2 boundaries to reach $p$ from $a$. The
numbers are the iteration in which the path segment is propagated. (b)
For general path-cost functions, a path may cross partition boundaries
several times.}\vspace{-5pt}
\label{\fullpaperid:f.bound}
\end{figure}

\section{Experimental results}

We implemented the PIFT as a client-server system, with a simple TCP
stream-based protocol for communication between the master
client that executes Algorithm~\ref{\fullpaperid:a.pift2} and the distributed
servers that execute Algorithm~\ref{\fullpaperid:a.pift}. In our implementation,
the image is always split in equal-sized partitions, using the $x$
coordinate to separate partitions (such as in Figure~\ref{\fullpaperid:f.pift1}a). We
chose 3 applications to evaluate the PIFT:

\begin{enumerate}
\item WS-BRAIN: Watershed-based segmentation of a 3D MR image of the
brain, using $f=f_{max}$ and $\CA=$6-neighborhood adjacency. Seeds
were selected interactively in the background and in the brain. The
gradient intensity was computed by a Gaussian enhacement filter
followed by morphological gradient computation~\cite{Bergo:Falcao04b}. The
size of the image is $356\times 356\times 241$, with a voxel size of
$0.70mm^3$ (Figure~\ref{\fullpaperid:f.app1}a--c).

\item EDT-RND: Euclidean distance transform of 1000 random points within a
$256^3$ volume, using $f=f_{euc}$ and $\CA=$26-neighborhood (Figure~\ref{\fullpaperid:f.app1}d).

\item EDT-BRAIN: Euclidean distance transform using the border of the
brain object (segmented in the first application) as seed
set ($|\CS|=355,556$). $f=f_{euc}$, $\CA=$26-neighborhood and volume size is $356\times
356\times 241$ (Figure~\ref{\fullpaperid:f.app1}e).
\end{enumerate}

\begin{figure}[!htb]\begin{center}\begin{tabular}{ccccc}
\includegraphics[height=1.97cm]{\fullpaperpdirectory/images/fig3a.png} &
\includegraphics[height=1.97cm]{\fullpaperpdirectory/images/fig3b.png} &
\includegraphics[height=1.97cm]{\fullpaperpdirectory/images/fig3c.png} &
\includegraphics[height=1.97cm]{\fullpaperpdirectory/images/fig3d.png} &
\includegraphics[height=1.97cm]{\fullpaperpdirectory/images/fig3e.png} \\
(a) & (b) & (c) & (d) & (e)
\end{tabular}\end{center}
\caption{Images from the evaluation applications: (a) Slice from the WS-BRAIN input image. 
(b) gradient intensity of (a). (c) 3D renderization of the WS-BRAIN result.
(d) Visualization of the discrete Voronoi diagram, result of the EDT-RND.
(e) Slice from the distance map computed in EDT-BRAIN.}\vspace{-5pt}
\label{\fullpaperid:f.app1}
\end{figure}

First, we measured the processing overhead of the PIFT as the number
of partitions ($N_P$) increases. We computed the 3 applications with
the PIFT, using from 1 to 10 partitions. Table~\ref{\fullpaperid:t.exp1} and
Figure~\ref{\fullpaperid:f.plot1} present the number of nodes processed in each case
and the upper bound for the speedup factor. These results indicate
that a 10-way parallel system may be able to offer a speedup factor of
6.60 to the EDT computation, and a factor of 2.34 to the Watershed
transform on these instances of problems.

The EDT computations required at most 4 iterations before
halting. PIFTs based on $f_{max}$ are less efficient, since they allow
free-form paths that can traverse several partitions. This can be
noticed by the irregularity and increased slope of the plot in
Figure~\ref{\fullpaperid:f.plot1}b, as compared to
Figure~\ref{\fullpaperid:f.plot1}a. The WS-BRAIN PIFTs required at
most 23 iterations to converge. The number of processed nodes grows
linearly with the number of partitions.  In real data with non-uniform
distributions (WS-BRAIN and EDT-BRAIN), bad choices of partition
boundaries may increase the number of processed nodes, such as in the
$N_P=4$ and $N_P=8$ cases of EDT-BRAIN and $N_P=5$ of WS-BRAIN.

\begin{table}[!htb]
\caption{Number of processed nodes and upper bound for the speedup
factor in each application, using up to 10 partitions.}\vspace{-5pt}
\label{\fullpaperid:t.exp1}
\begin{center}
\begin{tabular}{|c|rr|rr|rr|}
\hline
\multirow{2}{*}{\bf\fns $N_P$} & \multicolumn{2}{c|}{\bf\fns WS-BRAIN} & \multicolumn{2}{c|}{\bf\fns EDT-RND} & \multicolumn{2}{c|}{\bf\fns EDT-BRAIN} \\
& {\ssz Nodes} & {\ssz Speedup} & {\ssz Nodes} & {\ssz Speedup} & {\ssz Nodes} & {\ssz Speedup} \\
\hline
1  & {\fns  30.5\M} & {\fns 1.00} & {\fns 16.8\M} & {\fns 1.00} & {\fns 30.5\M} & {\fns 1.00} \\
2  & {\fns  48.6\M} & {\fns 1.25} & {\fns 17.3\M} & {\fns 1.94} & {\fns 31.5\M} & {\fns 1.93} \\
3  & {\fns  59.0\M} & {\fns 1.55} & {\fns 17.7\M} & {\fns 2.84} & {\fns 34.2\M} & {\fns 2.67} \\
4  & {\fns  62.7\M} & {\fns 1.94} & {\fns 18.2\M} & {\fns 3.69} & {\fns 39.6\M} & {\fns 3.08} \\
5  & {\fns  76.7\M} & {\fns 1.98} & {\fns 18.6\M} & {\fns 4.51} & {\fns 37.8\M} & {\fns 4.03} \\
6  & {\fns  75.1\M} & {\fns 2.43} & {\fns 19.2\M} & {\fns 5.25} & {\fns 38.7\M} & {\fns 4.72} \\
7  & {\fns  92.2\M} & {\fns 2.31} & {\fns 19.7\M} & {\fns 5.96} & {\fns 42.8\M} & {\fns 4.98} \\
8  & {\fns  98.1\M} & {\fns 2.48} & {\fns 20.1\M} & {\fns 6.68} & {\fns 46.8\M} & {\fns 5.21} \\
9  & {\fns 106.5\M} & {\fns 2.57} & {\fns 20.5\M} & {\fns 7.37} & {\fns 42.2\M} & {\fns 6.50} \\
10 & {\fns 130.2\M} & {\fns 2.34} & {\fns 21.0\M} & {\fns 8.00} & {\fns 46.2\M} & {\fns 6.60} \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[!htb]\begin{center}\begin{tabular}{cc}
\includegraphics[height=3.8cm]{\fullpaperpdirectory/images/plot1.png} &
\includegraphics[height=3.8cm]{\fullpaperpdirectory/images/plot2.png} \\
(a) & (b) \\
\end{tabular}\end{center}
\caption{Number of processed nodes vs. number of partitions for
(a) EDT-RND, EDT-BRAIN and (b) WS-BRAIN.}\vspace{-5pt}
\label{\fullpaperid:f.plot1}
\end{figure}

In a second set of experiments we used the PIFT to compute EDT-RND,
EDT-BRAIN and WS-BRAIN in two parallel systems: a PC with 2 CPUs
(Athlon MP 1800+@1150~MHz) and 2~GB of RAM, and a Compaq AlphaServer
GS140 6/525 with 10 CPUs (Alpha EV6@525 MHz) and 8~GB of
RAM. Table~\ref{\fullpaperid:t.exp2} presents the results. On the EDT applications,
we achieved speedup factors very close to the measured upper bounds
(Table~\ref{\fullpaperid:t.exp1}) for $N_P=2$ and $N_P=4$. On other hand, there was
little or no speedup for the watershed application. Our prototype
implementation uses a naive communication protocol with no data
compression. Besides that, the edge set transfers of lines 5 and 11 of
Algorithm~\ref{\fullpaperid:a.pift2} were implemented in a sequential way, and instances
with a large number of partitions and/or a large number of PIFT
iterations (such as WS-BRAIN with $N_P=10$) performed poorly because the
CPUs remained idle while waiting for the client to complete the
sequential edge set transfers.

\begin{table}
\caption{PIFT performance on two parallel computer systems. Times are given in seconds.}
\label{\fullpaperid:t.exp2}\vspace{-5pt}
\begin{center}
\begin{tabular}{|c|c|rr|rr|rr|}
\hline
\multirow{2}{*}{\bf\fns System} & \multirow{2}{*}{\bf\fns $N_P$} & \multicolumn{2}{c|}{\bf\fns WS-BRAIN} & \multicolumn{2}{c|}{\bf\fns EDT-RND} & \multicolumn{2}{c|}{\bf\fns EDT-BRAIN} \\
& & {\ssz Time} & {\ssz Speedup} & {\ssz Time} & {\ssz Speedup} & {\ssz Time} & {\ssz Speedup} \\
\hline
\multirow{2}{*}{{\fns Dual Athlon}} %
& {\fns 1} & {\fns 258.1} & {\fns 1.00} & {\fns 195.9} & {\fns 1.00} & {\fns 459.5} & {\fns 1.00} \\
& {\fns 2} & {\fns 242.9} & {\fns 1.06} & {\fns 106.2} & {\fns 1.84} & {\fns 246.3} & {\fns 1.87} \\
\hline
\multirow{5}{*}{{\fns 10-CPU GS140}} %
& {\fns 1} & {\fns 280.6} & {\fns 1.00} &  {\fns 228.8} & {\fns 1.00} & {\fns 611.4} & {\fns 1.00} \\
& {\fns 2} & {\fns 284.3} & {\fns 0.99} &  {\fns 126.6} & {\fns 1.81} & {\fns 324.2} & {\fns 1.89} \\
& {\fns 4} & {\fns 226.2} & {\fns 1.24} &  {\fns 73.0} & {\fns 3.13} & {\fns  274.3} & {\fns 2.23} \\
& {\fns 8} & {\fns 249.3} & {\fns 1.13} &  {\fns 49.3} & {\fns 4.64} & {\fns  214.1} & {\fns 2.86} \\
& {\fns 10} & {\fns 336.4} & {\fns 0.83} & {\fns 47.9} & {\fns 4.78} & {\fns  197.7} & {\fns 3.09} \\
\hline
\end{tabular}
\end{center}
\end{table}

\section{Conclusion and future works}

We introduced the Partitioned Image Foresting Transform, an algorithm
that computes minimum-cost path forests as a set of independent
DIFTs~\cites{Bergo:Falcao04a,Bergo:Falcao04b} in partitions of the input
image. The PIFT is useful for taking advantage of parallel computer
systems and for computing IFTs in computer systems with limited
memory, such as handhelds and embedded systems. The PIFT is applicable
to any IFT-based operator, and therefore can be readily employed to
parallelize morphological reconstructions~\cite{Bergo:Falcao01}, watershed
transforms~\cites{Bergo:Vincent91,Bergo:Beucher93,Bergo:Falcao04b,Bergo:Bergo06}, distance
transforms~\cites{Bergo:Danielsson80,Bergo:Falcao04a} and
skeletonizations~\cites{Bergo:Falcao02,Bergo:Torres04}, among other operators.  It
is a trend in microprocessor technology to compensate CPU speed
limitations by producing multi-core CPUs. The PIFT is an important
contribution that allows existing image processing applications to use
modern hardware efficiently with minimum effort.

We implemented a prototype PIFT system with a simple client-server
architecture built on top of TCP streams. Even with no data
compression and with some inneficient network operations, we achieved
speedup factors very close to the expected upper bounds for EDT
operations. PIFT-based watershed segmentation performed poorly due to
the inneficiency of edge set transfers in our prototype.  With a
better protocol, the PIFT should be able to reach speedup factors
closer to the upper bounds in Table~\ref{\fullpaperid:t.exp1}. 

Future works include: development of better protocols for
implementation of the PIFT in parallel systems, evaluation of the
speedup bounds for specific operators -- such as the watershed
transform -- and investigation of enhancements to the PIFT such as
partitioning schemes and iteration scheduling among nodes.

%
% BIBLIOGRAPHY
% for details, see ftp://ftp.ams.org/pub/tex/amsrefs/amsrdoc.pdf
%

\begin{bibsection}
	\begin{biblist}[\resetbiblist{9}]

\bib{Bergo:Ahuja93}{book}{
author = {Ahuja, R. K.},
author = {Magnanti, T. L.}
author = {Orlin, J. B.},
title = {Network Flows: Theory, Algorithms and Applications},
publisher = {Prentice-Hall},
date = {1993}
}

\bib{Bergo:Falcao04a}{article}{
author = {Falc{\~{a}}o, A. X.},
author = {Stolfi, J.},
author = {Lotufo, R. A.},
title  = {The Image Foresting Transform: Theory, Algorithms, and Applications},
journal = {IEEE Trans. on Pattern Analysis and Machine Intelligence},
volume = {26},
number = {1},
pages  = {19--29},
date   = {2004}
}

\bib{Bergo:Beucher93}{book}{
author = {Beucher, S.},
author = {Meyer, F.},
booktitle = {Mathematical Morphology in Image Processing},
title = {The morphological approach to segmentation: The watershed transformation},
publisher = {Marcel Dekker},
chapter   = {12},
pages = {433--481},
date = {1993}
}

\bib{Bergo:Vincent91}{article}{
author = {Vincent, L.},
author = {Soille, P.},
title  = {Watersheds in Digital Spaces: An Efficient Algorithm Based on Immersion Simulations},
journal = {IEEE Trans. on Pattern Analysis and Machine Intelligence},
volume = {13},
number = {6},
date = {1991-06}
}

\bib{Bergo:Bergo06}{article}{
author = {F. P. G. Bergo},
author = {A. X. Falc{\~{a}}o},
title = {Fast and automatic curvilinear reformatting of {MR} images of the brain for diagnosis of dysplastic lesions},
journal = {Proc. of 3rd Intl. Symp. on Biomedical Imaging},
publisher = {IEEE},
pages = {486--489},
date = {2006-04}
}

\bib{Bergo:Torres04}{article}{
author = {R. S. Torres},
author = {A. X. Falc{\~{a}}o},
author = {L. F. Costa},
title  = {A graph-based approach for multiscale shape analysis},
journal = {Pattern Recognition},
publisher  = {Elsevier},
pages      = {1163--1174},
volume     = {37},
number     = {6},
date       = {2004-06}
}

\bib{Bergo:Lotufo00}{article}{
 author    = {R. A. Lotufo},
 author    = {A. X. Falc{\~{a}}o},
 journal   = {Mathematical Morphology and its Applications to Image and Signal Processing},
 title     = {The ordered queue and the optimality of the watershed approaches},
 publisher = {Kluwer},
 volume    = {18},
 pages     = {341--350},
 date      = {2000-06}
}

\bib{Bergo:Falcao01}{article}{
 author    = {A. X. Falc{\~{a}}o},
 author    = {B. S. da Cunha},
 author    = {R. A. Lotufo},
 journal   = {Proc. of SPIE on Medical Imaging},
 title     = {Design of connected operators using the image foresting transform},
 volume    = {4322},
 pages     = {468--479},
 date      = {2001-02}
}

\bib{Bergo:Falcao02}{article}{
author     = {A. X. Falc{\~{a}}o},
author = {L. F. Costa},
author = {B. S. da Cunha},
title      = {Multiscale skeletons by image foresting transform and its applications to neuromorphometry},
journal    = {Pattern Recognition},
publisher  = {Elsevier},
pages      = {1571--1582},
volume     = {35},
number     = {7},
date       = {2002-04}
}

\bib{Bergo:Danielsson80}{article}{
  author = {P. E. Danielsson},
  title  = {Euclidean Distance Mapping},
  journal = {Computer Graphics and Image Processing},
  volume = {14},
  pages = {227--248},
  date = {1980}
}

\bib{Bergo:Felkel01}{article}{
author = {P. Felkel},
author = {M. Bruckschwaiger}
author = {R. Wegenkittl},
title = {Implementation and Complexity of the Watershed-from-Markers Algorithm Computed as a Minimal Cost Forest},
  number =       {3},
  volume =       {20},
  pages =        {C26--C35},
  date =         {2001},
  journal =    {Computer Graphics Forum (Eurographics)}
}

\bib{Bergo:Falcao04b}{article}{
author = {A. X. Falc{\~{a}}o},
author = {F. P. G. Bergo},
title     = {Interactive Volume Segmentation with Differential Image Foresting Transforms},
journal   = {IEEE Trans. on Medical Imaging},
volume    = {23},
number    = {9},
pages     = {1100--1108},
date      = {2004}
}


\bib{Bergo:Moga1998a}{article}{
author = {A. N. Moga},
author = {B. Cramariuc},
author = {M. Gabbouj},
title = {Parallel watershed transformation algorithms for image segmentation},
journal = {Parallel Computing},
volume = {24},
number = {14},
publisher = {Elsevier},
pages = {1981--2001},
date = {1998-12}
}

\bib{Bergo:Moga1998b}{article}{
author = {A. N. Moga},
author = {M. Gabbouj},
title = {Parallel Marker-Based Image Segmentation with Watershed Transformation},
journal = {Journal of Parallel and Distributed Computing},
volume = {51},
number = {1},
publisher = {Elsevier},
pages = {27--45},
date = {1998-05}
}

\bib{Bergo:Bruno04}{article}{
author = {O. M. Bruno},
author = {L. F. Costa},
title = {A parallel implementation of exact Euclidean distance transform based on exact dilations},
journal = {Microprocessors and Microsystems},
volume = {28},
number = {3},
date = {2004-04},
pages = {107--113},
publisher = {Elsevier}
}

\end{biblist}
\end{bibsection}


