%% 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}{Multi-level decomposition of Euclidean spheres}
%\title[ short-title ]{ full-title }
\title%
%	[Multi-level decomposition]% short title (maximum of 50 characters)
	{Multi-level decomposition of Euclidean spheres}% full title

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

\begin{Authors}
	\Author{Michael S. Vaz}% full name (without abbreviation)
%% Affiliation: Institution Name [(achronym)], [City], [State], Country
		\Affil{Barco Medical Imaging Systems, Beaverton, OR, USA \\*
			\email{michael.vaz@barco.com}}% institutional e-mail address
	\Author{Atilla P. Kiraly}% full name (without abbreviation)
		\Affil{Siemens Corporate Research, Princeton, NJ, USA\\*
			\email{atilla.kiraly@siemens.com}}% institutional e-mail address
%% If two or more authors share an affiliation, pass an optional
%% argument to the command \Affil when the affiliation is referred
%% for the first time. This argument should be a key consisting 
%% of one or more letters (only letters). Then use the command
%% \Affilref[<KEY>] after any subsequent authors that share the
%% affiliation.
%% Example:
	\Author{Russell M. Mersereau}% full name (without abbreviation)
     \Affil{Georgia Institute of Technology, Atlanta, GA, USA \\*% the affiliation
      	\email{rmm@ece.gatech.edu}}% institutional e-mail address
 \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}
	Mathematical Morphology (MM) is used in medical image analysis for applications such as segmentation and skeletonization.  The available efficient MM methods do not readily adapt to true Euclidean disc/sphere structuring elements (SE), which we are particularly interested in, without sacrificing accuracy for efficiency.  An efficient method for MM using convex/symmetric SE, including Euclidean discs/spheres, is presented.  Performance results for the proposed method are compared to the performance of a commercially available software package for a 512 x 512 x 418 chest CT dataset.  Increasing gains of the method for larger SE are demonstrated, making the method suitable for analysis of high resolution images.  The method is efficient for iso/anisotropic SE.
\end{Abstract}

\begin{Keywords}
	\Index{mathematical morphology},
	\Index{structuring element},
	\Index{sphere},
	\Index{disc},
	\Index{circle}.
\end{Keywords}

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

\section{Introduction}
\label{\fullpaperid:sec:Introduction}% you must prefix your labels {here sec:Introduction) with the string \fullpaperid: (this will be important when combining all the full papers for the final book)

Mathematical morphology (MM) is based on set theory and can be used to analyze image shape features \cites{Vaz:Serra:1982:SE,Vaz:Gonzalez:1992:SE,Vaz:Droogenbroeck:1996:SE}. MM is useful for a wide variety of applications including object recognition, image segmentation, and industrial inspection \cite{Vaz:Shih:1989:SE}.  In medical imaging, MM is used for such applications as brain segmentation from MR images \cite{Vaz:Dogdas:2002:SPIE} and airway/vessel segmentation from CT images \cites{Vaz:Kiraly:2002:SE,Vaz:Eiho:1997:SE}.

Dilation and erosion are the elementary operations of MM. Other operations such as morphological opening and closing may be formed by combining dilation and erosion in sequence \cites{Vaz:Serra:1982:SE,Vaz:Shih:1989:SE}. A structuring element (SE) is the morphological kernel that is translated over the image domain and compared with the overlapping image region. For flat (binary) SE, the comparison operation is a local maximum for dilation and a local minimum for erosion, where the local neighborhood is defined by the SE's shape.  These min/max comparisons for gray-level images reduce to logical AND/OR operations for erosion/dilation of binary images using flat SE.

Morphological analysis methods used in medical imaging generally require SE of different sizes and shapes. Since some of these methods require multiple larger sized SE, which require greater processing time, such methods can take several hours for a single image \cite{Vaz:Kiraly:2002:SE}.  Long processing times impede clinical utility of applications as well as application/algorithm development itself.  

%%%%%%%%PAGE 2%%%%%%%%%
%Reviwer read this incorrectly linear wrt mass, but polynomial wrt to radius 

The processing time for brute-force implementation of MM is proportional to the discrete mass (number of pixels/voxels) of the SE. As such, for 2D and 3D SE, which are used for medical image processing, the processing time can increase polynomially with respect to SE diameter and can become prohibitive. A multitude of methods have been proposed to accelerate MM and the majority of these involve some form of SE decomposition.  The HGW algorithm is regarded as the most efficient for a 1D straight line SE \cites{Vaz:Soille:1996:SE,Vaz:Soille:2001:SE}; logarithmic decomposition (LD) is also very efficient for 1D and as such may be applied to separable 2D and 3D SE \cites{Vaz:Droogenbroeck:1996:SE,Vaz:Bloomberg:BinMorph}.  However, only a subset of useful SE are separable.  Some methods are applicable only to binary images or require preprocessing/encoding of the image \cites{Vaz:Jones:1994:SE,Vaz:Nikopoulos:2000:SE}; others cater to particular constraints such as the limited region of support of specialized hardware \cites{Vaz:Park:1995:SE,Vaz:Hashimoto:2003:SE}.

