\section{Introduction}

Optical communication protocols using time-bin or pulse-position modulation encoding are a good fit for free-space communication systems. Temporal states having bin widths on the order of a few hundred pico-seconds are robust against errors caused by dispersion from the atmosphere~\cite{kral2005optical}, while achieving high data rates. One significant challenge of establishing a free-space optical link is the distortion of the spatial mode of the beam at the receiver due to atmospheric turbulence. A free-space quantum communication scheme, such as quantum key distribution (QKD) using a time-phase encoding~\cite{brougham2016information,islam2017provably}, requires a time-delay interferometer designed to support many spatial modes to perform measurements in the phase-basis.

Previous work has demonstrated classical and quantum interference with multi-mode time-delay interferometers that leverage various methods to enable a high quality interference~\cite{vallone2016interference,zeitler2016super,jin2018demonstration}. An optical beam with many spatial modes could potentially suffer degraded interference due to each mode having a slightly different wavevector $k$, and therefore a different propagation path length through the interferometer. Robust multi-mode interferometers are designed to minimize the path length difference for each mode to allow a complete interference for an arbitrary spatial mode.

One technique used to reduce the mode-dependent path length uses a 4\,-$f$\,-\,imaging system in each arm to relay the mode during propagation and therefore remove the $k$-vector dependence. This approach has demonstrated fringe visibility as high as $93\%$~\cite{zeitler2016super}. Another promising design for enabling high-quality multi-mode interference is a field-widened Michelson-type interferometer without relay optics~\cite{shepherd1985wamdii,gault1985optimization}. One such design has reported $91\,\%$ interference visibility and $85\,\%$ entanglement visibility~\cite{jin2018demonstration}. In this article we discuss the design, construction, and characterization of a multi-mode time-delay interferometer operating at $1550\,$nm. The design presented here uses a Michelson-type layout without relay optics to achieve robust multi-mode interference. Our design features different glass materials in each beam path to achieve improved spatial mode and thermal performance in a compact layout. We demonstrate the application of the interferometer to a single-photon phase-state measurement used in a time-phase QKD scheme.

\section{Multi-Mode Delay Interferometer Design}

A simplified interferometer layout is shown in Fig.~\ref{fig:Int}~\cite{gault1985optimization}. A ray incident from the left at an angle $\theta$ is split by a non-polarizing beamsplitter and each beam is directed to one of two pathways. One path has a short propagation length relative to the other. The short arm is made of material with index $n_{1}$ and thickness $d_{1}$, and the long arm is made of material with index $n_{2}$ and thickness $d_{2}$. Optionally, an air gap of thickness $d_{3}$ can be introduced between the end of the material in the long arm and the reflector. The two beams are reflected at the end of each arm and recombine at the beamsplitter.

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{Int}
\caption[Simplified schematic of a Michelson-type interferometer.]{Simplified schematic of a Michelson-type interferometer. This setup includes a small air gap at the end of the long arm ($d_{2}$).}
\label{fig:Int}
\end{center}
\end{figure}

The optical path length difference $\Delta$ between the two arms of the interferometer as a function of incident angle $\theta$ is given by~\cite[]{shepherd1985wamdii}

\begin{equation}\label{eq:OPD_1}
\begin{split}
\Delta = 2[n_{3}d_{3}(1-\sin^{2}\theta/n^{2}_{3})^{1/2}\\
+ n_{2}d_{2}(1-\sin^{2}\theta/n^{2}_{2})^{1/2}\\
- n_{1}d_{1}(1-\sin^{2}\theta/n^{2}_{1})^{1/2}],
\end{split}
\end{equation}

\noindent
where $n_{3}=1$ but is kept in the equation for completeness. Typically, a Michelson interferometer is designed for the chief ray to enter at normal incidence and small angular deviations from the chief ray are analyzed. In the design presented here, the input ray must enter at a finite angle $\theta_{0}$ to spatially separate each output beam from the input beam. Access to both outputs are required for the time-phase QKD protocol, discussed later. Unlike the analysis performed in~\cite{shepherd1985wamdii}, where Eq.~\ref{eq:OPD_1} is expanded about a normal incidence, we perform the expansion of Eq.~\ref{eq:OPD_1} about a finite input angle $\theta_{0}$. The first three terms of this expansion are given by

\begin{equation}\label{eq:OPD_2}
\begin{split}
\Delta = 2\left(d_{3}\beta_{3} +  d_{2}\beta_{2} -  d_{1}\beta_{1}\right)\\
 - 2\Phi(\sin\theta - \Phi)\left(\frac{d_{3}}{\beta_{3}} +  \frac{d_{2}}{\beta_{2}} -  \frac{d_{1}}{\beta_{1}}\right)\\
- (\sin\theta - \Phi)^{2}\left(\frac{n^{2}_{3}d_{3}}{\beta^{3}_{3}} +  \frac{n^{2}_{2}d_{2}}{\beta^{3}_{2}} -  \frac{n^{2}_{1}d_{1}}{\beta^{3}_{1}}\right),
\end{split}
\end{equation}

\noindent
where $\Phi = \sin\theta_{0}$ is the input angle and $\beta_{i}^{2} = n_{i}^{2} - \Phi^{2}$ is the effective index of refraction. The goal of our analysis is to satisfy three conditions that are derived from Eq.~\ref{eq:OPD_2}. A full mathematical derivation can be found in Ref.~\cite{cahall2019enabling}.

\section{Design Requirements}

The interferometer is designed to achieve the following; (1) an optical path-length difference (OPD) of 200\,ps (5\,GHz free spectral range), (2) satisfy the field-widened condition (discussed more below) for a large angular bandwidth, and (3) be thermally compensated. A successful design will feature high visibility near unity, long-term temperature stability over a modest temperature range ($\pm0.2^{\circ}$\,C), and have phase-tuning ability over a complete interference fringe.

The OPD of the chief ray is determined by setting $\sin\theta = \Phi$, leaving only the first term of Eq.~\ref{eq:OPD_2}. The path length difference $\Delta$ is in units of length, and therefore must be converted to a time delay $\Delta_{t}$ by dividing it by the speed of light in vacuum $c$ and given by
\begin{equation}\label{eq:TimeDelay}
\Delta_{t} = 2(d_{3}\beta_{3} +  d_{2}\beta_{2} -  d_{1}\beta_{1})/c.
\end{equation}

\noindent
A wide angular bandwidth is necessary for a multi-mode interferometer design because a multi-mode beam has a range of wavevectors and each will take a slightly different path through the interferometer. To achieve high visibility, the change in OPD as a function of input angle must be minimized. Field-widening is achieved by proper material choices and device geometry that make the second term in Eq.~\ref{eq:OPD_2} vanish. The field-widening condition is
\begin{equation}\label{eq:FieldWiden}
 \frac{d_{3}}{\beta_{3}} + \frac{d_{2}}{\beta_{2}} - \frac{d_{1}}{\beta_{1}} = 0.
 \end{equation}

\noindent
The result of satisfying Eq.~\ref{eq:FieldWiden} is that the OPD in Eq.~\ref{eq:OPD_2} becomes first-order insensitive to deviations from the input angle of the chief ray~\cite{gault1985optimization,shepherd1985wamdii}.

Lastly, thermal compensation provides long-term stability in an imperfectly controlled environment. Temperature fluctuations cause OPD changes due to expansion and contraction of the physical lengths of the glass material, as well as temperature-dependent changes in the index of refraction~\cite{shepherd1985wamdii}. The thermal drift is characterized by taking the partial derivative of Eq.~\ref{eq:OPD_2} with respect to temperature $T$, as

\begin{equation}\label{eq:ThermalDrift}
\begin{split}
\frac{\partial\Delta}{\partial T} = 2d_{2}\left[\beta_{2}\alpha_{2} + \frac{n_{2}}{\beta_{2}}\left(\frac{\partial n_{2}}{\partial T}\right)\right]\\
- 2d_{1}\left[\beta_{1}\alpha_{1} + \frac{n_{1}}{\beta_{1}}\left(\frac{\partial n_{1}}{\partial T}\right)\right],
\end{split}
\end{equation}

\noindent
where $\alpha_{i}=(1/d_{i})(\partial d_{i}/\partial T)$ is the coefficient of thermal expansion (CTE). This analysis assumes that the index and length change of the air gap with temperature is negligible and therefore those terms vanish when taking the derivative. The interferometer is made a-thermal by setting $\partial\Delta/\partial T = 0$, imposing a constraint on the path lengths $d_{1}$ and $d_{2}$, given by

\begin{equation}\label{eq:Athermal}
\left[\beta_{1}\alpha_{1} + \frac{n_{1}}{\beta_{1}}\left(\frac{dn_{1}}{dT}\right)\right]d_{1} = \left[\beta_{2}\alpha_{2} + \frac{n_{2}}{\beta_{2}}\left(\frac{dn_{2}}{dT}\right)\right]d_{2}.
\end{equation}

\section{Performance Simulation}

We conducted a survey of optical glass with the goal of identifying two glass types to construct the interferometer (See Appendix~\ref{appendix:material}) . The basic guidelines we used for making the glass choices are: (1) availability, (2) a large difference in index of refraction between the two materials, and (3) small relative changes in refractive index and length with temperature. A large difference in the index of refraction will enable smaller physical lengths to satisfy a given time-delay.

The finalized design uses N-SF66 and N-BK10 (Schott AG) for the long arm and short arm respectively. At the design wavelength of 1550\,nm, N-BK10 has an index of refraction of $n_{1}=1.4823$, $\partial n_{1}/\partial T = 1.4154\times10^{-6}\,\text{K}^{-1}$, and a CTE of $\alpha_{1}=5.8\times10^{-6}\,\text{K}^{-1}$. N-SF66 has an index of refraction of $n_{2}=1.8660$, $\partial n_{1}/\partial T = -1.5575\times10^{-6}\,\text{K}^{-1}$, and a CTE of $\alpha_{1}=5.9\times10^{-6}\,\text{K}^{-1}$. The index and change in index for each material are calculated at $\lambda=1550\,$nm using the Sellmeier equation and its derivative with respect to temperature. The glass rods are mounted to a substrate of Kovar nickel alloy. The CTE of Kovar is $\alpha_{\text{Kovar}}=5.86\times10^{-6}\,\text{K}^{-1}$ and is closely matched to the glasses used in our design. A good match between the glass and substrate, rather than using an ultra-low expansion (ULE) substrate with a much smaller CTE, is advantageous because matching CTE values between the glass and substrate will minimize the stress on the glass during temperature changes, which can cause birefringence and refraction.

