%% 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 | %% +-------------+ \renewcommand{\shorttitle}{Partitioned IFT} %\title[Partitioned IFT]{A partitioned algorithm for the image \\ foresting transform} \title{A partitioned algorithm for the image foresting transform} %% +---------+ %% | AUTHORS | %% +---------+ \begin{Authors} \Author{Felipe P.\,G. Bergo} \Affil[IC]{Laboratório de Informática Visual (LIV), Instituto de Computação (IC), Universidade Estadual de Campinas (Unicamp), 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 | %% +----------+ \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{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 Equations~\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}$ (Equation~\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}$ (Equation~\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. \vspace{2pt} \begin{nicealgo}{\fullpaperid:a.dift} \naTITLE{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. \vspace{5pt} \begin{figure}[!htb]\begin{center}\begin{tabular}{ccc} \includegraphics[width=3.3cm]{\fullpaperpdirectory/images/fig1a.png} & \includegraphics[width=3.3cm]{\fullpaperpdirectory/images/fig1b.png} & \includegraphics[width=3.3cm]{\fullpaperpdirectory/images/fig0.png} \\ (a) & (b) & (c) \end{tabular}\end{center} \vspace{-9pt} \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{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. \vspace{2pt} \begin{nicealgo}{\fullpaperid:a.pift2} \naTITLE{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 . \vspace{5pt} \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} \vspace{-8pt} \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.95cm]{\fullpaperpdirectory/images/fig3a.png} & \includegraphics[height=1.95cm]{\fullpaperpdirectory/images/fig3b.png} & \includegraphics[height=1.95cm]{\fullpaperpdirectory/images/fig3c.png} & \includegraphics[height=1.95cm]{\fullpaperpdirectory/images/fig3d.png} & \includegraphics[height=1.95cm]{\fullpaperpdirectory/images/fig3e.png} \\ (a) & (b) & (c) & (d) & (e) \end{tabular}\end{center} \vspace{-8pt} \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} \renewcommand{\baselinestretch}{1.2}% for tabular environment \small \begin{tabular}{crcrcrc} \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} \\ \cline{2-7} & \multicolumn{1}{c}{\ssz Nodes} & {\ssz Speedup} & \multicolumn{1}{c}{\ssz Nodes} & {\ssz Speedup} & \multicolumn{1}{c}{\ssz Nodes} & {\ssz Speedup} \\ \hline 1 & {\fns 30.5\BergoM} & {\fns 1.00} & {\fns 16.8\BergoM} & {\fns 1.00} & {\fns 30.5\BergoM} & {\fns 1.00} \\ 2 & {\fns 48.6\BergoM} & {\fns 1.25} & {\fns 17.3\BergoM} & {\fns 1.94} & {\fns 31.5\BergoM} & {\fns 1.93} \\ 3 & {\fns 59.0\BergoM} & {\fns 1.55} & {\fns 17.7\BergoM} & {\fns 2.84} & {\fns 34.2\BergoM} & {\fns 2.67} \\ 4 & {\fns 62.7\BergoM} & {\fns 1.94} & {\fns 18.2\BergoM} & {\fns 3.69} & {\fns 39.6\BergoM} & {\fns 3.08} \\ 5 & {\fns 76.7\BergoM} & {\fns 1.98} & {\fns 18.6\BergoM} & {\fns 4.51} & {\fns 37.8\BergoM} & {\fns 4.03} \\ 6 & {\fns 75.1\BergoM} & {\fns 2.43} & {\fns 19.2\BergoM} & {\fns 5.25} & {\fns 38.7\BergoM} & {\fns 4.72} \\ 7 & {\fns 92.2\BergoM} & {\fns 2.31} & {\fns 19.7\BergoM} & {\fns 5.96} & {\fns 42.8\BergoM} & {\fns 4.98} \\ 8 & {\fns 98.1\BergoM} & {\fns 2.48} & {\fns 20.1\BergoM} & {\fns 6.68} & {\fns 46.8\BergoM} & {\fns 5.21} \\ 9 & {\fns 106.5\BergoM} & {\fns 2.57} & {\fns 20.5\BergoM} & {\fns 7.37} & {\fns 42.2\BergoM} & {\fns 6.50} \\ 10 & {\fns 130.2\BergoM} & {\fns 2.34} & {\fns 21.0\BergoM} & {\fns 8.00} & {\fns 46.2\BergoM} & {\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} \renewcommand{\baselinestretch}{1.2}% for tabular environment \small \begin{tabular}{ccrcrcrc} \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} \\ \cline{3-8} & & {\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: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: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: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} } \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: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: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: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: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: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: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: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: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: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} } \end{biblist} \end{bibsection}