We are particularly interested in 2D Euclidean disk and 3D Euclidean sphere SE that incur no further approximation beyond quantization itself, with no restriction placed on the isotropy of the quantization. Efficient methods that have already been proposed for convex/symmetric SE \cites{Vaz:Li:1990:SPIE,Vaz:Zhuang:1992:DeSe} are usually discussed and demonstrated in 2D and don't readily lend themselves to a true 3D extension.  Although the sphere falls into convex/symmetric category, as far as we know, these methods may not be used to obtain the decomposition of a sphere.  Moreover, the definition of convexity used in these methods is not consistent and as such, Euclidean disks/spheres are outside the scope of some such methods.  Methods for decomposition of a disc SE usually sacrifice accuracy for efficiency and still do not lend themselves to the 3D sphere \cites{Vaz:Adams:1993:RadDisk,Vaz:Soille:2001:SE, Vaz:Droogenbroeck:1996:SE,Vaz:Zhuang:1992:DeSe}. Some methods refer to gray-level disk SE as spheres \cites{Vaz:Adams:1993:RadDisk,Vaz:Droogenbroeck:1996:SE}; these SE are not the focus of this paper and are not to be confused with binary 3D spheres.  Although the method using local histograms described in \cite{Vaz:Droogenbroeck:1996:SE} may be applied to spheres, it does not readily lend itself to an efficient vectorized implementation, which is important for certain types of computation platforms.

SE decomposition as a union of partitions is described in \cite{Vaz:Anelli:1998:DeArBiSE} where genetic algorithms are used to compute the decomposition.  A deterministic method for decomposing a 2D disc SE into a union of partitions is introduced in \cite{Vaz:Vaz:2006:EfMatMor}.  However, the decomposition is described analytically and may be difficult to understand and implement.  We present a reworking of this method and describe the decomposition using morphological primitives.  As such, it is now easier to understand/implement.  The new morphological decomposition naturally extends itself to all 2D and 3D convex symmetric SE, such as Euclidean spheres.  We present timing results and comparison for a large CT lung volume dataset. 


%%%%%%%%%%%%%%%%START BELOW%%%%%%%KOKO

%%I took out the subsubsection here since it is a single one. Please avoid subsubsection when they are by themselves

In the remainder of this paper we refer to the following notations and definitions and introduce additional elements as necessary.
 
\begin{tabbing}
$\cap : Intersection$  \hspace{5ex} \=  $\ominus : Erosion$   \hspace{5ex} \=  $\circ : Morphological$ $opening$ \kill
$\cup : Union$  \> $\oplus : Dilation$  \> $\bullet : Morphological$ $closing$ \\
$\cap : Intersection$ \>  $\ominus : Erosion$ \>  $\circ : Morphological$ $opening$ 
\end{tabbing}


\vspace{-4mm}
\begin{definition}[Convex SE]
	\label{\fullpaperid:convexSE}
If the discrete SE is equal to the set of all voxels that fall inside its Euclidean convex hull \cite{Vaz:Soille:2003:MatMedImg} it is considered to be convex.  Continuous domain convex shapes that are discretized satisfy this definition of digital convexity.	\qed
\end{definition}

\begin{definition}[Symmetric SE]
	\label{\fullpaperid:symmetricSE}
A symmetric SE has a clearly defined center.  When the center of a symmetric SE is translated to the origin of a standard Cartesian coordinate system, the translated SE would be symmetric about the $x = 0, y = 0,$ and $z = 0$ planes. \qed
\end{definition}

\begin{definition}[Sparse SE]
	\label{\fullpaperid:sparseSE}
In general, a 3D SE that is not face-connected, or a 2D SE that is not edge-connected, is said to be sparse. \qed
\end{definition}