The lengths $d_{1}$ and $d_{2}$ are optimized for an input angle $\theta_{0}=3^{\circ}$ and satisfy the OPD, field-widening, and a-thermal conditions. However, satisfying these three conditions simultaneously leads to large physical dimensions of the interferometer ($>50$\,mm). The negative consequence of large glass pathways are: (1) large glass path lengths increase the thermal mass to be stabilized and can increase the settling time for phase tuning: (2) as the path length increases so does the clear aperture required for a beam input at a fixed input angle; and (3) longer path-lengths will be more sensitive to overlap misalignment at the beamsplitter. Misalignment is especially sensitive for multi-mode interference because the spatial frequency of the mode can easily be on the order of the alignment tolerance, whereas in single-mode interference, the mode size is large compared to the alignment tolerance.

The goal of the design presented here is to minimize the physical dimensions of the glass beam paths. Minimizing the glass length will allow for a more robust and manufacturable setup. To meet this goal, we relax the thermal requirements of the design so that the a-thermal condition is satisfied only approximately, while the OPD and field-widening conditions are satisfied exactly. We find a unique solution that satisfies Eq.~\ref{eq:TimeDelay} and Eq.~\ref{eq:FieldWiden} by eliminating the air gap $d_{3}$. The OPD of the chief ray and field-widening condition for an interferometer without an air gap $d_{3}$ are given by

\begin{equation}\label{eq:TimeDelay2}
\Delta_{t} = 2(d_{2}\beta_{2} -  d_{1}\beta_{1})/c,
\end{equation}
and
\begin{equation}\label{eq:FieldWiden2}
\frac{d_{2}}{\beta_{2}} - \frac{d_{1}}{\beta_{1}} = 0,
\end{equation}
respectively. The a-thermal condition remains unchanged when removing the air gap.

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{OPDangle}
\caption[Angular bandwidth about the chief ray of the multi-mode interferometer.]{Simulated change in the optical path difference as a function of deviations of the input about the chief ray of the multi-mode interferometer.}
\label{fig:OPDangle}
\end{center}
\end{figure}

The absence of an air gap leaves two equations and two free parameters, $d_{1}$ and $d_{2}$. The path lengths that satisfy Eq.~\ref{eq:TimeDelay2} and Eq.~\ref{eq:FieldWiden2} are $d_{1} = 34.5\,$mm and $d_{2} = 43.5\,$mm. At an input angle of $\theta_{0}=3^{\circ}$ and a beam diameter $\sim2\,$mm (discussed in more detail below), the required clear aperture is $\sim10\,$mm. A clear aperture of this size allows the square cross-section of the glass rods to be $12.7\,\text{mm}\times12.7\,\text{mm}$, and leads to manageable total thermal mass. Additionally, a 10\,mm clear aperture is $<80\,\%$ of a 12.7\,mm surface, which allows each surface to be polished easily with standard glass processing capabilities.

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{OPDdrift}
\caption[Simulated thermal performance of the multi-mode interferometer.]{Simulated thermal performance of the multi-mode interferometer. (a) Thermal stability for small temperature drifts of the complete system. (b) Interference phase tuning via control over the temperature of the short arm $d_{1}$.}
\label{fig:OPDdrift}
\end{center}
\end{figure}

The angular bandwidth of the interferometer is characterized by plotting the OPD in Eq.~\ref{eq:OPD_2} for a range of input angles $\theta$ around the chief ray $\Phi = \sin3^{\circ}$, shown in Fig.~\ref{fig:OPDangle}. Recall that satisfying the field-widening condition given in Eq.~\ref{eq:FieldWiden2} implies that the first term in the angular dependence of the OPD (Eq.~\ref{eq:OPD_2}) is proportional to $\sin^{2}\theta$.

The thermal stability of the setup is estimated by calculating Eq.~\ref{eq:ThermalDrift} for small temperature changes $\Delta T$, shown in Fig.~\ref{fig:OPDdrift}(a). The simulation shows that thermal drift in the OPD of $\leq\lambda/50$ is achievable for modest temperature fluctuations of $\pm0.1^{\,\circ}$C. The phase of the interference is adjusted by controlling the temperature in one of the glass arms, the short arm for example, to $\pm1^{\,\circ}$C. Temperature control over this range will result in a change in the OPD from $-\lambda/2$ to $\lambda/2$, allowing for a tuning range of $0$ to $2\pi$ for the phase $\phi$, as shown in Fig.~\ref{fig:OPDdrift}(b). Thermal tuning is preferred over mechanical tuning because mechanical tuning techniques such as piezo-driven mirrors can be expensive, require an extra control and feedback system, and have large thermal drift.

\section{Interference Characterization}

The performance of the interferometer is verified using optical beams with single-mode (SM) and multi-mode (MM) intensity profiles. Atmospheric turbulence is emulated by mode mixing in a perturbed multi-mode fiber~\cite[]{ursin2007entanglement,jin2018demonstration}. The experimental setup, shown schematically in Fig.~\ref{fig:System}, consists of a continuous-wave laser diode (Fitel F0L15DCWC) operating at 1550\,nm with a laser controller (Arroyo 6310 Combo Source) that controls the diode current and the temperature of the package. The wavelength of the output for diodes of this type is a function of the diode current and temperature. The controller has an analog voltage input for manual adjustment of the current supplied to the laser diode. A controllable voltage source is used as the scan voltage input for diode current sweeps. Adhesive heater tape (Omega KHLVA-101/5-P) is used on the short arm to provide phase tuning of the interferometer. The light from the laser is directed to the interferometer via a single or multi-mode optical fiber. The output of the interferometer is coupled in to an optical fiber where the output power is measured with an optical power meter (Newport Optical Power Meter 1803C with Model 818-IG sensor head), or the intensity profile of the beam can be viewed with a beam profiling camera (Ophir-Spricon SP907-1550).

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{System}
\caption{Schematic of the experimental setup used to characterize the interference quality and thermal performance of the interferometer (gray box). SM: single-mode, MM: multi-mode, TEC: thermoelectric cooler, ADC: analog-to-digital converter.}
\label{fig:System}
\end{center}
\end{figure}

In a practical setting, the temperature of the interferometer is meant to be passively controlled by modest isolation from the surrounding lab environment. For the purposes of thermal characterization, the experimental setup includes precise temperature control with thermo-electric cooling elements (TEC, Laird 64975-502) and a TEC controller (Thorlabs TED200C) with proportional-integral-derivative (PID) feedback control. Three TEC elements are placed between the bottom of the interferometer base plate and a large thermal bath, which is the optical table in this case. The temperature of the base plate is sensed via high-accuracy thermistors (Vishay NTCALUG02A103F161) placed at various locations on the plate. A thermistor near the center of the plate is used for the TEC feedback control. The analog output from the TEC controller is digitized with an analog-to-digital converter (ADC, National Instruments NI-9239) and recorded. A representative example of the settling time and temperature stability of the interferometer is shown in Fig.~\ref{fig:TempSet}.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=1.0\linewidth]{TempSet}
\caption{Temperature control of the interferometer bottom plate. (a) Settling time of the plate temperature after a step change of $2^{\,\circ}$C. The target temperature is stable after $\sim3$ minutes. (b) Long-term temperature stability at a fixed set point. The temperature is controlled with an accuracy of $\pm0.0015^{\,\circ}$C.}
\label{fig:TempSet}
\end{center}
\end{figure}

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{PhaseTune}
\caption{Interference phase tuning with short arm temperature control. The output from the interferometer is plotted as a function of the square of the voltage supplied to the heater element on the short arm showing a tuning range over a full interference fringe.}
\label{fig:PhaseTune}
\end{center}
\end{figure}

\subsection{Fringe Visibility}

The phase of the interference fringe is scanned by heating one of glass arms and the result is shown in Fig.~\ref{fig:PhaseTune}. To maintain repeatability and execute measurements on faster timescales than thermal settling allows, interference fringes are recorded by scanning the input wavelength rather than scanning the phase of the interferometer. The wavelength of the laser diode is adjusted with an analog voltage input to the laser controller, which controls the current supplied to the diode. As the laser current is scanned, the wavelength changes but so does the total output power.  The output power as a function of the laser current is first calibrated by measuring the power directly from the laser. The change in optical power as a function of diode current is compensated for by dividing the interference fringe measurements by the calibration measurement. The power output of the interferometer is recorded as a function of the input scan voltage value. The quality of the interference is quantified by the visibility $\mathcal{V}$, given by

\begin{equation}\label{eq:Vis}
\mathcal{V}=\frac{\mathcal{P}_{+}-\mathcal{P}_{-}}{\mathcal{P}_{+}+\mathcal{P}_{-}}
\end{equation}

\noindent
where $\mathcal{P}_{+}$ is the maximum and $\mathcal{P}_{-}$ is the minimum value of the interference fringe.

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{Fringes}
\caption[Interference fringes for SM and MM ($50\,\mu$m core fiber) input beams.]{Interference fringes for SM and MM ($50\,\mu$m core fiber) input beams as a function of the input voltage to scan the laser diode current. Each data point (black) is the average of 10 measurements and error bars indicate the standard deviation of those 10 measurements. Red lines show a fit to the data using the function given by Eq.~\ref{eq:Pout}. (a) SM interference fringe. The visibility of the interference is $\mathcal{V}_{\mathrm{SM}}=99.02\pm0.05\,\%$. The high order fringe structure in the constructive interference could be due to reflections from multiple surfaces within the interferometer. (b) Zoomed view of the SM fringe minimum. (c) MM interference fringe. The visibility of the interference is $\mathcal{V}_{\mathrm{MM}}=98.38\pm0.01\,\%$. (d) Zoomed view of the MM fringe minimum.}
\label{fig:Fringes}
\end{center}
\end{figure}

The power $P_{\mathrm{out},\pm}$ measured at the $+/-$ port of the interferometer as a function of small changes of the wavelength $\delta\lambda\ll\Delta$, when the phase of the interference is set to $\phi$, is given by~\cite[]{islam2017robust}

\begin{equation}\label{eq:Pout}
\frac{P_{\mathrm{out},\pm}}{P_{0}}=\frac{a}{2}\left[1\mp\sin\left(k \delta\lambda + \phi\right)\right] + b,
\end{equation}

\noindent
where $P_{0}$ is the input optical power, $a\in\left\{0,1\right\}$ is the insertion loss of the interferometer, $b$ accounts for imperfect destructive interference in the minimum of the fringe, and $k=2\pi/\lambda$ is the wavenumber of the light. The output power is recorded as a function of the laser diode scan voltage, which serves as the wavelength adjustment parameter. The fringe measurement is fit with a curve given by Eq.~\ref{eq:Pout} with $a$, $b$, and $\phi$ as fitting parameters. The visibility of the normalized fringe is extracted from this fit where

\begin{equation}\label{eq:Vis2}
\mathcal{V}=\frac{a}{a + 2b}
\end{equation}

The measurement results from a single-mode input fiber (Corning SMF-28-Ultra) and multi-mode input fiber ($50\,\mu$m core, 0.1\,NA, Coastal Connections M-FUFU-50sL) are shown in Fig.~\ref{fig:Fringes}. The interferometer demonstrates single-mode interference visibility $\mathcal{V}_{\mathrm{SM}}=99.02\pm0.05\,\%$ and multi-mode interference visibility $\mathcal{V}_{\mathrm{MM}}=98.38\pm0.01\,\%$ extracted from the fitted functions. The single-mode results are close to those reported by high-performing commercial interferometers such as those used in a recent QKD demonstration~\cite{islam2017robust,islam2017provably}. Additionally, the multi-mode results presented here are significantly better than those recently reported for multi-mode interferometers~\cite{zeitler2016super,jin2018demonstration}.