\vspace{-3mm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%METHODS SECTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

		
\section{Method}
\label{\fullpaperid:sec:Method}

\subsection{Overview of the proposed decomposition}
	
	The proposed method decomposes the SE into a union of partitions (Equation \ref{\fullpaperid:eqn:1a}), where each partition, $P_i$, consists of two factors; the first is the largest cube, $C_i$, that can morphologically open the partition without change, and the second is a sparse factor $S_i$ (Equation \ref{\fullpaperid:eqn:1b}).  This method will decompose any 2D or 3D convex/symmetric SE.

\vspace{-3mm}
	  \begin{subequations}
    \begin{eqnarray}
    SE & = & P_{1} \cup P_2 \cup P_{3} \ldots  \label{\fullpaperid:eqn:1a} \\
    SE & = & ( C_1 \oplus S_1 ) \cup (C_2 \oplus S_2 ) \cup ( C_3 \oplus S_3 ) \ldots  \label{\fullpaperid:eqn:1b} \\
    SE & = & ( C_1 \oplus S_1 ) \cup ((C_1 \oplus L_2) \oplus S_2) \cup ((C_1 \oplus L_2 \oplus L_3) \oplus S_3) \ldots \label{\fullpaperid:eqn:1c} \\ 
    SE & = & C_1 \oplus (S_1 \cup L_2 \oplus ( S_2 \cup (L_3 \oplus S_3) \ldots )) \label{\fullpaperid:eqn:1d}
    \end{eqnarray}
    \label{\fullpaperid:eqn:1}
  \end{subequations}

\vspace{-3mm}
As a consequence of convexity of the underlying SE, $C_1$ is a factor of $C_2$ and $C_2$ is a factor of $C_3$ and so on (Equation \ref{\fullpaperid:eqn:1c}).  The $L_i$ factors that relate the different $C_i$ factors are described below and illustrated in Figure \ref{\fullpaperid:fig:3}.  This allows us to further decompose $C_i$ and reuse morphological comparisons across the cubic factors (Equation \ref{\fullpaperid:eqn:1d}).

Figure \ref{\fullpaperid:fig:1} illustrates Equation \ref{\fullpaperid:eqn:1a} and shows the three partitions obtained by using the proposed method to decompose the example 2D SE shown on the right hand side.  The cubic and sparse factors of each of these partitions are shown in Figure \ref{\fullpaperid:fig:2}, which illustrates Equation \ref{\fullpaperid:eqn:1b}.

		
%%%%%%%%%%%%%%%FIGURE 1
	\begin{figure}[ht]
	\centering
%		\includegraphics[0mm,0mm][89.6789mm,23mm]{\fullpaperpdirectory/images/FIG1.png}
%		\includegraphics[bb = 0mm 0mm 89.6789mm 23mm]{\fullpaperpdirectory/images/FIG1.png}
		\includegraphics[width=0.76\hsize]{\fullpaperpdirectory/images/FIG1.png}
		 % \hspace{0.02\hsize}
	\caption{The proposed method decomposes the example SE (right) into a union of partitions $P_{1}$, $P_{2}$, and $P_{3}$.  The origin of each partition is the center.  The partitions overlap; this is possible due to the idempotency property of comparison operations (OR, AND, min, max) used to combine the partitions.}
	\label{\fullpaperid:fig:1}
%	\vspace{-4mm}
\end{figure}
	
%%%%%%%%%%%%%%%FIGURE 2	 MORE DETAILED CAPTION WOULD BE NICE
	\begin{figure}[ht]
	\centering
		%\includegraphics[width=0.9\hsize]{\fullpaperpdirectory/images/FIG2.png}}	
%		\includegraphics[0mm,0mm][68.5689mm,74.4mm]{\fullpaperpdirectory/images/FIG2.png}
		\includegraphics[width=0.57\hsize]{\fullpaperpdirectory/images/FIG2.png}
		 % \hspace{0.02\hsize}
	\caption{Each partition $P_{i}$ has a cubic factor $C_{i}$ and a sparse factor $S_{i}$, where $i = 1, 2, 3$ in this example.  The origin of each factor is at its center.}
	\label{\fullpaperid:fig:2}
%	\vspace{-4mm}
\end{figure}
	

\subsection{Usage of the decomposition}
	
\subsubsection{Binary image dilation and erosion} %%Some text can help clarify this section
	
\vspace{-2mm}
\begin{subequations}
    \begin{eqnarray}
    I \oplus SE & = & (I \oplus P_{1}) \cup (I \oplus P_2) \cup (I \oplus P_{3}) \ldots  \label{\fullpaperid:eqn:2a} \\
    I \oplus SE & = & ( I \oplus C_1 \oplus S_1 ) \cup (I \oplus C_2 \oplus S_2 ) \cup (I \oplus C_3 \oplus S_3 ) \ldots  \label{\fullpaperid:eqn:2b} \\ 
    I \oplus SE & = & I \oplus C_1 \oplus (S_1 \cup L_2 \oplus ( S_2 \cup (L_3 \oplus S_3) \ldots )) \label{\fullpaperid:eqn:2d}
    \end{eqnarray}
    \label{\fullpaperid:eqn:2}
  \end{subequations}

\vspace{-4mm}
Substituting $\oplus$ with $\ominus$ and $\cup$ with $\cap$ yields the equations for erosion.

\subsubsection{Gray-level image dilation and erosion}

\begin{equation}
I \oplus SE  =  \mathrm{Supermum} (  (I \oplus P_1) , (I \oplus P_2) , (I \oplus P_3) \ldots  )
\label{\fullpaperid:eqn:4}
\end{equation}

\vspace{-5mm}
\begin{equation}
I \ominus SE  =  \mathrm{Infimum} (  (I \ominus P_1) , (I \ominus P_2) , (I \ominus P_3) \ldots  )
\label{\fullpaperid:eqn:5}
\end{equation}


\subsubsection{Efficient MM for the cubic factors $C_i$}

The cubic factors of each partition may be decomposed further.  Logarithmic decomposition (LD) \cite{Vaz:Droogenbroeck:1996:SE} is simple and can yield efficiency through computation reuse across partitions. The logarithmic factors for the considered example are illustrated in Figure \ref{\fullpaperid:fig:3}.  The figure shows how the dilation or erosion for all three cubic factors can be obtained with merely 9 comparison operations (ops) per output pixel or voxel.  This can be accomplished by using a cascaded implementation as is suggested in Equation \ref{\fullpaperid:eqn:2d}.

While we favor the cascaded implementation with LD as described above, the HGW algorithm \cite{Vaz:Soille:1996:SE} may be used to efficiently compute the morphology result for the cubic factors.  However this algorithm is more complex and does not gain from computation reuse that is possible due to the fact that $C_i$ is a factor of $C_{i+1}$.   If the user prefers parallelization, using the proposed decomposition in the form shown in Equation \ref{\fullpaperid:eqn:2b} with the HGW algorithm might be a good choice.

			%%%%%%%%%%%%%%%%START ABOVE%%%%%%%
	%%%%%%%%%DONE%%%%%%%%%%%%%
	
	