\subsection{Spatial Mode Performance}

A robust multi-mode interferometer is required to maintain high visibility independent of the input spatial mode structure. The multi-mode results presented above are optimized and performed on a single beam profile. To investigate the robustness of the performance, we characterize the visibility of the interference fringe as the spatial profile is changed.  The visibility is maximized for the first spatial profile by tuning the beam alignment such that there is good overlap at the beamsplitter. Then, interference fringes are recorded for subsequent spatial profiles after the mode has been changed. The difference between spatial mode profiles is quantified by the two dimensional Pearson correlation coefficient $R$ defined as

\begin{equation}\label{eq:2dcorr1}
R = \frac{\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\left[I_{ij}\left(t_{1}\right)-\mu\left(t_{1}\right)\right]\left[I_{ij}\left(t_{2}\right)-\mu\left(t_{2}\right)\right]}{\sigma(t_{1})\sigma(t_{2})},\\[1.0em]
\end{equation}

\noindent
where $I_{ij}(t_{n})$ is the optical intensity collected at time $t_{n}$ (a frame), indexed by a discrete array of positions $(i,j)$ that correspond to the pixels of the camera, $\mu\left(t_{n}\right)$ is the mean intensity, and $\sigma(t_{n})$ is the standard deviation of the intensity distribution. The correlation coefficient $R$ has a value between $+1$ and $-1$, where $+1 (-1)$ indicates perfect positive (negative) correlation, and values close to $0$ indicate no correlation.

\begin{figure}[th]
\begin{center}
\includegraphics[width=1.0\linewidth]{MMode}
\caption{Spatial mode characterization from a $50\,\mu$m and $105\,\mu$m core fiber. (a) Schematic showing the setup used to adjust the mode profile of the beam. The spatial mode is changed by inducing bends in the fiber controlled by a positioning stage and a micrometer. (b) Spatial mode examples for different micrometer settings from a $50\,\mu$m core fiber. (c) Spatial mode examples for different micrometer settings from a $105\,\mu$m core fiber.}
\label{fig:MMode}
\end{center}
\end{figure}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=1.0\linewidth]{VisCorr}
\caption{Multi-mode interference as a function of the input spatial mode. Correlation coefficient between the first frame and subsequent frames as a function of the micrometer setting for $50\,\mu$m (circles) and $105\,\mu$m (triangles) core fibers is shown against the left vertical axis. Also plotted is the visibility of the interference fringe over the same micrometer range for each fiber, shown in the right vertical axis. Error in the fringe visibility is $\leq5\times10^{-4}$.}
\label{fig:VisCorr}
\end{center}
\end{figure}

Optical intensity profiles are collected by the beam profiling camera and the correlation coefficient is processed in MATLAB. The spatial mode of the beam is changed by introducing a controlled bend in the optical fiber with a micrometer and positioning stage. A schematic depicting the setup to induce changes in the mode is shown in Fig.~\ref{fig:MMode}(a). The fiber is anchored at three different points with the middle point on the positioning stage. As the position of the stage is adjusted, the fiber is bent at two different points that leads to the change in the mode profile. Special care is taken to ensure that all other points in the fiber are kept stationary throughout the experiment. When the fiber is kept secure apart from the micrometer adjustment, the spatial modes can be reproduced as a function of the micrometer position.

The profile of the beam is measured at the output of the fiber (without collimation) using the beam profiler. Examples of the intensity patterns from $50\,\mu$m and $105\,\mu$m core (Thorlabs FG105LVA) fibers for different micrometer positions are shown in Fig.~\ref{fig:MMode}(b) and (c). Note that the experiment is carried out by first recording the mode profiles with the beam profiler at each micrometer position, then performing the interference measurements at each micrometer position. Between the beam profiler measurements and the interference measurements, the output of the fiber must be moved from the beam profiler to the interferometer input. The movement of the fiber will unavoidably change the actual spatial mode that travels through the interferometer from those in the recorded images in Fig.~\ref{fig:MMode}(b) and (c). However, the intensity profile and the \emph{correlation} between profiles at different micrometer positions recorded by the beam profiler are representative of the actual mode passing through the interferometer. Therefore, the mode analysis using the beam profiler adequately represents the spatial mode at the input of the interferometer.

The correlation coefficient and interference visibility as a function of micrometer position are shown in Fig.~\ref{fig:VisCorr}. The correlation between the first measurement and subsequent measurements decays to a value $<0.5$ as the spatial mode becomes more distinguishable. As the spatial mode is changed, the interferometer maintains excellent fringe visibility for each measurement. The quality of the interference using a $50\,\mu$m core fiber has an average value $\mathcal{V}_{\mathrm{avg}}=98.05\pm0.01\%$ and a maximum value $\mathcal{V}_{\mathrm{max}}=98.38\pm0.01\%$. Additionally, the quality of the interference using a $105\,\mu$m core fiber has an average value $\mathcal{V}_{\mathrm{avg}}=97.55\pm0.04\%$ and a maximum value $\mathcal{V}_{\mathrm{max}}=98.14\pm0.02\%$. The quality of these results indicate an interferometer design that is robust for arbitrary input spatial modes.

\subsection{Thermal Performance}

The thermal performance of the interferometer is characterized by measuring the visibility of the interference fringe and the path length shift as a function of temperature. The temperature of the bottom plate is set and maintained by the TEC controller that supplies current to the TEC elements in contact with the substrate of the interferometer. When the setpoint is changed, the system is allowed to come to thermal equilibrium over a period of $\sim30$ minutes. As thermal equilibrium is reached, the path length will drift as a function of time, causing a change in the measured output power. Therefore, when a new setpoint is investigated, the output of the interferometer is monitored as equilibrium is reached. We infer that the interferometer has come to thermal equilibrium once the output power is stabilized, shown in Fig.~\ref{fig:EqTime}. During setpoint changes, the phase of the interference is set such that it lies near the point of the steepest slope to have the maximum sensitivity in determining the equilibrium.

\begin{figure}[thb]
\begin{center}
\includegraphics[width=1.0\linewidth]{EqTime}\\[1.1em]
\caption{Example measurement of the output power from the interferometer while the system reaches thermal equilibrium. This particular set of data shows a temperature change from $21.0^{\,\circ}$C to $21.5^{\,\circ}$C}
\label{fig:EqTime}
\end{center}
\end{figure}

\begin{figure}[thb]
\begin{center}
\includegraphics[width=1.0\linewidth]{TempScan}
\caption{Multi-mode interference as a function of temperature. (a) Fitted functions to the measured interference fringe at each temperature setting. The temperature settings for the curves moving from left to right are $21.0, 21.5, 21.9, 22.0, 22.5,$ and $23.0\,^{\circ}$C. The phase of the interference at $22.0\,^{\circ}$C is set at data point marked with an `X.' The path length shift is calculated using points where the fringe intersects the dotted line. (b) Experimentally measured (black) and theoretically predicted (red) path length change as a function of temperature. (c) Interference visibility as a function of temperature.}
\label{fig:TempScan}
\end{center}
\end{figure}

Once the system has come to thermal equilibrium, a full interference fringe is scanned and fit with the function given by Eq.~\ref{eq:Pout} in the same manner as discussed above. A fringe is recorded at each temperature setting and plotted as a function of the scanning voltage applied to the diode controller, corresponding to the change in the laser wavelength. The results of the curve fitting are shown in Fig.~\ref{fig:TempScan}. Path length drift brought on by the temperature change causes the fringe to move horizontally. While operating at a fixed value of $\phi$ (shown by the `X' on the yellow curve in the figure) a horizontal shift of the fringe manifests in a change in the power measured at the output port of the interferometer. The change in output power describes change in path length.

The path length change and visibility are measured for a temperature range of $\pm1.0\,^{\circ}$C centered at $22.0\,^{\circ}$C. The results of the measurements are shown in Fig.~\ref{fig:TempScan}(b) and (c). A high fringe visibility of $\sim98\,\%$ is maintained across the entire temperature range. The experimentally determined path-length shift as a function of temperature is 130\,nm/$\,^{\circ}$C. This value matches well with the theoretically predicted performance. Therefore, stability on the order of $\lambda/50$ can be achieved with moderate temperature control in a lab environment on the order of $0.25\,^{\circ}$C.

\section{Application to Time-Phase QKD}

Experimental realizations of QKD systems are often hindered by non-idealalities in the equipment used to generate and measure the photons. Rate limitations are dominated by detector saturation and dead-time, and so it is advantageous to increase the channel capacity by encoding bits in a high dimensional Hilbert space. We use a time-phase encoding~\cite{brougham2016information} because it can support a high dimensional encoding. An encoding scheme having $d$-dimensions is often referred to as a qudit protocol, where $d=2$ is the well-known qubit protocol. Qudits can transmit up to $\log_{2}d$ bits per photon, whereas qubit schemes, such as BB84~\cite{bennett1984quantum}, are limited to a single bit per photon. Additionally, high-dimensional schemes can tolerate a higher bit-error-rate, which enable longer communication distances~\cite{scarani2009security}. Recently, a high-dimensional QKD scheme using a $d=4$ time-phase encoding has achieved a record high maximum secret-key-rate (SKR) of 26.2\,Mbit/s, and an SKR of 1.07\,Mbit/s at an equivalent fiber distance of 83\,km~\cite{islam2017provably}.

A time-phase QKD protocol consist of two mutually unbiased bases that are referred to as the temporal (or time) basis and the phase (or frequency) basis. Data is encoded in the temporal basis and transmitted from Alice to Bob, and the phase basis is used to check for the presence of an eavesdropper, Eve. In general, $d-1$ cascaded interferometers are needed to perform a complete measurement of the phase basis~\cite{brougham2013security}. A full description of such a protocol, including a demonstration utilizing a high-dimensional encoding, can be found in Ref.~\cite{islam2017provably}. Here, we use a $d=2$ encoding to demonstrate the performance of our multi-mode interferometer.

The states for a $d=2$ time-phase protocol are shown in Fig.~\ref{fig:2dTimePhase}. Temporal states are constructed by intensity modulation of a continuous-wave (CW) laser source to create short duration wavepackets. The states are highly attenuated to produce weak coherent states with a mean photon number per pulse of $\mu\approx1$. Phase states are constructed using intensity modulators to create the superposition of time states and an additional phase modulator to adjust the phase of the wavepackets in each time bin. The phase states are attenuated so that the probability amplitudes of the entire frame sum to unity. It is important to use a modulation scheme with a high extinction ratio to prevent errors caused by background photons leaking through the modulators.

A measurement in the temporal basis is performed with a single-photon detector, where the electrical signal generated from a detection event is recorded with a high resolution time-tagger. The single-photon detector and processing electronics, including amplifiers and time-tagging electronics, must have low timing uncertainty (detection jitter) to localize the time-of-arrival of the photon with high accuracy~\cite{marsili2013detecting,cahall2018scalable}. The arrival time of the single-photon detection event is compared to the shared global clock to declare the result of the measurement.

Phase-state measurements require a time-delay interferometer, shown schematically in Fig.~\ref{fig:2dTimePhase}(c). A phase state entering the interferometer is split by a 50/50 beamsplitter with half of the wavepacket traveling through a long path and half through a short path. The difference in path lengths is chosen such that the time-delay is equal to $\tau$, or one time bin width. Therefore, the wavepackets recombining at the second beamsplitter are shifted by one time bin relative to one another. An additional phase adjustment $\phi$ is available to tune the phase of the interference fringe, where $\phi=0$ in the $d=2$ scheme.

A phase state that passes through the interferometer will result in constructive or destructive interference in the shaded time-bin shown in Fig.~\ref{fig:2dTimePhase}(e). There is a one-to-one mapping between the phase-state input and the detector in which constructive interference is observed. Therefore, constructive interference is detected at a unique detector, with destructive interference at all other detectors. In the $d=2$ protocol, a photon detection event in the shaded time bin of detector $D0$ ($D1$) indicates a measurement outcome of $\ket{f_{0}}$ ($\ket{f_{1}}$). The lobes in the interference pattern lying outside the shaded time bin are considered losses and do not contribute to the outcome of the measurement. The phase-basis measurement succeeds at a rate of $1/d$, or $50\%$ for the $d=2$ scheme. The security of the channel is assessed at the end of all communication by calculating the visibility of the constructive and destructive interference fringes. The visibility of the interference is given by

\begin{equation}\label{eq:Vis3}
\mathcal{V}=\frac{\mathcal{P}_{0}-\mathcal{P}_{1}}{\mathcal{P}_{0}+\mathcal{P}_{1}}
\end{equation}

\noindent
where $\mathcal{P}_{0}$ ($\mathcal{P}_{1}$) is the summation of counts recorded by detector $D0$ ($D1$) in the shaded time bin. The visibility of interference has a direct impact on the length of the secret key that is distilled from the raw key. For the $d=2$ time-phase scheme discussed here, the quantum bit-error-rate (QBER) of a phase-basis measurement is given by~\cite{scarani2009security}

\begin{equation}\label{eq:QBER}
QBER = \frac{1-\mathcal{V}}{2}.
\end{equation}

\noindent
As the visibility of the interference increases the QBER decreases and therefore a larger secret key can be distilled from the raw key.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1.0\linewidth]{2dTimePhase}
\caption{Time-phase protocol with $d=2$. (a) Temporal and (b) phase basis states. (c) Time-delay interferometer and it's action on a single wavepacket. (d) Interferometer setup for a $d=2$ phase basis measurement. (e) Interference patterns for a measurement on the $\ket{f_{0}}$ state.}
\label{fig:2dTimePhase}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[width=1.0\linewidth]{25um}
\caption{Interference characterization and spatial mode for phase-state measurement. (a) Output from the plus (minus) port of the interferometer shown in blue (black) with Eq.~\ref{eq:Pout} fitted to each curve shown in red. The visibility between the plus and minus port with the phase set at the maximum of the plus port is $98.7\,\%$. (b) Example of the spatial mode under test measured at the output of a $25\,\mu$m core fiber.}
\label{fig:25um}
\end{center}
\end{figure}

Experimentally, the generation of the states shown in Fig.~\ref{fig:2dTimePhase}(a) and (b) are carried out in a similar manner to the method discussed in Ref.~\cite{islam2017robust}. A continuous-wave laser operating at $1550\,$nm  is intensity modulated with electro-optic modulators to create short duration ($\sim80$\,ps) wavepackets. The modulators are driven by a 10\,GHz transceiver from a field-programmable gate array (FPGA). The mean photon number for each wavepacket is controlled using a variable optical attenuator (VOA). The modulation and attenuation are performed in single-mode fiber, then the light is coupled to a multi-mode fiber with a $25\,\mu$m core (Thorlabs FG025LJA) to produce wavepackets with multiple spatial modes. Mode mixing is achieved by stressing the fiber in the same manner as shown in Fig.~\ref{fig:MMode}(a). A $25\,\mu$m core fiber is chosen over the $50\,\mu$m and $105\,\mu$m core fibers mentioned previously because of limitations in our single-photon detection system. The light is sent through the interferometer, where each output is collected with a $25\,\mu$m core fiber. The fibers are then coupled to our single-photon detection system.

Our detection system is comprised of amorphous superconducting nanowire single-photon detectors (SNSPDs) produced by Quantum Opus. SNSPDs are high speed and high efficiency single photon counting devices (see Ref.~\cite{dauler2014review} for a review) that are widely used in quantum communication~\cite{boroson2014overview,islam2017provably} and quantum optics experiments~\cite{shalm2015strong}. Though large area devices~\cite{wollman2017uv} and pixel arrays~\cite{verma2014four,allman2015near} designed to support free-space or multi-mode fiber coupling exist, our devices have an active area  with a diamter of $\sim15\,\mu$m and are designed to couple to a single-mode telecommunications fiber (eg. Corning SMF-28)~\cite{miller2011compact}. However, in this experiment, we couple the $25\,\mu$m core fibers from the output of the interferometer. There will be some amount of aperture loss at the detector due to the mismatch in the size of the detector and the fiber core. We choose to suffer the extra loss at the detector for the purpose of demonstrating multi-mode, single-photon interference. The fringes from the plus and minus port of the interferometer and an example of the spatial mode profile of the beam under test at the output of the fiber are shown in Fig.~\ref{fig:25um}.

\begin{figure}[thb]
\begin{center}
\includegraphics[width=1.0\linewidth]{PhaseMeas}
\caption{Phase-state $\ket{f_{0}}$ preparation and measurement using a $25\,\mu$m core fiber. (a) Single-photon data (blue) with a sum of Gaussian functions fitted to the data (red) for the $\ket{f_{0}}$ state preparation. Zeros above the data denote the phase of each bin, indicating a $\ket{f_{0}}$ state. Time bin boundaries are shown with black dotted lines. (b) Measurement of the $\ket{f_{0}}$ state using the multi-mode interferometer showing the plus $(+)$ output port collected at detector D0 (blue) and minus $(-)$ output port collected at detector D1 (black). Time bin boundaries are shown with dotted lines. The resulting visibility between the bright and dark fringe in the middle time bin is $97.37\pm0.01\,\%$.}
\label{fig:Meas}
\end{center}
\end{figure}

The preparation of the $\ket{f_{0}}$ state is shown with single-photon detection events in Fig.~\ref{fig:Meas}(a). The error rate of the preparation, calculated from the overlap of each wavepacket in to the opposite time-bin, is $0.4\,\%$ and is largely due the limited bandwidth of the modulators used to create the wavepackets. The result of a measurement on the $\ket{f_{0}}$ state using the time-delay interferometer is shown in Fig.~\ref{fig:Meas}(b). As expected, constructive interference is observed in the middle time-bin of detector D0, with destructive interference observed in the middle time-bin of detector D1. Additionally, the heights of the lobes to the left and right of the middle time-bin are approximately $1/4$ the height of the constructive interference fringe, and match well with the predicted measurement outcome. The visibility between the two detectors in the middle time-bin is determined by summing all the counts falling in the time-bin window and using the formula given by Eq.~\ref{eq:Vis3}. The resulting visibility is $97.37\pm 0.01\,\%$. This value is within $1\,\%$ of the visibility reported for a similar phase-state measurement using commercial single-mode time-delay interferometers~\cite{islam2017robust}.

\section{Conclusion}

In this article, we detail the design methodology, performance simulations, construction and characterization of a time-delay interferometer with a 5\,GHz free-spectral-range (200\,ps time-delay). The interferometer is set-up in a Michelson-type layout with glass beam paths in each arm. A wide field-of-view and excellent thermal stability can be achieved with a proper material choice and optical design, while satisfying the target time-delay. The performance results are highlighted by a maximum single-mode visibility of $\mathcal{V}_{\mathrm{SM}}=99.02\pm0.05\,\%$, comparable to commercial single-mode products, and multi-mode visibility of $\mathcal{V}_{\mathrm{MM}}=98.38\pm0.01\,\%$, several percent higher than recent reports. Additionally, multi-mode visibility $\mathcal{V}_{\mathrm{MM}}\approx98\,\%$ is maintained independent of spatial mode profile for both $50\,\mu$m and $105\,\mu$m core fibers. High-quality interference is demonstrated over a modest temperature range of $\pm1^{\,\circ}$C where the path-length shift of the interferometer is 130\,nm/$^{\,\circ}$C, enabling good thermal stability with modest temperature control.

Finally, we demonstrate the application of our interferometer to a time-phase QKD protocol, where the presence of an eavesdropper is quantified by the interference visibility between the two output ports of a time-delay interferometer. The high visibility of our interferometer enables phase-state measurements with a low QBER, allowing efficient extraction of a secret key. We demonstrate a phase-state measurement using weak coherent states with visibility of $97.37\pm 0.01\,\%$ for $1550\,$nm photons from a $25\,\mu$m core fiber.

Further improvements can be made to the design and construction of the interferometer to achieve better performance.
In our estimation the largest contributing factor to the degradation of the multi-mode interference visibility is an error in the length of the glass rods during polishing. Each piece of glass is undersized by $400-600\,\mu$m. The implication of this error is that the actual dimensions of the glass no longer satisfy the field-widening condition and therefore result in degraded interference, especially for multi-mode beam profiles. We were able to compensate for this error by a coarse adjustment of the position of one of the reflecting mirrors. Using a more precise method for measuring the position of the mirror, we expect the quality of multi-mode visbility could be improved to over $99\,\%$.

\begin{acknowledgments}
This work is supported by the Office of Naval Research (ONR) (N00014-13-1-0627) and the National Aeronautics and Space Administration (NASA) (NNX13AP35A).
\end{acknowledgments}


\section{Introduction}

Upsampling is an essential stage for most dense prediction tasks using deep convolutional neural networks (CNNs). The frequently used upsampling operators include transposed convolution~\cite{zeiler2014visualizing,long2015fully}, unpooling~\cite{badrinarayanan2017segnet}, \textit{periodic shuffling}~\cite{shi2016real} (also known as depth-to-space), and naive interpolation~\cite{lin2017refine,chen18v3} followed by convolution. These operators, however, are not general-purpose designs and often have different behaviors in different tasks.

The widely-adopted operator in semantic segmentation or depth estimation is bilinear interpolation, rather than unpooling. A reason is that the feature map generated by unpooling is too sparse, while bilinear interpolation is likely to generate the feature map that depicts semantically-consistent regions. This is particularly true for semantic segmentation and depth estimation where pixels in a region often share the same class label or have similar depth. However, bilinear interpolation performs much worse than unpooling in boundary-sensitive tasks such as image matting. A fact is that the leading deep image matting model~\cite{xu2017deep} largely borrows the design from the SegNet~\cite{badrinarayanan2017segnet}, where unpooling is introduced. When adapting other state-of-the-art segmentation models, such as DeepLabv3+~\cite{chen18v3} and RefineNet~\cite{lin2017refine}, to this task, unfortunately, we observe both DeepLabv3+ and RefineNet fail to recover boundary details (Fig.~\ref{fig:alpha_mattes}), compared to SegNet. This makes us to ponder over \textit{what is missing} in these encoder-decoder models. After making a thorough comparison between different architectures and conducting ablative studies (Section~\ref{subsec:ablation_study}), the answer is finally made clear---\textit{indices matter}.

\begin{figure}[t]
	\captionsetup{font=small,singlelinecheck=true}
	\begin{center}
		\includegraphics[width=\linewidth]{mattes_intro.pdf}
	\end{center}
	\caption{Alpha mattes of different models. From left to right, Deeplabv3+~\cite{chen18v3}, RefineNet~\cite{lin2017refine}, Deep Matting~\cite{xu2017deep} and Ours. %
	Bilinear upsampling fails to recover subtle details, but unpooling and our learned upsampling operator can produce much clear mattes with good local contrast.}
	\label{fig:alpha_mattes}
\end{figure}

Compared to the bilinearly upsampled feature map, unpooling uses max-pooling indices to guide upsampling. Since boundaries in the shallow layers usually have the maximum responses, indices extracted from these responses record the boundary locations. The feature map projected by the indices thus shows improved boundary delineation. Above analyses reveal a fact that, different upsampling operators have different characteristics, and we expect a specific behavior of the upsampling operator when dealing with specific image content in a certain visual task.

It would be interesting to pose the question: \textit{Can we design a generic operator to upsample feature maps that better predict boundaries and regions simultaneously?} A key observation of this work is that max unpooling, bilinear interpolation or other \textit{upsampling operators are some forms of index functions}. For example, the nearest neighbor interpolation of a point is equivalent to allocating indices of one to its neighbor and then map the value of the point. In this sense, indices are models~\cite{kraska2018case}, therefore indices can be modeled and learned. In this work, we model indices as a function of the local feature map and learn an index function to perform upsampling within deep CNNs. In particular, we present a novel index-guided encoder-decoder framework, which naturally generalizes SegNet. Instead of using  max-pooling and unpooling, we introduce \mbox{indexed pooling} and indexed upsampling operators where downsampling and upsampling are guided by learned indices. The indices are generated dynamically conditioned on the feature map and are learned using a fully convolutional network, termed IndexNet, without supervision. IndexNet is a highly flexible module, which can be used as a plug-in applying to any off-the-shelf convolutional networks that have coupled downsampling and upsampling stages. Compared to the fixed $\max$ function, learned index functions show potentials for simultaneous boundary and region delineation.

We demonstrate the effectiveness of IndexNet on natural image matting as well as other visual tasks. In image matting, the quality of learned indices can be visually observed from predicted alpha mattes. By visualizing learned indices, we show that the indices automatically learn to capture the boundaries and textural patterns.
We further investigate alternative ways to design IndexNet, and show through extensive experiments that IndexNet can effectively improve deep image matting  both qualitatively and quantitatively. In particular, we observe that our best MobileNetv2-based~\cite{sandler2018mobilenetv2} model exhibits at least $16.1\%$ improvement against the previous best deep model, i.e., the VGG-16-based model in \cite{xu2017deep}, on the Composition-1k matting dataset.  We achieve this with using less training data, and a much more compact model, therefore significantly faster inference speed.

\section{Related Work}

We review existing widely-used upsampling operators and the main application of IndexNet---deep image matting.

\paragraph{Upsampling in Deep Networks}
Upsampling is an essential stage for almost all dense prediction tasks.
It has been intensively studied about \textit{what is the principal way to recover the resolution of the downsampled feature map} (decoding). The \textit{deconvolution} operator, also known as transposed convolution, was initially used  in~\cite{zeiler2014visualizing} to visualize convolutional activations and latter introduced to semantic segmentation \cite{long2015fully}.
To avoid checkerboard artifacts, a follow-up suggestion is the ``resize+convolution'' paradigm, which has currently become the standard configuration in state-of-the-art semantic segmentation models~\cite{chen18v3,lin2017refine}. Aside from these, \textit{perforate}~\cite{osendorfer2014image} and \textit{unpooling}~\cite{badrinarayanan2017segnet} are also two operators that generate sparse indices to guide upsampling. The indices are able to capture and keep boundary information, but the problem is that two operators induce sparsity after upsampling. Convolutional layers with large filter sizes must follow for densification. In addition, \textit{periodic shuffling} ($\mathcal{PS}$) was introduced in~\cite{shi2016real} as a fast and memory-efficient upsampling  operator for image super-resolution. $\mathcal{PS}$ recovers resolution by rearranging the feature map of size $H\times W\times Cr^2$ to $rH\times rW\times C$.

Our work is primarily inspired by the unpooling operator~\cite{badrinarayanan2017segnet}. We remark that, it is important to keep  the spatial information before loss of such information occurred in feature map downsampling, and more importantly, to use stored information during upsampling. Unpooling shows a simple and effective case of doing this, but we argue there is much room to improve. In this paper, we illustrate that the unpooling operator is a special form of index function, and we can learn an index function beyond unpooling.

\begin{figure*}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=.879\textwidth,angle=0]{indexnet_framework.pdf}

	\caption{Index-guided encoder-decoder framework. The proposed \mbox{IndexNet} dynamically predicts indices for individual local regions, conditional on the input local feature map itself. The predicted indices are further utilized to guide the downsampling in the encoding stage and the upsampling in corresponding decoding stage.}
	\label{fig:framework}