%%%%%%%%%%%%%%%FIGURE 3
	\begin{figure}[ht]
	\centering
%		\includegraphics[0mm,0mm][94.0285mm,61mm]{\fullpaperpdirectory/images/FIG3.png}
		\includegraphics[width=0.82\hsize]{\fullpaperpdirectory/images/FIG3.png}	
		 % \hspace{0.02\hsize}
	\caption{$C_1$ is a factor of $C_2$ and $C_2$ is a factor of $C_3$; as such, we may reuse computation across cubic factors.  We use logarithmic decomposition (LD) to simply and efficiently perform the dilation/erosion of Ci.  The origin of each $C_i$ should be in its center.  This can be satisfied either by assigning the origin of each logarithmic factor $L_{ix}$ accordingly, or by updating the origin on each $C_i$.  This is possible due to translational invariance property of dilation.}
	\label{\fullpaperid:fig:3}
\end{figure}
	
%	\vspace{-8mm}
\vspace{-6mm}
\subsection{Decomposing a SE using the proposed method}
	
	Figure \ref{\fullpaperid:fig:4} illustrates the steps for decomposing a SE.   The proposed decomposition follows the steps below:

	\begin{enumerate}
	\item Initialization: set CSE (current SE) to be equal to SE.
	\vspace{-2mm}
	\item Find the largest cube with which CSE might be morphologically opened without change.  This is the cubic factor $C_i$.
	\vspace{-2mm}
	\item Find the corresponding sparse factor $S_i$. (See discussion below and Figure \ref{\fullpaperid:fig:5}).  Now that $C_i$ and $S_i$ have been found we have essentially decomposed a partition from CSE.  %% THIS IS NOT A STEP
	\vspace{-2mm}
	\item Update CSE with the subset that remains to be decomposed (RSE: remaining SE) as shown in Figure \ref{\fullpaperid:fig:4}.
	\vspace{-2mm}
	\item Go to Step 2 and repeat until RSE is NULL.
	\vspace{-3mm}
	\end{enumerate}


%%%%%%%%%%%%%%%FIGURE 4
	\begin{figure}[ht]
	\centering
%		\includegraphics[0mm,0mm][94.0285mm,61mm]{\fullpaperpdirectory/images/FIG3.png}
		\includegraphics[width=0.74\hsize]{\fullpaperpdirectory/images/FIG4.png}	
		 % \hspace{0.02\hsize}
	\caption{Decomposing the SE.  At each iteration the current SE (CSE) is updated to reflect the subset of the SE that remains to be decomposed (RSE).  $C_i$ is the largest cube with which $CSE_i$ can be morphologically opened without change.  We can then determine $S_i$ which completes the task of determining $P_i$ and then determine the subset of $CSE_i$ that remains to be decomposed $RSE_i$.  The next iteration begins with updating $CSE_{i+1}$ by setting it equal to $RSE_i$.}
	\label{\fullpaperid:fig:4}
%	\vspace{-5mm}
\vspace{-4mm}
\end{figure}


	
\subsubsection{Determining $S_i$}
Figure \ref{\fullpaperid:fig:5} illustrates the process to determine $S_2$ for the example in Figures \ref{\fullpaperid:fig:1} - \ref{\fullpaperid:fig:4}.
	\begin{enumerate}
	\item Once $C_i$ is available, obtain candidate sparse factors $ST_1$, $ST_2$, $ST_3$.
	\vspace{-2mm}
	\item Test each $ST_j$ for sparseness (see sub-section below).
	\vspace{-2mm}
	\item Select the $ST_j$ that is sparse and has the largest number of voxels (largest discrete mass) and assign it to $S_i$.
%	\vspace{-2mm}
	\end{enumerate}
		% the user is free means that the method is not automatic! You must give a precise method, any variations can be mentioned in the discussion
		
		

\subsubsection{Testing for sparseness}

We enforce that our sparse factor is not face-connected for 3D.  For 2D, it should not be edge-connected.  There are a few exceptions to this rule, such as if there is connectivity with the central voxel (local origin) during the first iteration of decomposition, i.e., during the process of determining $S_1$.

%%%%%%%%%%%%%%%FIGURE 5
	\begin{figure}[ht]
	\centering
%		\includegraphics[0mm,0mm][95.0202mm,67mm]{\fullpaperpdirectory/images/FIG5.png}
		\includegraphics[width=0.74\hsize]{\fullpaperpdirectory/images/FIG5.png}	
		 % \hspace{0.02\hsize}
	\caption{Determining $S_2$ for the considered 2D example.  Obtain test cubic factors $CT_j$, by dilating $C_i$ with $T_j$, where j = 1, 2, and 3 for 2D SE.  As illustrated above, we can then obtain a set of candidate $RSE$,  $RT_j$, and candidate sparse factors $ST_j$.  Determine the subset of $ST_j$  that is indeed sparse and from this subset pick the "best" one and assign it to $S_i$.  We apply a greedy criterion for "best", which is to pick the sparse $ST_j$  that has the largest number of pixels/voxels.  The implementation tests each $ST_j$ for sparseness.  Of course, once a particular $ST_j$ is selected, we can simply assign its corresponding $RT_j$ to $RSE_i$.}
	\label{\fullpaperid:fig:5}