\end{figure*}

\paragraph{Deep Image Matting}

In the past decades, image matting methods have been extensively studied from a low-level view~\cite{aksoy2017designing,chen2013knn,chen2013image,chuang2001bayesian,guan2006easy,he2011global,lee2011nonlocal,levin2008closed,sun2004poisson}; and particularly, they have been designed to solve the matting equation. Despite being theoretically elegant, these methods heavily rely on the color cues, rendering failures of matting in general natural scenes where colors cannot be used as reliable cues.

With the tremendous success of deep CNNs in high-level vision tasks~\cite{girshick2014rich,krizhevsky2012imagenet,long2015fully}, deep matting methods are emerging. Some initial attempts appeared in~\cite{cho2016natural} and~\cite{shen2016deep}, where classic matting approaches, such as closed-form matting~\cite{levin2008closed} and KNN matting~\cite{chen2013knn}, are still used as the backends in deep networks. Although the networks are trained end-to-end and can extract powerful features, the final performance is limited by the conventional backends. These attempts may be thought as semi-deep matting. Recently fully-deep image matting was proposed~\cite{xu2017deep}. In \cite{xu2017deep} the authors presented the first deep image matting approach based on SegNet~\cite{badrinarayanan2017segnet} and significantly outperformed other competitors. Interestingly, this SegNet-based architecture becomes the standard configuration in many recent deep matting methods~\cite{chen2018tom,chen2018human,wang2018propagation}.

SegNet is effective in matting but also computation-expensive and memory-inefficient. For instance, the inference can only be executed on CPU when testing high-resolution images, which is practically unattractive.
We show that, with our proposed IndexNet, even a lightweight backbone such as MobileNetv2-based model can surpass the VGG-16 based method in \cite{xu2017deep}.

\section{An Indexing Perspective of Upsampling}

With the argument that upsampling operators are index functions, here we offer an unified index perspective of upsampling operators. The unpooling operator is straightforward. We can define its index function in a $k\times k$ local region as an indicator function
\begin{equation}
I_{max}(x)=\mathbbm{1}(x=\max (\mat X))\,,x\in\mat X \,,
\end{equation}
where $\mat X\in\mathbb{R}^{k\times k}$. Similarly, if one extracts indices from the average pooling operator, the index function takes the form
\begin{equation}
I_{avg}(x)=\mathbbm{1}(x\in\mat X)\,.
\end{equation}
If further using $I_{avg}(x)$ during upsampling, it is equivalent to the nearest neighbor interpolation. Regarding the bilinear interpolation and deconvolution operators, their index functions have an identical form
\begin{equation}
I_{bilinear/dconv}(x)=\mat W\otimes\mathbbm{1}(x\in\mat X)\,,
\end{equation}
where $\mat W$ is the weight/filter of the same size as $\mat X$, and $\otimes$ denotes the element-wise multiplication. The difference is that, $\mat W$ in deconvolution is learned, while $\mat W$ in bilinear interpolation stays fixed. Indeed, bilinear upsampling has been shown to be a special case of deconvolution~\cite{long2015fully}. Notice that, in this case, the index function generates soft indices. The sense of index for the $\mathcal{PS}$ operator~\cite{shi2016real} is even much clear, because the rearrangement of the feature map per se is an indexing process. Considering $\mathcal{PS}$ a tensor $\mathcal{Z}$ of size $1\times1\times r^2$ to a matrix $\mat Z$ of size $r\times r$, the index function can be expressed by the one-hot encoding
\begin{equation}
I^l_{ps}(x)=\mathbbm{1}(x= \mathcal{Z}_l)\,,l=1,...,r^2\,,
\end{equation}
such that $\mat Z_{m,n}=\mathcal{Z}[I^l_{ps}(x)]$, where $m=1,...,r$, $n=1,...,r$, and $l=(r-1)*m+n$. $\mathcal{Z}_l$ denotes the $l$-th element of $\mathcal{Z}$. A similar notation applies to $\mat Z_{m,n}$.

Since upsampling operators can be unified by the notion of index function, in theory it is possible to learn an index function that adaptively captures local spatial patterns.

\section{\mbox{Index-Guided Encoder-Decoder Framework}}

Our framework is a natural generalization of SegNet, as schematically illustrated in Fig.~\ref{fig:framework}. For ease of exposition, we assume the downsampling and upsampling rates are $2$, and the pooling operator has a kernel size of $2\times2$. At the core of our framework is the IndexNet module that dynamically generates indices given the feature map. The proposed indexed pooling and indexed upsampling operators further receive generated indices to guide the downsampling and upsampling, respectively. In practice, multiple such modules can be combined and used analogues to the max pooling layers. We provide details as follows.

\subsection{Learning to Index, to Pool, and to Upsample}

\noindent\textbf{IndexNet} models the index as a function of the feature map $\mathcal{X}\in\mathbb{R}^{H\times W\times C}$. It generates two index maps for downsampling and upsampling given the input $\mathcal{X}$. An important concept for the index is that an index can either be represented in a natural order, e.g., 1, 2, 3, ..., or be represented in a logical form, i.e., 0, 1, 0, ..., which means an index map can be used as a mask. In fact, this is how we use the index map in downsampling and upsampling. The predicted index shares the same physical notation of the index in computer science, except that we generate \textit{soft} indices for smooth optimization, i.e., for any index $i$, $i\in[0,1]$.

IndexNet consists of a predefined index block and two index normalization layers. An index block can  simply be a heuristically defined function, e.g., a $\max$  function, or more generally, a neural network. In this work, the index block is designed to use a fully convolutional network. According to the shape of the output index map, we investigate two families of index networks: \textit{holistic index networks (HINs)} and \textit{depthwise (separable) index networks (DINs)}. Their conceptual differences are shown in Fig.~\ref{fig:holistic_depthwise}. HINs learn an index function $I(\mathcal{X}):\mathbb{R}^{H\times W\times C}\rightarrow\mathbb{R}^{H\times W\times1}$. In this case, all channels of the feature map share a holistic index map. In contrast, DINs learn an index function $I(\mathcal{X}):\mathbb{R}^{H\times W\times C}\rightarrow\mathbb{R}^{H\times W \times C}$, where the index map is of the same size as the feature map. We will discuss concrete design of index networks in Sections~\ref{subsec:holistic_networks} and~\ref{subsec:depthwise_networks}.

Note that the index map sent to the encoder and decoder are normalized differently. The decoder index map only goes through a \textit{sigmoid} function such that for any predicted index $i\in(0,1)$. As for the encoder index map, indices of a local region $L$ are further normalized by a \textit{softmax} function such that $\sum_{i\in{L}}i=1$. The reason behind the second normalization is to guarantee the magnitude consistency of the feature map after downsampling.

\vspace{5pt}
\noindent\textbf{Indexed Pooling} ($\mathcal{IP}$) executes downsampling using generated indices. Given a local region $E\in\mathbb{R}^{k\times k}$, $\mathcal{IP}$ calculates a weighted sum of activations and corresponding indices over $E$ as $\mathcal{IP}(E)=\sum_{x\in E}I(x)x$, where $I(x)$ is the index of $x$. It is easy to infer that max pooling and average pooling are both special cases of $\mathcal{IP}$. In practice, this operator can be easily implemented with an element-wise multiplication between the feature map and the index map, an average pooling layer, and a multiplication of a constant, as instantiated in Fig.~\ref{fig:framework}.

\vspace{5pt}
\noindent\textbf{Indexed Upsampling} ($\mathcal{IU}$) is the inverse operator of $\mathcal{IP}$. $\mathcal{IU}$ upsamples $d\in\mathbb{R}^{1\times 1}$ that spatially corresponds to $E$ taking the same indices into account. Let $I\in\mathbb{R}^{k\times k}$ be the local index map formed by $I(x)$s, $\mathcal{IU}$ upsamples $d$ as $\mathcal{IU}(d)=I\otimes D$, where $\otimes$ denotes the element-wise multiplication, and $D$ is of the same size as $I$ and is upsampled from $d$ with the nearest neighbor interpolation. An important difference between deconvolution and $\mathcal{IU}$ is that, deconvolution applies a fixed kernel to all local regions, even if the kernel is learned, while $\mathcal{IU}$ upsamples different regions with different kernels (indices).

\begin{figure}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=3in,angle=0]{holistic_depthwise.pdf}

	\caption{Conceptual differences between holistic index and depthwise index.}
	\label{fig:holistic_depthwise}
\end{figure}

\subsection{Holistic Index Networks}\label{subsec:holistic_networks}

Here we instantiate two types of HINs. Recall that HINs learn an index function $I(\mathcal{X}):\mathbb{R}^{H\times W\times C}\rightarrow\mathbb{R}^{H\times W\times1}$. A naive design choice is to assume a linear relationship between the feature map and the index map.

\vspace{5pt}
\noindent\textit{Linear Holistic Index Networks.} An example is shown in Fig.~\ref{fig:holistic_networks}(a). The network is implemented in a fully convolutional way. It first applies $2$-stride $2\times2$ convolution to the feature map of size $H\times W\times C$, generating a concatenated index map of size $H/2\times W/2\times 4$. Each slice of the index map ($H/2\times W/2\times 1$) is designed to correspond to the indices of a certain position of all local regions, e.g., the top-left corner of all $2\times 2$ regions. The network finally applies a $\mathcal{PS}$-like shuffling operator to rearrange the index map to the size of $H\times W\times 1$.

In many situations, assuming a linear relationship is not sufficient. An obvious fact is that a linear function even cannot fit the $\max$ function. Naturally the second design choice is to add nonlinearity into the network.

\vspace{5pt}
\noindent\textit{Nonlinear Holistic Index Networks.} Fig.~\ref{fig:holistic_networks}(b) illustrates a nonlinear HIN where the feature map is first projected to a map of size $H/2\times W/2\times 2C$, followed by a batch normalization layer and a ReLU function for nonlinear mappings. We then use point-wise convolution to reduce the channel dimension to an indices-compatible size. The rest transformations follow its linear counterpart.

\vspace{5pt}
\noindent\textbf{Remark 1}. Note that, the holistic index map is shared by all channels of the feature map, which means the index map should be expanded to the size of $H\times W\times C$ when feeding into $\mathcal{IP}$ and $\mathcal{IU}$. Fortunately, many existing packages support implicit expansion over the singleton dimension. This index map could be thought as a collection of local attention maps~\cite{mnih2014recurrent} applied to individual local spatial regions. In this case, the $\mathcal{IP}$ and $\mathcal{IU}$ operators can also be referred to ``attentional pooling'' and ``attentional upsampling''.

\begin{figure}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=2.4in,angle=0]{holistic_index_networks.pdf}\vspace{-10pt}
	\caption{Holistic index networks. (a) a linear index network; (b) a nonlinear index network.}
	\label{fig:holistic_networks}
\end{figure}

\begin{figure}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=3.2in,angle=0]{depthwise_index_networks.pdf}\vspace{-8pt}
	\caption{Depthwise index networks. $N=C$ for the O2O assumption, and $N=1$ for the M2O. The masked modules are invisible to linear networks.}
	\label{fig:depthwise_networks}
\end{figure}

\subsection{Depthwise Index Networks}\label{subsec:depthwise_networks}

In DINs, we find $I(\mathcal{X}):\mathbb{R}^{H\times W\times C}\rightarrow\mathbb{R}^{H\times W \times C}$, i.e., each spatial index corresponds to each spatial activation. This family of networks further has two high-level design strategies that correspond to two different assumptions.

\vspace{3pt}
\noindent\textbf{One-to-One (O2O) Assumption} assumes that each slice of the index map only relates to its corresponding slice of the feature map. It can be denoted by a local index function $l(\mathcal{X}): \mathbb{R}^{k\times k\times1}\rightarrow\mathbb{R}^{k\times k\times1}$, where $k$ denotes the size of local region. Similar to HINs, DINs can also be designed to have linear/nonlinear modeling ability. Fig.~\ref{fig:depthwise_networks} shows an example when $k=2$. Note that, different from HINs, DINs follow a multi-column architecture. Each column predicts indices specific to a certain spatial location of all local regions. The O2O assumption can be easily satisfied in DINs with grouped convolution.

\vspace{3pt}
\noindent\textit{Linear Depthwise Index Networks.}  As per Fig.~\ref{fig:depthwise_networks}, a feature map goes through four parallel convolutional layers with the same kernel size of $2\times2\times C$, a stride of $2$, and $C$ groups, leading to four downsampled feature maps of size $H/2\times W/2 \times C$. The final index map is composed from the four feature maps by shuffling and rearrangement. Note that the parameters of four convolutional layers are not shared.

\vspace{3pt}
\noindent\textit{Nonlinear Depthwise Index Networks.} Nonlinear DINs can be easily modified from linear DINs by inserting four extra convolutional layers. Each of them is followed by a BN layer and a ReLU unit, as shown in Fig.~\ref{fig:depthwise_networks}. The rest remains the same as the linear DINs.

\vspace{3pt}
\noindent\textbf{Many-to-One (M2O) Assumption} assumes that each slice of the index map relates with all channels of the feature map. The local index function is defined as $l(\mathcal{X}): \mathbb{R}^{k\times k\times C}\rightarrow\mathbb{R}^{k\times k\times1}$. Compared to O2O DINs, the only difference in implementation is the use of standard convolution instead of group convolution, i.e., $N=1$ in Fig.~\ref{fig:depthwise_networks}.

\vspace{3pt}
\noindent\textbf{Learning with Weak Context}. A desirable property of IndexNet is that it can predict indices even from a large local feature map, e.g., $l(\mathcal{X}): \mathbb{R}^{2k\times 2k\times C}\rightarrow\mathbb{R}^{k\times k\times1}$. An intuition behind this idea is that, if one identifies a local maximum point from a $k\times k$ region, its surrounding $2k\times2k$ region can further support whether this point is a part of a boundary or just an isolated noise point. This idea can be easily implemented by enlarging the convolutional kernel and is also applicable to HINs.

\vspace{3pt}
\noindent\textbf{Remark 2}. Both HINs and DINs have merits and drawbacks. It is clear that DINs have higher capacity than HINs, so DINs may capture more complex local patterns but also be at a risk of overfitting. By contrast, the index map generated by HINs is shared by all channels of the feature map, so the decoder feature map can reserve its expressibility without forcibly reducing its dimensionality to fit the shape of the index map during upsampling. This gives much flexibility for decoder design, while it is not the case for DINs.

\subsection{Relation to Other Networks}

If considering the dynamic property of IndexNet, \mbox{IndexNet} shares a similar spirit with some recent networks.

\vspace{3pt}
\noindent\textbf{Spatial Transformer Networks (STNs)}~\cite{jaderberg2015spatial}. The STN learns dynamic spatial transformation by regressing desired transformation parameters $\theta$ with a localized network. A spatially-transformed output is then produced by a sampler parameterized by $\theta$. Such a transformation is holistic for the feature map, which is similar to HINs. The differences between STN and IndexNet are that their learning targets have different physical definitions (spatial transformations vs. spatial indices), and that, STN is designed for global transformation, while IndexNet predicts local indices.

\vspace{3pt}
\noindent\textbf{Dynamic Filter Networks (DFNs)}~\cite{jia2016dynamic}. The DFN dynamically generates filter parameters on-the-fly with a so-called filter generating network. Compared to conventional filter parameters that are initialized, learned, and stayed fixed during inference, filter parameters in DFN are dynamic and sample-specific. The main difference between DFN and IndexNet lies in the motivation of the design. Dynamic filters are learned for adaptive feature extraction, but learned indices are used for dynamic downsampling and upsampling.