%	\vspace{-5mm}
\vspace{-4mm}
\end{figure}
	
\subsection{Method for 3D SE}
	
\begin{itemize}
\item Once the decomposition is available for a 3D SE, the usage is identical to the 2D case discussed above.
\vspace{-2mm}
\item Decomposition is slightly more elaborate in the stage of selecting a sparse factor.  Whereas for the 2D case we have three $ST_j$ sparse factor candidates to consider, for the 3D case we have seven.  This is a consequence of now having $T_j$, where $j = 1, 2, \ldots 7$.  These $T_j$ are the bar-3 and the square-3, each in the horizontal, vertical and axial orientations as well as the cube-3 (cube with discrete diameter 3).
\end{itemize}

\subsubsection{3D sphere example}	 %% ONLY ONE SUBSUB SECTION HERE, PERHAPS THIS SHOULD BE REMOVED?
Figures \ref{\fullpaperid:fig:6}, \ref{\fullpaperid:fig:7} and \ref{\fullpaperid:fig:8} illustrate the cascaded implementation of the proposed method for a radius-5.5 Euclidean sphere obtained under isotropic unit quantization.  The method decomposes this SE into 3 partitions, as illustrated in the figures.  Only the upper half of the SE is illustrated without loss of generality due to symmetry.  The cubic factors $C_1$, $C_2$, $C_3$ are equal to cube-3, cube-5, cube-9 respectively.  The origin is marked in dark gray (center of slice 6 of 11) in each of the figures and is not part of $S_1$ or $S_2$, which are depicted in medium gray; $S_3$, however, is just the single voxel at the origin.  The medium and light gray voxels together depict a particular partition.  The gray and white voxels together in Figure \ref{\fullpaperid:fig:7} illustrate the union of $P_1$ and $P_2$, where the white pixels indicate the portion of the SE that is exclusively covered by $P_1$.  The union of the gray and white voxels in Figure \ref{\fullpaperid:fig:8} depicts the union of $P_1$, $P_2$ and $P_3$ and is thus the SE itself.  The white voxels in Figure \ref{\fullpaperid:fig:8} indicate the region that is covered exclusively by $P_1 \cup P_2$.

%\setcounter{totalnumber}{4}
%%%%%%%%%%%%%%%FIGURE 6
	\begin{figure}[tp]
	\centering
%		\includegraphics[0mm,0mm][66.5271mm,42mm]{\fullpaperpdirectory/images/FIG6.png}
		\includegraphics[width=0.57\hsize]{\fullpaperpdirectory/images/FIG6.png}
	%	\hspace{0.02\hsize}
%	\caption{Upper half of $P_1$ (of 3) for a radius-5.5 Euclidean sphere at unit quantization. Dark Gray - origin, Light \& Medium Gray - $P_1$, Medium Gray - $S_1$.  $C_1$ is cube-3.}
	\caption{Upper half of $P_1$ for a radius-5.5 Euclidean sphere at unit quantization. The origin is the dark gray voxel in the center of slice 6, $P_1$ is the union of the light and medium gray voxels, medium gray - $S_1$.  $C_1$ is cube-3.}
	%\hspace{-50mm}
%	\vspace{-6mm}
\vspace{-1mm}
	\label{\fullpaperid:fig:6}
\end{figure}
%%%%%%%%%%%%%%%FIGURE 7
	\begin{figure}[tp]
	\centering
%		\includegraphics[0mm,0mm][66.5271mm,42mm]{\fullpaperpdirectory/images/FIG7.png}
		\includegraphics[width=0.57\hsize]{\fullpaperpdirectory/images/FIG7.png}	
	%	\hspace{0.02\hsize}
	\caption{Foreground is $P_1 \cup P_2$ for the radius-5.5 sphere (only upper half of the sphere is shown). Dark gray - origin, light and medium gray together - $P_2$, medium gray - $S_2$, white - region of $P_1$ that does not overlap with $P_2$.  $C_2$ is cube-5.}
		%\hspace{-10mm}
	\label{\fullpaperid:fig:7}
	\vspace{-6mm}
\end{figure}
%%%%%%%%%%%%%%%FIGURE 8
	\begin{figure}[tp]
	\centering
%		\includegraphics[0mm,0mm][66.5271mm,42mm]{\fullpaperpdirectory/images/FIG8.png}
		\includegraphics[width=0.57\hsize]{\fullpaperpdirectory/images/FIG8.png}	
%		\hspace{0.02\hsize}
	\caption{Dark gray - origin as well as $S_3$, light gray and origin together - $P_3$ = $C_3$ = cube-7. White - region of $P_1 \cup P_2$ that does not overlap with $P_3$.  The union of all foreground voxels is the union of all three partitions, and is the radius-5.5 Euclidean isotropic sphere.}
	\label{\fullpaperid:fig:8}
%	\vspace{-4mm}
\end{figure}