\vspace{3pt}
\noindent\textbf{Deformable Convolutional Networks (DCNs)}~\cite{dai2017deformable}. The DCN introduces deformable convolution and deformable RoI pooling. The key idea is to predict offsets for convolutional and pooling kernels, so DCN is also a dynamic network. While these convolution and pooling operators concern spatial transformations, they are still built upon standard max pooling and are not designed for upsampling purposes. By contrast, index-guided $\mathcal{IP}$ and $\mathcal{IU}$ are fundamental operators and may be integrated into RoI pooling.

\vspace{3pt}
\noindent\textbf{Attention Networks}~\cite{mnih2014recurrent}. Attention networks are a broad family of networks that adopt attention mechanisms. The mechanisms introduce multiplicative interactions between inferred attention maps and feature maps. In Computer Vision, these mechanisms often refer to spatial attention~\cite{wang2017residual}, channel attention~\cite{hu2018squeeze} or both~\cite{woo2018cbam}. As aforementioned, $\mathcal{IP}$ and $\mathcal{IU}$ in HINs can be viewed as attentional operators to some extent, which means indices are attention. In a reverse sense, attention is also indices. For example, max-pooling indices are a form of hard attention. Indices offer a new perspective to understand attention. It is worth noting that, despite IndexNet in its current implementation closely relates to attention, it has a distinct physical definition and specializes in upsampling rather than refining feature maps.

\section{Results and Discussions}

We evaluate our framework and IndexNet on the task of image matting. This task is particularly suitable for visualizing the quality of learned indices. We mainly conduct experiments on the Adobe Image Matting dataset~\cite{xu2017deep}. This is so far the largest publicly available matting dataset. The training set has 431 foreground objects and ground-truth alpha mattes.\footnote{The original paper reported that there were 491 images, but the released dataset only includes 431 images.
As a result, we
use fewer
training data
than the original paper.}
Each foreground is composited with 100 background images randomly chosen from MS COCO~\cite{lin2014microsoft}. The test set termed Composition-1k includes 100 unique objects. Each of them is composited with 10 background images chosen from Pascal VOC~\cite{everingham2010pascal}. Overall, we have 43100 training images and 1000 testing images. We evaluate the results using widely-used Sum of Absolute Differences (SAD), Mean Squared Error (MSE), and perceptually-motivated Gradient (Grad) and Connectivity (Conn) errors~\cite{rhemann2009perceptually}. The evaluation code implemented by~\cite{xu2017deep} is used. In what follows, we first describe our modified MobileNetv2-based architecture and training details. We then perform extensive ablation studies to justify choices of model design, make comparisons of different index networks, and visualize learned indices. We also report performance on the $\tt alphamatting.com$ online benchmark~\cite{rhemann2009perceptually} and extend IndexNet to other visual tasks.

\subsection{Implementation Details}

Our implementation is based on PyTorch~\cite{paszke2017automatic}. Here we describe the network architecture used and some essential training details.
\begin{figure}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=2.8in,angle=0]{network_architecture.pdf}\vspace{-5pt}
	\caption{Customized MobileNetv2-based encoder-decoder network architecture. Our modifications are boldfaced.}
	\label{fig:network_architecture}