\section{Results and analysis}
\label{\fullpaperid:sec:Results}

%%%%%%%%%%%%%%%FIGURE 9 PROBLAMATIC IN SPACING
	\begin{figure}[!ht]
	\centering
%		\includegraphics[0mm,0mm][88.6787mm,43mm]{\fullpaperpdirectory/images/FIG9.png}
		\includegraphics[width=0.9\hsize]{\fullpaperpdirectory/images/FIG9.png}	
		 % \hspace{0.02\hsize}
	\caption{Processing times for gray-level dilation (binary SE) using sphere SE for a $512 \times 512 \times 418$ chest CT image.  Results obtained on and Intel Pentium IV Xeon Dual CPU 2.00 GHz platform with 2.00 GB of RAM.}
	\label{\fullpaperid:fig:9}
	\vspace{-4mm}
\end{figure}


The proposed method was implemented and compared to the commercially available SDC morphology toolbox \cite{Vaz:SDC:toolbox}.  Dilation on a gray-level chest CT image, sized $512 \times 512 \times 418$, was performed using various spherical SE. The standard "brute-force" direct implementation was also compared. Timing results were obtained for spherical SE with discrete radii 1-12. A dual-CPU Intel 2.0 GHz $Xeon$\texttrademark\ computer with 2 GB of RAM was used. The proposed method was implemented as a single-threaded executable. As is evident in Figure \ref{\fullpaperid:fig:9}, the processing time for our own implementation of the direct method becomes prohibitive for larger SE.  The algorithmic advantage of the proposed method facilitates gains over the SDC Morphology Toolbox that increase as SE size increases.


	In \cite{Vaz:Park:1995:SE}, the four factors on the right of Figure \ref{\fullpaperid:fig:10} are offered as an optimal decomposition for the convex/symmetric SE shown to their left ("Iteration 2").  The proposed method will decompose the SE into 2 partitions whose cubic factors are square-3 and square-7 and obtain the MM result in 13 ops per output pixel (4 ops for $C_1$, 4 ops for $L_2$, and 5 ops for the sparse factors and the union), while the 4 factor decomposition requires 16 ops.

%%%%%%%%%%%%%%%FIGURE 10
	\begin{figure}[ht]
	\centering
%		\includegraphics[0mm,0mm][77mm,28mm]{\fullpaperpdirectory/images/FIG10.png}
		\includegraphics[width=0.67\hsize]{\fullpaperpdirectory/images/FIG10.png}	
		 % \hspace{0.02\hsize}
	\caption{[Left] 2-partition decomposition of a 2D SE using the proposed method.  Cubic factors are square-3 and square-7.  Dark gray - origin (center): not part of $S_1$, however $S_2$ is the pixel at the origin.  Light and medium gray together - the partition, medium gray (only in iteration 1) - $S_1$, white - region of $P_1$ that does not overlap with $P_2$.  [Right] 4-factor decomposition for the same SE from Example 1 in \cite{Vaz:Park:1995:SE}.  The proposed method is more efficient.}
	\label{\fullpaperid:fig:10}
		\vspace{-4mm}
\end{figure}

The number of comparison ops required per output voxel for dilation or erosion using discrete Euclidean sphere SE of radii 1-40 was calculated for the direct implementation, the proposed method, and two variations of the proposed method as shown in Figure \ref{\fullpaperid:fig:11}.  The op count for the direct implementation is one less than the number of voxels in the SE. One flavor of the proposed method was to consider the sphere SE as the collection of its 2D slices, implement the MM for each of the unique slices separately using the proposed method and combine the results.  This implementation is less efficient and requires more memory.  It illustrates that a true 3D method is superior to a 2D method directly applied to 3D SE, as would be required in order to make use of most of the applicable methods available in the literature.

The majority of comparison ops of the proposed method occur by way of the sparse factors of the partitions: for a radius-20 sphere SE, the proposed 2-canvas implementation (Equation \ref{\fullpaperid:eqn:2d}) would incur 600 comparisons per output voxel, 36 of which would be due to the cubic factors while the remaining 564 would be due to the sparse factors. By decomposing the sparse factors, this could be reduced to 292, thus reducing the total number to 328 ops per output voxel.  While this could be considered an optimization, it comes at the expense of a third copy of the image volume for interim work, which may be undesirable for computation platforms with memory capacity constraints.

%%%%%%%%%%%%%%%FIGURE 11
	\begin{figure}[htp]
	\centering
%		\includegraphics[0mm,0mm][78.2047mm,43mm]{\fullpaperpdirectory/images/FIG11.png}
		\includegraphics[width=0.9\hsize]{\fullpaperpdirectory/images/FIG11.png}	
		 % \hspace{0.02\hsize}
	\caption{Number of comparison ops required per output voxel for MM using sphere SE.  Note, op count presented on a log scale.  Data obtained analytically.  Proposed implementation follows Equation \ref{\fullpaperid:eqn:2d} and implements the sparse factors directly.  A possible optimization would be to decompose the sparse factors, but this would require another interim copy of the image volume.}
	\label{\fullpaperid:fig:11}
	\vspace{-2.5mm}
\end{figure}

\vspace{0mm}
\section{Discussion and conclusion}
\label{\fullpaperid:sec:Conclusion}
	\vspace{-1.5mm}