\end{figure}
\begin{table*}[!t] \small
	\captionsetup{font=small,singlelinecheck=true}
	\centering
	\addtolength{\tabcolsep}{3pt}
	\renewcommand\arraystretch{1.0}
	\begin{tabular}{llccccc|ccccc}
			\hline
			No. & Architecture					& Backbone		& Fusion 	& Indices 	& Context 	& OS 	& SAD 	& MSE 	& Grad & Conn \\
			\hline
			B1	& DeepLabv3+~\cite{chen18v3}	& MobileNetv2	& Concat	& No		& ASPP		& 16	& 60.0	& 0.020	& 39.9 & 61.3 \\
			B2 	& RefineNet~\cite{lin2017refine}& MobileNetv2	& Skip 		& No 		& CRP 		& 32 	& 60.2	& 0.020	& 41.6 & 61.4 \\
			\rowcolor{mygray}
			B3 	& SegNet~\cite{xu2017deep}		& VGG16			& No		& Yes		& No		& 32	& \textbf{54.6} 	& \textbf{0.017} & 36.7 & 55.3 \\
			\rowcolor{mygray}
			B4 	& SegNet						& VGG16			& No		& No		& No		& 32	& 122.4 & 0.100 & 161.2 & 130.1 \\
			\rowcolor{mygray}
			B5 	& SegNet		 				& MobileNetv2	& No 		& Yes 		& No 		& 32	& 60.7	& 0.021	& 40.0 & 61.9 \\
			\rowcolor{mygray}
			B6 	& SegNet		 				& MobileNetv2	& No 		& No 		& No 		& 32	& 78.6  & 0.031 & 101.6& 82.5 \\
			B7 	& SegNet		 				& MobileNetv2	& No 		& Yes 		& ASPP 		& 32	& 58.0	& 0.021	& 39.0 & 59.5 \\
			B8 	& SegNet		 				& MobileNetv2	& Skip 		& Yes 		& No 		& 32	& 57.1 	& 0.019	& 36.7 & 57.0 \\
			B9 	& SegNet		 				& MobileNetv2	& Skip 		& Yes 		& ASPP 		& 32	& 56.0	& \textbf{0.017}	& 38.9 & 55.9 \\
			B10 & UNet	 						& MobileNetv2	& Concat 	& Yes 		& No 		& 32 	& 54.7	& \textbf{0.017}	& 34.3 & \textbf{54.7} \\
			B11 & UNet	 						& MobileNetv2	& Concat 	& Yes 		& ASPP 		& 32	& 54.9	& \textbf{0.017} & \textbf{33.8} & 55.2 \\
			%
		%
			\hline
		\end{tabular}
	%
	\vspace{-8pt}
	\caption{Ablation study of design choices. Fusion: fuse encoder features; Indices: max-pooling indices (when Indices is `No', bilinear interpolation is used for upsampling); CRP: chained residual pooling~\cite{lin2017refine}; ASPP: atrous spatial pyramid pooling~\cite{chen18v3}; OS: output stride.
	%
	The lowest errors are boldfaced.}
	\label{tab:architecture}
	%
\end{table*}

\vspace{3pt}
\noindent\textbf{Network Architecture}. We build our model based on MobileNetv2~\cite{sandler2018mobilenetv2} with only slight modifications to the backbone. An important reason why we choose MobileNetv2 is that this lightweight model allows us to infer high-resolution images on a GPU, while other high-capacity backbones cannot. The basic network configuration is shown in Fig.~\ref{fig:network_architecture}. It also follows the encoder-decoder paradigm same as SegNet. We simply change all 2-stride convolution to be 1-stride and attach 2-stride $2\times2$ max pooling after each encoding stage for downsampling, which allows us to extract indices. If applying the IndexNet idea, max pooling and unpooling layers can be replaced with $\mathcal{IP}$ and $\mathcal{IU}$, respectively. We also investigate alternative ways for low-level feature fusion and whether encoding context (Section~\ref{subsec:ablation_study}). Notice that, the matting refinement stage~\cite{xu2017deep} is not considered in this paper.

\vspace{3pt}
\noindent\textbf{Training Details}. To enable a direct comparison with deep matting~\cite{xu2017deep}, we follow the same training configurations used in~\cite{xu2017deep}. The 4-channel input concatenates the RGB image and its trimap. We follow exactly the same data augmentation strategies, including $320\times320$ random cropping, random flipping, random scaling, and random trimap dilation. All training samples are created on-the-fly. We use a combination of the alpha prediction loss and the composition loss during training as in~\cite{xu2017deep}. Only losses from the unknown region of the trimap are calculated. Encoder parameters are pretrained on ImageNet~\cite{deng2009imagenet}. Note that, the parameters of the $4$-th input channel are initialized with zeros. All other parameters are initialized with the improved Xavier~\cite{he2015delving}. The Adam optimizer~\cite{kingma2015adam} is used. We update parameters with $30$ epochs (around $90,000$ iterations). The learning rate is initially set to $0.01$ and reduced by $10\times$ at the $20$-th and $26$-th epoch respectively. We use a batch size of $16$ and fix the BN layers of the backbone.

\subsection{Adobe Image Matting Dataset}
\label{subsec:ablation_study}

\vspace{3pt}
\noindent\textbf{Ablation Study on Model Design}. Here we investigate strategies for fusing low-level features (no fusion, skip fusion as in ResNet~\cite{he2016deep} or concatenation as in UNet~\cite{ronneberger2015u}) and whether encoding context for image matting. $11$ baselines are consequently built to justify model design. Results on the Composition-1k testing set are reported in Table~\ref{tab:architecture}. B3 is cited from~\cite{xu2017deep}. We can make the following observations: \textbf{i}) Indices are of great importance. Matting can significantly benefit from only indices (B3 vs.\ B4, B5 vs.\ B6); \textbf{ii}) State-of-the-art semantic segmentation models cannot be directly applied to image matting (B1/B2 vs.\  B3); \textbf{iii}) Fusing low-level features help, and concatenation works better than the skip connection but at a cost of increased computation (B5 vs.\ B8 vs.\ B10 or B7 vs.\ B9 vs.\ B11); \textbf{iv}) Our intuition tells that the context may not help a low-level task like matting, while results show that encoding context is generally encouraged (B5 vs.\  B7 or B8 vs.\ B9 or B10 vs.\  B11). Indeed, we observe that the context sometimes can help to improve the quality of the background; \textbf{v}) A MobileNetv2-based model can work as well as a VGG-16-based one with appropriate design choices (B3 vs.\ B11).

For the following experiments, we now mainly use B11.

\vspace{3pt}
\noindent\textbf{Ablation Study on Index Networks}. Here we compare different index networks and justify their effectiveness. The configurations of index networks used in the experiments follow Figs.~\ref{fig:holistic_networks} and~\ref{fig:depthwise_networks}. We primarily investigate the $2\times2$ kernel with a stride of $2$. Whenever the weak context is considered, we use a $4\times4$ kernel in the first convolutional layer of index networks. To highlight the effectiveness of HINs, we further build a baseline called \textit{holistic max index (HMI)} where max-pooling indices are extracted from a squeezed feature map $\mathcal{X}'\in\mathbb{R}^{H\times W\times1}$. $\mathcal{X}'$ is generated by applying the max function along the channel dimension of $\mathcal{X}\in\mathbb{R}^{H\times W\times C}$. We also report the performance when setting the width multiplier of MobileNetV2 used in B11 to be $1.4$ (B11-1.4). This allows us to justify whether the improved performance is due to increased model capacity.
Results on the \mbox{Composition-1k} testing dataset are listed in Table~\ref{tab:index_function}. We observe that, except the most naive linear HIN, all index networks consistently reduce the errors. In particular, nonlinearity and the context generally have a positive effect on deep image matting. Compared to HMI, the direct baseline of HINs, the best HIN (``Nonlinear+Context'') has at least $12.3\%$ relative improvement. Compared to B11, the baseline of DINs, M2O DIN with ``Nonlinear+Context'' exhibits at least $16.5\%$ relative improvement. Notice that, our best model even outperforms the state-of-the-art DeepMatting~\cite{xu2017deep} that has the refinement stage, and is also computationally efficient with less memory consumption---the inference can be performed on the GTX 1070 over $1920\times1080$ high-resolution images. Some qualitative results are shown in Fig.~\ref{fig:composition-1k-results}. Our predicted mattes show improved delineation for edges and textures like hair and water drops.

\begin{table}[!t] \footnotesize
	\captionsetup{font=small,singlelinecheck=true}
	\centering
	\renewcommand\arraystretch{1.0}
	\begin{tabular}{>{\centering}p{.3cm}|>{\centering}p{.3cm}|c|c|c|c|c|c}
			\hline
			\multicolumn{2}{c|}{Method} 				& \#Param. & GFLOPs & SAD & MSE & Grad & Conn\\
			\hline
			\multicolumn{2}{l|}{B3~\cite{xu2017deep}} 	& 130.55M & 32.34 & 54.6 & 0.017 & 36.7 & 55.3 \\
			\multicolumn{2}{l|}{B11} 					& 3.75M	  & 4.08  & 54.9 & 0.017 & 33.8 & 55.2 \\
			\multicolumn{2}{l|}{B11-1.4} 				& 8.86M	  & 7.61  & 55.6 & 0.016 & 36.4 & 55.7 \\
			\multicolumn{2}{l|}{HMI} 					& 3.75M   & 4.08  & 56.5 & 0.021 & 33.0 & 56.4 \\
			\hline
			NL & C & $\Delta$ \\
			\hline
			& & \multicolumn{6}{c}{HINs}\\
			\hline
			& 									 		& +4.99K  & 4.09  & 55.1 & 0.018 & 32.1 & 55.2\\
			& \checkmark                                & +19.97K & 4.11  & 53.5 & 0.018 & 31.0 & 53.5\\
			\checkmark & 						 		& +0.26M  & 4.22  & 50.6 & 0.015 & 27.9 & 49.4\\
			\checkmark & \checkmark 			 		& +1.04M  & 4.61  & 49.5 & 0.015 & \textbf{25.6} & 49.2\\
			\hline
			& & \multicolumn{6}{c}{O2O DINs}\\
			\hline
			&          					                & +4.99K  & 4.09  & 50.3 & 0.015 & 33.7 & 50.0\\
			& \checkmark                                & +19.97K & 4.11  & 47.8 & 0.015 & 26.9 & 45.6\\
			\checkmark &            	                & +17.47K & 4.10  & 50.6 & 0.016 & 26.5 & 50.3\\
			\checkmark & \checkmark 	                & +47.42K & 4.15  & 50.2 & 0.016 & 26.8 & 49.3\\
			\hline
			& & \multicolumn{6}{c}{M2O DINs}\\
			\hline
			&          				             		& +0.52M & 4.34   & 51.0 & 0.015 & 33.7 & 50.5\\
			& \checkmark                                & +2.07M & 5.12   & 50.6 & 0.016 & 31.9 & 50.2\\
			\checkmark &                         		& +1.30M & 4.73   & 48.9 & 0.015 & 32.1 & 47.9\\
			\checkmark & \checkmark              		& +4.40M & 6.30   & \textbf{45.8} & \textbf{0.013} & 25.9 & \textbf{43.7}\\
			\hline
			\hline
			\multicolumn{4}{l|}{Closed-Form~\cite{levin2008closed}} & 168.1 & 0.091 & 126.9 & 167.9 \\
			\multicolumn{4}{l|}{DeepMatting w. Refinement~\cite{xu2017deep}} & 50.4 & 0.014 & 31.0 & 50.8 \\
			\hline
		\end{tabular}
	%

	\caption{Results on the Composition-1k testing set. GFLOPs are measured on a $224\times224\times4$ input. NL: Non-Linearity; C: Context. The lowest errors are boldfaced.}
	\label{tab:index_function}
	%

\end{table}

\begin{figure*}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=\linewidth,angle=0]{composed_1k.pdf}\vspace{-8pt}
	\caption{Qualitative results on the Composition-1k testing set. From left to right, the original image, trimap, ground-truth alpha matte, closed-form matting~\cite{levin2008closed}, deep image image~\cite{levin2008closed}, and ours (M2O DIN with ``nonlinear + context'').

	}
	\label{fig:composition-1k-results}

\end{figure*}

\begin{table*}[!t] \scriptsize
	\captionsetup{font=small,singlelinecheck=false}
	\centering
	\addtolength{\tabcolsep}{-3pt}
	\renewcommand\arraystretch{1.0}
	\begin{tabular}{l|c|ccc|ccc|ccc|ccc|ccc|ccc|ccc|ccc|ccc}
		\hline
		\multirow{2}{*}{Gradient Error} & \multicolumn{4}{c|}{Average Rank}	& \multicolumn{3}{c|}{Troll}	& \multicolumn{3}{c|}{Doll} & \multicolumn{3}{c|}{Donkey} & \multicolumn{3}{c|}{Elephant}	& \multicolumn{3}{c|}{Plant} 	& \multicolumn{3}{c|}{Pineapple} 	& \multicolumn{3}{c|}{Plastic Bag} 	& \multicolumn{3}{c}{Net} \\
		& Overall & S & L & U & S & L & U & S & L & U & S & L & U & S & L & U & S & L & U & S & L & U & S & L & U & S & L & U\\
		\hline
		IndexNet Matting & \textbf{9} & \textbf{7.3} & \textbf{7.6} & \textbf{12.3} & \textbf{0.2} & \textbf{0.2} & \textbf{0.2} & \textbf{0.1} & \textbf{0.1} & 0.3 & 0.2 & 0.2 & \textbf{0.2} & \textbf{0.2} & \textbf{0.2} & \textbf{0.4} & 1.7 & 1.9 & 2.5 & 1 & 1.1 & \textbf{1.3} & 1.1 & 1.2 & 1.2 & 0.4 & 0.5 & \textbf{0.5}\\
		\rowcolor{mygray}
		AlphaGAN~\cite{lutz2018alphagan} & 13.2 & 12 & 10.8 & 16.8 & 0.2 & 0.2 & 0.2 & 0.2 & 0.2 & 0.3 & 0.2 & 0.3 & 0.3 & 0.2 & 0.2 & 0.4 & 1.8 & 2.4 & 2.7 & 1.1 & 1.4 & 1.5 & 0.9 & 1.1 & 1 & 0.5 & 0.5 & 0.6\\
		Deep Matting~\cite{xu2017deep} & 14.3 & 10.8 & 11 & 21 & 0.4 & 0.4 & 0.5 & 0.2 & 0.2 & \textbf{0.2} & \textbf{0.1} & \textbf{0.1} & 0.2 & 0.2 & 0.2 & 0.6 & \textbf{1.3} & \textbf{1.5} & \textbf{2.4} & \textbf{0.8} & \textbf{0.9} & 1.3 & \textbf{0.7} & \textbf{0.8} & \textbf{1.1} & \textbf{0.4} & \textbf{0.5} & 0.5\\
		\hline
	\end{tabular}
	%
	\caption{Gradient errors (top 3) on the $ \tt alphamatting.com$ online benchmark.
	%
	The lowest errors are boldfaced.}
	\label{tab:alphamatting.com}
	%
\end{table*}

\begin{figure}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=3.2in,angle=0]{indices_visualization.pdf}\vspace{-8pt}
	\caption{Visualization of the randomly initialized index map (left) and the learned index map (right) of HINs (top) and DINs (bottom). Best viewed by zooming in.}
	\label{fig:indices_visualization}
	\vspace{-12pt}
\end{figure}

\vspace{3pt}
\noindent\textbf{Index Map Visualization}. It is interesting to see what indices are learned by IndexNet. For the holistic index, the index map itself is a 2D matrix and is easily to be visualized. Regarding the depthwise index, we squeeze the index map along the channel dimension and calculate the average responses. Two examples of learned index maps are visualized in Fig.~\ref{fig:indices_visualization}. We observe that, initial random indices have poor delineation for edges, while learned indices automatically capture the complex structural and textual patterns, e.g., the fur of the dog, and even air bubbles in the water.

\subsection{alphamatting.com Online Benchmark}

We also report results on the $ \tt alphamatting.com $ online benchmark~\cite{rhemann2009perceptually}. We directly test our best model trained on the Adobe Image Dataset, without fine-tuning. Our approach (IndexNet Matting) ranks the first in terms of the gradient error among published methods, as shown in Table~\ref{tab:alphamatting.com}. According to the qualitative results in Fig.~\ref{fig:alphamatting}, our approach produces significantly better mattes on hair.

\subsection{Extensions to Other Visual Tasks}

We further evaluate IndexNet on other three visual tasks. For image classification, we compare three classification networks (LeNet~\cite{lecun1998gradient}, MobileNet~\cite{howard2017mobilenets} and VGG-16~\cite{simonyan2014very}) on the CIFAR-10 and CIFAR-100 datasets~\cite{krizhevsky2009learning} with/without IndexNet. For monocular depth estimation, we attach IndexNet upon a recent ResNet-50 based baseline~\cite{hu2019revisiting} and report the performance on the NYUDv2 dataset~\cite{silberman2012indoor}. On the task of scene understanding, we evaluate SegNet~\cite{badrinarayanan2017segnet} with/without IndexNet on the SUN-RGBD dataset~\cite{song2015sun}. Results show that IndexNet consistently improves the performance in all three tasks. We refer readers to the Supplement for quantitative and qualitative results.

\begin{figure*}[!tb]
	\captionsetup{font=small,singlelinecheck=true}
	\setlength{\abovecaptionskip}{10pt}
	\centering
	\includegraphics[width=.7\textwidth,angle=0]{alphamatting.pdf
	\caption{Qualitative results on the $\tt alphamatting.com$ dataset. From left to right, the original image, deep image matting, ours.}
	\label{fig:alphamatting}
\end{figure*}

\section{Conclusion}
Inspired by an observation in image matting, we delve deep into the role of indices and present an unified perspective of upsampling operators using the notion of index function. We show that an index function can be learned within a proposed index-guided encoder-decoder framework. In this framework, indices are learned with a flexible network module termed IndexNet, and are used to guide downsampling and upsampling using two operators called $\mathcal{IP}$ and $\mathcal{IU}$. IndexNet itself is also a sub-framework that can be designed depending on the task at hand. We instantiated, investigated three index networks, compared their conceptual differences, discussed their properties, and demonstrated their effectiveness on the task of image matting, image classification, depth prediction and scene understanding. We report state-of-the-art performance on image matting with a modified MobileNetv2-based model on the Composition-1k dataset. We believe that
IndexNet is an important step towards the design of generic upsampling operators.

Our model is simple with much room for improvement. %
It may be used as a strong baseline for future research. We plan to explore the applicability of IndexNet to other dense prediction tasks.

{\bf Acknowledgments}
We would like to thank Huawei Technologies for the donation of GPU cloud computing resources.

{\small
\bibliographystyle{ieee}