An efficient method for MM using Euclidean disk and sphere SE is presented.  These SE do not lend themselves well to factorization methods that can usually be used for convex/symmetric SE.  The proposed method can be applied to all 2D/3D convex/symmetric flat (binary) SE.  The method may be used for binary as well as gray-scale images and does not require any prior analysis or encoding of the image: its gain is not based on image content.  The low memory overhead and efficient looping as well as vectorized implementations promote the utility of the method for development (scripting) as well as application platforms. 

We believe the challenge of decomposing SE into partitions to be an instance of an NP-complete problem.  The proposed method uses locally optimal (greedy) criteria to select a sparse factor during each iteration of the decomposition.  As such, it is not guaranteed to be optimal.  We are investigating optimality relevant to the proposed method and expect to report on how the proposed greedy method compares to an optimal decomposition.  This will of course require a clear definition of optimality.

	The method is robust to scale: it works for iso/anisotropic SE.  As such, it facilitates the development of scale-robust applications/algorithms.  For medical image segmentation, SE may be defined in the context of anatomy instead of using discrete parameters.  The chest CT image used to obtain the results shown in Figure \ref{\fullpaperid:fig:5} was quantized to a $mm$ scale of $0.67 \times 0.67 \times 0.80$.  With the proposed method, an algorithm could define a 3 $mm$ sphere SE and this SE would be discretized according to the quantization scale of the CT image at run-time.

	A comparison of proposed method with others is underway.  An empirical assessment of the gains of the method in the context of real medical applications will be performed.  Since the method lends itself to vectorization, even hardware acceleration by means of a graphics processing unit (GPU) implementation would be possible.




%----------------------------------------

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

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

		\bib{Vaz:Serra:1982:SE}{book}{
			author={Serra, J.},
			title={Image Analysis and Mathematical Morphology},
			publisher={Academic Press},
			address={London},
			date={1982}
		}
		
		\bib{Vaz:Gonzalez:1992:SE}{book}{
			author={Gonzalez, R.},
			author={Woods, R.},
			title={Digital Image Processing},
			publisher={Addison-Wesley Publishing Company},
			address={Reading},
			date={1992},
			edition={1},
			pages={518,519,549}  %%maybe this should go as a note
		}
		
		\bib{Vaz:Droogenbroeck:1996:SE}{article}{
			author={Droogenbroeck, V.},
			author={Talbot, H.},
			date={1996},
			title={Fast computation of morphological operation with  arbitrary structuring elements},
			journal={Pattern Recognition Letters},
			volume={17},
			number={14},
			pages={1451--1460}		
		}
		
		\bib{Vaz:Shih:1989:SE}{article}{
			author={Shih, F. Y.},
			author={Mitchell, O. R.},
			date={1989},
			title={Threshold decomposition of gray-scale morphology into binary morphology},
			journal={IEEE Trans. Pattern Analysis \& Mach. Intel.},
			volume={11},
			number={1},
			pages={31--42}
		}
	
	%%%?SPIE
		\bib{Vaz:Dogdas:2002:SPIE}{article}{
			author={Dogdas, B.},
			author={Shattuck, D. W.},
			author={Leahy, R. M.},
			title={Segmentation of the skull in 3D human MR images using mathematical morphology},
			booktitle={Medical Imaging 2002: Image Processing},
			address={Bellingham},
			publisher={SPIE},
 			book={
  			editor={Sonka, M.},
 				editor={Fitzpatrick, J. M.},
 				publisher={SPIE},
				address={Bellingham},
				volume={4684},
				date={2002},
			},
			conference={
				title={SPIE Medical Imaging},
				date={2002-02-{23--28}},
				address={San Diego, USA},
			},
 	 		pages={1553--1562}
  	}
	
		\bib{Vaz:Kiraly:2002:SE}{article}{
			author={Kiraly, A. P.},
			author={Higgins, W. E.},
			author={McLennan, G.},
			author={Hoffman, E. A.},
			author={Reinhardt, J. M.},
			date={2002},
			title={3D human airway segmentation for clinical virtual bronchoscopy},
			journal={Academic Radiology},
			volume={9},
			number={10},
			pages={1153--1168}
		}
	
	%%? volume %%KOKO 7
		\bib{Vaz:Eiho:1997:SE}{article}{
			author={Eiho, S.},
			author={Qian, Y.},
			date={1997},
			title={Detection of coronary artery tree using morphological operator},
			 book={
			 title={Computers in Cardiology 1997},
  			editor={Murray, A.},
 				editor={Swiryn, S.},
 				volume={24},
 				publisher={IEEE},
				address={Piscataway},
				date={1997},
			},
			conference={
				title={SPIE Medical Imaging},
				date={1997-09-{07--10}},
				address={Lund, Sweden},
			},
			pages={525--528}
		}
	
	
		
	
	
		\bib{Vaz:Soille:1996:SE}{article}{
			author={Soille, P.},
			author={Breen, E. J.},
			author={Jones, R.},
			date={1996},
			title={Recursive implementation of erosions and dilations along discrete lines at arbitrary angles},
			journal={IEEE Transactions on Pattern Analysis and Machine Intelligence},
			volume={18},
			number={5},
			pages={562--567}
		}
		
		
		\bib{Vaz:Soille:2001:SE}{article}{
			author={Soille, P.},
			author={Talbot, H.},
			date={2001},
			title={Directional Morphological Filtering},
			journal={IEEE Trans. Pattern Analysis and Machine Intelligence},
			volume={23},
			number={11},
			pages={1313--1329}
		}
		%ref 10
		\bib{Vaz:Bloomberg:BinMorph}{book}{
			author={Bloomberg, D.},
			title={Implementation Efficiency of Binary Morphology},
			note={Available from: <\url{http://www.leptonica.com}> Access in: 2006-01-05.}
		}
		
		\bib{Vaz:Jones:1994:SE}{article}{
			author={Jones, R.},
			author={Svalbe, I.},
			date={1994},
			title={Morphological filtering as template matching},
			journal={IEEE Trans. Pattern Analysis and Machine Intelligence},
			volume={16},
			number={4},
			pages={438--443}
		}
	
		\bib{Vaz:Nikopoulos:2000:SE}{article}{
			author={Nikopoulos, N.},
			author={Pitas, I.},
			date={2000},
			title={A fast implementation of 3D binary morphological transformations},
			journal={IEEE Transactions on Image Proc.},
			volume={20},
			number={2},
			pages={283--286}
		}
	
		\bib{Vaz:Park:1995:SE}{article}{
			author={Park, H.},
			author={Chin, R. T.},
			date={1995},
			title={Decomposition of arbitrarily shaped morphological structuring elements},
			journal={IEEE Trans. Pattern Analysis and Machine Intelligence},
			volume={17},
			number={1},
			pages={2--15}
		}
		
		
			\bib{Vaz:Hashimoto:2003:SE}{article}{
			author={Hashimoto, R. F.},
			author={Barrera, J.},
			date={2003},
			title={A greedy algorithm for decomposing convex structuring elements},
			journal={Journal of Mathematical Imaging and Vision},
			volume={18},
			pages={269--289},
		}
	
	%%%FIX ?  15 KOKOK
		\bib{Vaz:Li:1990:SPIE}{article}{
			author={Li, D.},
			author={Ritter, G. X.},
			title={Decomposition of separable and symmetric convex templates},
			booktitle={Image Algebra and Morph. Image Proc.},
			address={Bellingham},
			publisher={SPIE},
 			book={
  			editor={Gader, P. D.},
 				publisher={SPIE},
				address={Bellingham},
				volume={1350},
				date={1990},
			},
			conference={
				title={Image Algebra and Morph. Image Proc.},
				date={1990-7-{10--12}},
				address={San Diego, USA},
			},
 	 		pages={408--418}
  	}
	
		%%%FIX ?		
		\bib{Vaz:Zhuang:1992:DeSe}{article}{
			author={Zhuang, X.},
			date={1992},
			title={Morphological structuring function decomposition},
			journal={Proc. IEEE Conf. Computer Vision and Pattern Recognition},
			volume={},
			number={},
			conference={
				title={Computer Vision and Pattern Recognition},
				date={1992-6-{15--18}},
				address={San Diego, USA},
			},
			pages={566--571},
%			month={6}
		}



		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%COMPLETE%%%%%%%%%%%%%%%%%
		\bib{Vaz:Adams:1993:RadDisk}{article}{
			author={Adams, R.},
			date={1993},
			title={Radial decomposition of disks and spheres},
			journal={Computer Vision, Graphics, and Image Processing},
			volume={55},
			number={5},
			pages={325--332},
%			month={9}
		}
		
		\bib{Vaz:Anelli:1998:DeArBiSE}{article}{
			author={Anelli, G.},
			author={Broggi, A.},
			date={1998},
			title={Decomposition of arbitrary shaped binary morphological structuring elements using genetic algorithms},
			journal={IEEE Trans. Pattern Analysis and Machine Intelligence},
			volume={20},
			number={2},
			pages={217--224}
		}
		
		\bib{Vaz:Vaz:2006:EfMatMor}{article}{
			author={Vaz, M. S.},
			author={Kiraly, A. P.},
			title={An efficient method for computing mathematical morphology for medical imaging},
			booktitle={Medical Imaging 2006: Image Processing},
			address={Bellingham},
			publisher={SPIE},
 			book={
  			editor={Reinhardt, J. M.},
 				editor={Pluim, J. P.},
 				publisher={SPIE},
				address={Bellingham},
				volume={6144},
				date={2006},
			},
			conference={
				title={SPIE Medical Imaging},
				date={2006-02-{11--16}},
				address={San Diego, USA},
			},
 	 		pages={1894--1902}
  	}
		
		\bib{Vaz:Soille:2003:MatMedImg}{book}{
			author={Soille, P.},
			date={2003},
			title={Morphological Image Analysis: Principles and Applications},
			edition={2},
			address={New York},
			publisher={Springer Verlag Inc.},
			pages={44--45}
		}
		
		\bib{Vaz:SDC:toolbox}{book}{
			author={SDC Information Systems},
			title={SDC Morphological Toolbox for C++},
			note={Available from: <\url{http://www.morph.com}> Access in: 2006-03-10.}
		}
		
	\end{biblist}
\end{bibsection}