text:
\section{Introduction}
\subsection{Research context}
Electro-Optical (EO) visible range sensors play a significant role in maritime target tracking~\cite{EOsurvey2017, zardoua_ais}. A recent survey~\cite{EOsurvey2017} shows that contemporary maritime target tracking algorithms require significant enhancements. Recent works indicate that reliable sea target detection algorithms should be tailored to handle unique maritime conditions~\cite{EOsurvey2017,liu2016infrared,zhang2020integrated}. The study in~\cite{EOsurvey2017} pointed out that ship detection and tracking systems are composed of three components: (1) sea horizon detection, (2) background subtraction, and (3) foreground segmentation. The development of each of these components is not a trivial problem~\cite{EOsurvey2017}; we must carefully study one component at a time to achieve significant progress. In this paper, we focus on detecting the sea horizon line.

The sea horizon line is defined as a linear shape separating the sea region and the region right above it (see Fig.~\ref{fig_horizon_examples_and_representation}). The horizon line is a valuable feature that has been incorporated in commercial systems such as ASV\textregistered (Automatic Sea Vision)~\cite{asv} and AIVS3 (Automated Intelligent Video Surveillance System for Ships)~\cite{aiv3s}. Detection of the sea horizon avoids unnecessary computations by reducing the search region of maritime targets. Moreover, the tilt $\phi$ and position $Y$ of the horizon line allow straightforward stabilization of video frames. Such stabilization enables convenient video visualization and useful temporal information of tracked targets. After localizing an obstacle on the image coordinates, the horizon line can prevent collision risks as it is involved in estimating the distance to the detected obstacles~\cite{distanceestimation1, distanceestimation2}. Unlike infrared sensors, visible range EO sensors capture, in the daytime, richer data with more color and texture changes. Despite the simple shape of the horizon, we will see that such changes make robust horizon detection harder. Robustness is not the only constraint; light and real-time computations are essential in increasing the scale and autonomy of target tracking networks.

\begin{figure}[!h]
	\centering
	\subfigure[]{
		\includegraphics[width = 0.47\linewidth, keepaspectratio]{fig_sea_horizon.png}
		\label{fig_sea_horizon}}
	\subfigure[]{
		\includegraphics[width = 0.47\linewidth, keepaspectratio]{fig_coastal_horizon.png}
		\label{fig_coastal_horizon}}

	\subfigure[]{
		\includegraphics[width = 0.47\linewidth, keepaspectratio]{fig_horizon_rep.png}
		\label{fig_horizon_rep}}
	\caption{\footnotesize{The horizon line separates the sea from the sky (a) and the sea from the coast (b); (c) position $Y$ and tilt $\phi$ of the horizon.}}
	\label{fig_horizon_examples_and_representation}
\end{figure}
\subsection{Literature review}
In this section, we succinctly discuss the horizon detection literature. Readers may refer to our previous paper for an in-depth survey~\cite{zardoua2021survey}. Gershikov et al.~\cite{gershikov(1)} detect one horizon point per column as the pixel with the maximum vertical edge response. The final horizon is fit on these points using the least-squares technique, which is known for its high sensitivity to outliers. Better algorithms project typical edge maps to another space favoring the appearance of the horizon as a peak. The first attempt in this context was conducted by Bao et al.~\cite{bao}. They detected image edges using LoG (Laplacian of Gaussian) and extracted the horizon using the global maximum of the transformed space. The horizon may not correspond to the global maximum due to image noises. Therefore, Zhang et al.~\cite{zhang-hai(1)} enhanced the work in~\cite{bao} by analyzing three geometric features of local peaks of the transformed space. Yet, such analysis assumes a high scatter of noisy edges and is thus not suitable against sea-clutter forming lined-up edges such as waves and wakes. Similar to~\cite{bao}, Schwendeman and Thomson~\cite{smtj2015} detected the horizon as the global peak of the Hough space and suggested a contrast-based quality metric to remove faulty global peaks.

Many methods applied different noise removal techniques to facilitate the eventual discrimination of the horizon. Shen et al.~\cite{shen(1), shen(2)} proposed an adaptive size of the Gaussian kernel to avoid excessive filtering. Authors of~\cite{gershikov(1), lipschutz(1)} found that morphological erosion preserves horizon edges better than Gaussian blurs. A more sophisticated morphological filter has been recently proposed by Li et al.~\cite{li2021sea}, which computes the reconstruction by erosion of the edge response. The results indicate a good performance against sea clutter appearing as blobs, but the iterative nature of such a filter slows it down. The multi-scale median filter is another way that produced the best state-of-the-art results; Prasad et al.\cite{prasadmuscowert} applied five median scales and transformed the result into a Hough space where the voting rule is modified to favor longer edges. Such a process is resource-greedy and takes tens of seconds per frame. A similar and faster alternative applies ten median scales~\cite{prasad2016mscm} and detects horizon candidates from each scale using vertical edge response as in ~\cite{gershikov(1)} and the  standard Hough transform. The final candidate is selected using a goodness score based on the strength and colinearity of edges corresponding to the considered candidate. Instead of detecting candidate lines on each median scale, Jeong et al.\cite{jeong(roi)} produced one weighted edge map computed by averaging edges of each median scale. The latter proved to be more accurate and is easily 90 times faster than~\cite{prasad2016mscm, prasadmuscowert}. Overall, multi-scale median filtering performed the best against multiple sea clutter. However, weak horizon edges usually get suppressed by the median scales.

The methods we mentioned so far rely primarily on edge information. Other works incorporate regional properties such as color and texture, which are usually fused with edge-based features to improve the accuracy and computational load; for instance, Jeong et al.~\cite{jeong(roi)} process only a region of interest extracted through analysis of color distribution difference of multiple sub-images. Other methods exploit the color properties of the sky to detect it and assume that the horizon is the linear boundary right below it~\cite{dumble, cornall(2)}. Such assumption easily breaks on images depicting coastal regions, where the horizon is the boundary separating the coast from the sea rather than the sky class. Ettinger~\cite{ettinger} and Fefilatyev et al.~\cite{fgsl} consider the horizon as the line maximizing the intra-class variance of the two regions split by that line, which requires expensive computations even on low-resolution images.

Faster methods relied on limiting the number of pixels to process; a recent algorithm proposed by Liang and Liang~\cite{liang} evaluated Hough-based candidates using color and texture information in two small patches sliding along the considered candidate. Simple assumptions about the color and texture of maritime semantic classes lead to failures even for machine-learning-based algorithms, such as SVM (Support Vector Machine), decision trees~\cite{fsgh}, and GMM (Gaussian Mixture Model)~\cite{kristan(ssm), liu2021real}. Jeong et al.~\cite{jeong(ncnn)} proposed to detect the horizon using a CNN (Convolutional Neural Network) to classify one edge point per forward pass. To reduce the computational load, they processed the input image with the multi-scale median filter proposed in~\cite{jeong(roi)} with two more scales, outputting thus a small number of candidate edge points. Unfortunately, this brings back the issue of eliminated horizon edges even when they have a decent magnitude. Zou et al.~\cite{wsl_zou2020} suggested a more robust filter by avoiding image smoothing techniques; they removed temporally non-persistent line segments and refined the result using an epipolar geometry constraint. However, the results they reported show that this method requires prohibited computations and does not suit our application.

\subsection{Main contributions and paper structure}
\label{sec_controbutions}
Previous and recent works show that the edge information is an inherent and substantial property of the horizon line; our analysis indicates that the current literature lacks an effective edge filter that satisfies two criteria: high robustness against weak horizon edges; and real-time execution through efficient and light computations. Our first contribution focuses on meeting these two constraints.

The remainder of this paper is structured as follows: Section~\ref{sec_alg} describes the proposed algorithm, Section~\ref{sec_exp} compares our method with contemporary algorithms, and Section~\ref{conclusion} concludes the work.
\section{Algorithm}
\label{sec_alg}
We will provide an overview of the proposed algorithm before detailing its stages. We show the overall algorithm pipeline in Fig.~\ref{flow}. The main ideas we employ in stages 2 and 3 (see Fig.~\ref{flow}) have been suggested after the extensive literature analysis we published in~\cite{zardoua2021survey}. Stage 1 extracts the red channel\footnote{We applied the LSD on 9 color channels on a wide range of sea images, visually analyzed the results, and concluded that the red channel produced the best results} of the input image $I_{rgb}(x,y)$ and applies a Line Segment Detector (LSD) algorithm to detect the set of line segments $S^a$. The key process is to select from $S^a$ candidate horizon segments (CHSs), which are more likely to correspond to the horizon line. We denote CHSs by set $S^f$. We rely on two filtering stages in establishing set $S^f$: the length-slope filter (LSF) and the region of interest filter (ROIF) (see Fig.~\ref{flow}).

The LSF selects from $S^a$ two sets of segments, $S^c$ and $S^d$, which contain the longest segments with a slope less than a threshold $\alpha_{th}$. We directly consider segments of set $S^c$ as CHS, i.e., $S^c \subset S^f$, because they are all longer than segments of set $S^d$. The ROIF establishes one tight region of interest (ROI) that tightly encompasses each segment of set $S^c$. Subsequently, a segment of set $S^d$ is selected as an additional CHS only if it is encompassed by at least one of the defined ROIs. We denote segments satisfying such a condition by $S^e$. We denote the set containing all CHSs by $S^f$, which is equal to $S^c \cup S^e$. Stage 4 establishes the filtered edge map $E(x,y)$ as pixels along each segment of set $S^f$. Finally, stage 5 infers the horizon using line fitting techniques and temporal information.

\begin{figure}[h!]
	\centering
	\includegraphics[width = 1\linewidth, keepaspectratio]{flowchart_v6.png}
	\caption{\footnotesize{Pipeline of our method}}
	\label{flow}
\end{figure}

\subsection{Detection of line segments}
Our method detects and evaluates image line segments instead of full lines to create a filtered edge map $E(x, y)$ with relevant horizon edges. Usage of line segments provides numerous benefits. First, we gain more robustness when the sea horizon is partially occluded with sea objects or structures. Such occlusion significantly changes horizon properties in occluded regions. As a result, even the latest CNN-based full line detectors fail at detecting the sea horizon~\cite{zardoua2021survey}. Second, line segments are inherently efficient scale-invariant descriptors and can hence be used by other algorithms of the overall sea monitoring system. For instance, exploiting the angles between line segments provides robust scale-invariant descriptors that do not only allow detecting objects but also to predict their size in cluttered backgrounds~\cite{xiao2014scale}. Third, as we will show in Sections~\ref{sec_lsf} and~\ref{sec_roif}, line segments exhibit relevant properties on cluttered maritime images. Forth, the detection of line segments is an old problem, and the current literature offers several fast and robust algorithms with sub-pixel accuracy.

In this paper, we detect line segments using the fast Line Segment Detector (LSD) developed by Von Gioi et al.~\cite{von(lsd)}, which combines the advantages of its antecedents while significantly limiting their inconveniences. Its properties can be exploited to create a robust and fast horizon edge filter: it is not misled by faulty segments induced by highly textured regions thanks to the gradient angle; a significant decrease of the false positives and negatives related to fixed detection thresholds; adopts a novel criterion that allows final selection of both long and short segments; its linear-time execution is satisfying for a large range of maritime images. As we show in Fig.~\ref{flow}, we extract the red channel and downsample it by $\kappa \in ]0, 1[$. Applying the LSD on the downsampled red channel $I'_{red}$ reduces the computational time and mitigates the negative effect of high-frequency components as well. Throughout this paper, we denote an unknown set of segments $u$ by $S^u$, its number of elements by $N^u$, and the $v$-th segment of set $S^u$ by $S^{u,v}$. Thus, set of segments $u$ can be expressed as $S^u = \{S^{u,v}\}_{v=1:N^u}$. Given this notation, we refer to segments output by the LSD by $S^a = \{S^{a,i}\}_{i=1:N^a}$. Figure~\ref{fig_filtering_steps_of_segments_SA} shows that despite the weak edge response of the horizon (see Fig.~\ref{fig_filtering_steps_of_segments_org}), set $S^a$ contains enough horizon segments to allow its detection. This positive result is explained by the fact that horizon pixels have almost equal gradient orientations, which allows the LSD~\cite{von(lsd)} to grow line segments on the horizon even when it is blurred. We favor such a possibility by setting for the LSD a small gradient magnitude threshold\footnote{To be precise, we talk about the lowest hysteresis threshold}.

\subsection{Length-Slope Filter (LSF)}
\label{sec_lsf}
\subsubsection{Description and motivation}
\label{sec_lsf_1}

The LSF (stage 2 in Fig.~\ref{flow}) filters segments of set $S^a$ using their slope and length. In our case, the EO sensor is mounted on terrestrial moving platforms whose deviation is much smaller compared to aerial-based platforms. Therefore, we suppress too tilted line segments, which may lead to detecting false sea horizons when they are remarkably long. Too tilted segments can be induced by several factors, such as the vehicle's wakes depicted by the rear camera or the long vertical edges (e.g., antennas) of nearby ships. Therefore, the LSF first filters out segments of set $S^a$ according to the slope condition of equation~\ref{eq_slope_filter}:
\begin{equation}\label{eq_slope_filter}
S^{a,i} \in S^b \text{ if } |\alpha^{a,i}| \leq \alpha_{th}
\end{equation}
where $|~.~|$ is the absolute value, $\alpha^{u, v}$ denotes the slope of segment $S^{u,v}$, $\alpha_{th}$ is a scalar threshold, and $S^b$ is the set of segments satisfying the slope condition. It is tempting to fix $\alpha_{th}$ to a very small value (e.g., $\alpha_{th} = 0.09$), but this can remove valuable horizon segments because, as reported in previous works~\cite{fgsl, kristan(ssm),EOsurvey2017}, cameras onboard small sea platforms (e.g., buoys, unmanned surface vehicles) may incur non-negligible deviations when facing sea winds and waves. We further filter segments of set $S^b$ by removing relatively short segments. Concretely, segment $S^{b, j}$ will survive the filter only if it is among the longest $N^c$ segments of set $S^b$. This condition is formally expressed in equation~\ref{eq_segs_c}:

\begin{equation}\label{eq_segs_c}
S^{b,j} \in S^c \text{ if } L^{b,j} \geq L^{b,N^c}_{srt}
\end{equation}
where $S^c$ is the set containing survived segments, $L^{u,v}$ is the length of the $v$-th segment $S^{u,v}$, $N^c$ is a scalar threshold specifying the number of the longest segments to select\footnote{$N^c$ is also the number of elements (segments) in set $S^c$.}, and $L^{u,v}_{srt}$ is the length (in pixels) of the $v$-th-longest segment of set $S^u$\footnote{For instance, $L^{b,1}_{srt}$ is the length value of the longest segment in set $S^b$}. The underlying assumption of the filtering condition of equation~\ref{eq_segs_c} is that pixels of the horizon line will have, ideally, the same gradient orientation value, which allows the LSD to group horizon pixels into the lengthiest segments. In contrast, pixels of the most common sea clutter (e.g., sea waves, sunglints, wakes) will correspond to scattered gradient orientations. In light of this, Fig.~\ref{fig_filtering_steps_of_segments_SA} shows that despite a large number of line segments on the sea surface, most of them are shorter than the few segments detected on the horizon. Hence, we perform no further filtering on segments of set $S^c$ and consider them as horizon segments, i.e., $S^c \subset S^f$ (see output of stage 3 in Fig.~\ref{flow}). Figure~\ref{fig_filtering_steps_of_segments_SC} shows segments of set $S^c$ for $N^c = 15$. We explain in Section~\ref{sec_roif} that settling for set $S^c$ as the only horizon segments could lead to missing valuable horizon segments due to some factors that affect gradient orientations of horizon pixels. Therefore, we select more segments by taking the next longest $N^d$ segments from set $S^b$ as expressed in equation~\ref{eq_segs_d}:

\begin{equation}\label{eq_segs_d}
S^{b,j} \in S^d \text{ if }  L^{b,N^{c}}_{srt} > L^{b,j} \geq L^{b,N^{d'}}_{srt}
\end{equation}
where $N^{d'} = N^c + N^d$ and $N^d$ is a scalar threshold that correspond to the number of segments in set $S^d$ ($N^d \gg N^c$; e.g., $N^d = N^c \times 10$). Unlike set $S^c$, we do not directly consider segments in set $S^d$ as horizon segments because their length is relatively close to that of noisy segments, as shown in Fig.~\ref{fig_filtering_steps_of_segments_SD}. The next stage captures additional horizon segments from set $S^d$ using information extracted from set $S^c$.
\begin{figure}[!h]
	\centering
	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_org.png}
		\label{fig_filtering_steps_of_segments_org}}
	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_SA.png}
		\label{fig_filtering_steps_of_segments_SA}}
	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_SC.png}
		\label{fig_filtering_steps_of_segments_SC}}

	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_SD.png}
		\label{fig_filtering_steps_of_segments_SD}}
	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_SE.png}
		\label{fig_filtering_steps_of_segments_SE}}
	\subfigure[]{
		\includegraphics[width = \figwidthb, keepaspectratio]{fig_filtering_steps_of_segments_SF.png}
		\label{fig_filtering_steps_of_segments_SF}}

	\caption{\footnotesize{Intermediate results of line segments filtering: (a) original image; (b) segments of set $S^a$; (c) segments of set $S^c$; (d) segments of set $S^d$; (e) segments of set $S^e$; (f) segments of set $S^f$}}
	\label{fig_filtering_steps_of_segments}
\end{figure}
\subsubsection{LSF: implementation details}
\label{sec_lsf_2}

We denote Cartesian coordinates of the starting point and ending point\footnote{We arbitrarily consider that the starting point is always on the left of the ending point.} of the $v$-th segment $S^{u,v}$ of set $S^u$ by $p^{u,v}_{s}:~(x^{u,v}_{s}, y^{u,v}_{s})$ and $p^{u,v}_{e}:~(x^{u,v}_{e}, y^{u,v}_{e})$, respectively. Thus, we denote $x^u_s$, $x^u_e$, $y^u_s$, and $y^u_e$ as vectors containing coordinates of segments of set $S^u$:

\begin{equation}\label{eq_xy_segs}
\begin{split}
x^u_s&= [x^{u,1}_s, x^{u,2}_s, \dots x^{u,N^u}_s]^T \in \mathbb{R}^{N^u \times 1}\\
x^u_e&= [x^{u,1}_e, x^{u,2}_e, \dots x^{u,N^u}_e]^T \in \mathbb{R}^{N^u \times 1}\\
y^u_s&= [y^{u,1}_s, y^{u,2}_s, \dots y^{u,N^u}_s]^T \in \mathbb{R}^{N^u \times 1}\\
y^u_e&= [y^{u,1}_e, y^{u,2}_e, \dots y^{u,N^u}_e]^T \in \mathbb{R}^{N^u \times 1}
\end{split}
\end{equation}

There are two processes to vectorize in the LSF: the slope-based (equation~\ref{eq_slope_filter}) and length-based (equations~\ref{eq_segs_c} and~\ref{eq_segs_d}) filtering conditions. The former starts by computing vector $\alpha^a$ according to equation~\ref{eq_alpha_a}:
\begin{equation}\label{eq_alpha_a}
\begin{split}
\alpha^a & = [\alpha^{a,1}, \alpha^{a,2},\dots, \alpha^{a,N^a}]^T \in \mathbb{R}^{N^a\times 1}\\
& = (y^a_e - y^a_s)\oslash(x^a_e - x^a_s)
\end{split}
\end{equation}
where $\alpha^u$ is the vector corresponding to the slopes of segments in $S^u$ and $\oslash$ is the Hadamard (component-wise) division. All filtered sets of segments in this paper, i.e., $S^b$, $S^c$, $S^d$, $S^e$, and $S^f$, represent subsets of set $S^a$ and are thus established by a vectorized indexing; concretely and generally, we use the different filtering conditions (e.g., equations~\ref{eq_slope_filter} and~\ref{eq_segs_c}) to establish vector $I_{u',u} \in \mathbb{R}^{N^{u'}\times 1}$ whose values are indices that would obtain the filtered set $S^{u'}$ by indexing set $S^u$. We will see in Section~\ref{sec_roif_3} that vector $I_{u',u}$ also allows obtaining attributes of a given filtered set $S^{u'}$ from homologous attributes of one of its supersets $S^{u}$\footnote{For instance, we can get the slope vector $\alpha^{u'}$ by indexing $\alpha^{u}$ with vector $I_{u',u}$.}, avoiding thus redundant computations. Using discussed indexing, we get segments of set $S^b$ by establishing vector $I_{b,a} \in \mathbb{R}^{N^b \times 1}$ whose values are indices of slopes $\alpha^{a,i}$ (see equation~\ref{eq_alpha_a}) satisfying the condition expressed in equation~\ref{eq_slope_filter}. Eventually, we index set $S^a$ by vector $I_{b,a}$ to produce the desired set $S^b$. As set $S^c$ and $S^d$ are established based on the length of segments of set $S^b$, we first compute a length vector according to equation~\ref{eq_Len_b}:

\begin{equation}\label{eq_Len_b}
\begin{split}
L^{b} &= [L^{b,1}, L^{b,2},\dots, L^{b,N^b}]^T \in \mathbb{R}^{N^b \times 1} \\
&= \left((x^b_e - x^b_s)^{\odot 2}-(y^b_e - y^b_s)^{\odot 2}\right) ^{\odot \nicefrac{1}{2}}
\end{split}
\end{equation}
where $L^{u, v}$ is the length of the $v$-th segment $S^{u, v}$ and $^{\odot}$ is the Hadamard (element-wise) power. Then, we establish vector $I^{L^b}_{srt} \in \mathbb{R}^{N^b\times 1}$ whose values are indices that would sort values of $L^b$ from highest to lowest. Thus, vector $I_{c,b}$, whose elements are indices of set $S^b$ satisfying equation~\ref{eq_segs_c}, is established by slicing the first $N^c$ elements of vector $I^{L^b}_{srt} \in \mathbb{R}^{N^b\times 1}$. Similarly, we get vector $I_{d,b}$, whose elements are indices of set $S^b$ satisfying equation~\ref{eq_segs_d}, by slicing $I^{L^b}_{srt}$ from the $(N^c + 1)$-th element to the $(N^c + 1 + N^d)$-th element. Eventually, we get set $S^c$ and $S^d$ by indexing set $S^b$ using vectors $I_{c,b}$ and $I_{d,b}$, respectively.

\subsection{Region Of Interest Filter (ROIF)}
\label{sec_roif}
\subsubsection{Motivation and assumption}
\label{sec_roi_0}
We mentioned in Section~\ref{sec_lsf} that taking the longest $N^c$ segments from set $S^b$ could miss valuable line segments. This issue is fundamentally caused by certain types of noise, such as the Gaussian noise of low quality sensors and the poor horizon contrast. Such noise disturb the gradient orientation of horizon pixels, hindering thus the LSD to grow long segments on the horizon line. Fig.~\ref{fig_filtering_steps_of_segments_SD} shows that the low contrast condition led to growing short segments on the horizon. The same disturbed property, i.e., gradient orientation, provides the key to mitigating this issue; we assume that, unlike noisy segments, short horizon segments have a much higher colinearity with longer horizon segments. Figure~\ref{fig_synthetic_assumption_roi_1} depicts a synthetic example of this assumption. We show the result of using this assumption on the real image in Fig.~\ref{fig_filtering_steps_of_segments_org}; Fig.~\ref{fig_filtering_steps_of_segments_SE} shows that additional horizon segments, denoted as set $S^e$, are successfully captured from $S^d$. Selecting more horizon segments from $S^d$ often leads to capturing more noisy segments as well, as shown in Fig.~\ref{fig_filtering_steps_of_segments_SE}. However, our experiments indicate that having more horizon segments allows more performance.

\begin{figure}[!h]
	\centering
	\subfigure[]{
		\includegraphics[width = \figwidtha, keepaspectratio]{fig_synthetic1.png}
		\label{fig_synthetic_assumption_roi_1}}
	\subfigure[]{
		\includegraphics[width = \figwidtha, keepaspectratio]{fig_synthetic2.png}
		\label{fig_synthetic_defined_rois}}
	\subfigure[]{
		\includegraphics[width = \figwidtha, keepaspectratio]{fig_synthetic3.png}
		\label{fig_synthetic_assumption_survived_segs}}
	\caption{\footnotesize{A synthetic example showing the assumption made in the region of interest filter and its usage; the green dotted line is the horizon line; (a) segments of set $S^c$ (in red) and $S^d$ (in blue); (b) the filtering regions $\mathcal{R}^{c,1}$ and $\mathcal{R}^{c,2}$ (in pink) encompassing the two segments of set $S^c$; (c) segments of set $S^d$ that survived the two filtering regions (in yellow)}}
	\label{fig_synthetic_assumption_roi}
\end{figure}

\subsubsection{The filtering condition: initial formulation}
\label{sec_roif_1}
We aim in this Section at formulating an initial expression of the filtering condition based on the colinearity assumption we mentioned in Section~\ref{sec_roi_0} and illustrated in Fig.~\ref{fig_synthetic_assumption_roi_1}. To this end, we define one region of interest, denoted as $\mathcal{R}^{c,k}$, to encompass the $k$-th segment $S^{c,k}$. Subsequently, we consider that the $l$-th segment $S^{d,l}$ is an additional horizon segment only if both of its endpoints fall in at least one and the same region $\mathcal{R}^{c,k}$. Figure~\ref{fig_synthetic_defined_rois} shows in pink texture an example of the filtering regions that would correspond to Fig.~\ref{fig_synthetic_assumption_roi_1}. Figure~\ref{fig_synthetic_assumption_survived_segs} shows in yellow the segments of set $S^d$ satisfying the condition we just mentioned. We initially formulate this condition in equation~\ref{eq_roi_cond_1}:

\begin{equation}\label{eq_roi_cond_1}
S^{d,l} \in S^e\text{ if }\exists~k~\text{:}~\left((p^{d,l}_{s} \in \mathcal{R}^{c,k})\land~(p^{d,l}_{e} \in \mathcal{R}^{c,k})\right) = \mathsf{1}
\end{equation}
where $S^e$ is the set grouping additional horizon segments selected from set $S^d$, $p^{d,l}_{s}$ and $p^{d,l}_{e}$ are the starting and ending point of segment $S^{d,l}$, respectively, $\in$ is the inclusion operator , $\land$ is the logical \textit{and}, and $\mathsf{1}$ is the true Boolean.

\subsubsection{Mathematical development of equation~\ref{eq_roi_cond_1}}
\label{sec_roif_2}
Implementing the filtering condition in equation~\ref{eq_roi_cond_1} first requires further mathematical development. This breaks down to verifying that a given endpoint of segment $S^{d,l}$ is included in a given region $\mathcal{R}^{c,k}$. We execute this task in four steps: (1) set the width $W_{\mathcal{R}}$ of all regions $\mathcal{R}^{c,k}$ to an arbitrary small value; (2) define function $g^{c,k}(x)$ as the linear curve crossing both endpoints $p_s^{c,k}$ and $p_e^{c,k}$ of segment $S^{c,k}$; (3) compute $\forall h \in \{s, e\}$ the quantity $\mathbbm{d}^{k,l}_{h}$, which is the normal distance from point $p^{d,l}_{h}$ to line $g^{c,k}(x)$; (4) consider that $p^{d,l}_{h} \in \mathcal{R}^{c,k}$ only if $\mathbbm{d}^{k,l}_{h} < t_{roi}$, where $t_{roi}$ is the parameter controlling the width $W_{\mathcal{R}}$: $t_{roi} = 2 \times W_{\mathcal{R}}$. We eventually express the developed filtering condition equivalent to equation~\ref{eq_roi_cond_1} in equation~\ref{eq_roi_cond_2}, where quantity $\mathbbm{d}^{k,l}_{h}$ is computed using equations~\ref{eq_Dy} and~\ref{eq_beta}:

\begin{equation}\label{eq_roi_cond_2}
S^{d,l} \in S^e\text{ if }\exists~k~\text{:}~\left((\mathbbm{d}^{k,l}_{s} \leq t_{roi})\land~(\mathbbm{d}^{k,l}_{e} \leq t_{roi})\right) = \mathsf{1}
\end{equation}

\begin{equation}\label{eq_Dy}
\begin{split}
\mathbbm{d}^{k,l}_{h} &= \frac{1}{\sqrt{1 + (\alpha^{c,k})^2}}\times|g^{c,k}(x^{d,l}_h) - y^{d,l}_h|;~\forall h \in \{s, e\}\\
&= \frac{1}{\sqrt{1 + (\alpha^{c,k})^2}}\times|\alpha^{c,k} \times x^{d,l}_h+ \beta^{c,k} - y^{d,l}_h|
\end{split}
\end{equation}

\begin{equation}\label{eq_beta}
\begin{split}
\beta^{c,k} &= y^{c,k}_s - \alpha^{c,k} \times x^{c,k}_s\\
&= y^{c,k}_e - \alpha^{c,k} \times x^{c,k}_e
\end{split}
\end{equation}
where $\alpha^{c,k}$ and $\beta^{c,k}$ correspond to the slope and intercept of curve $g^{c,k}(x)$, respectively. Note that $\alpha^{c,k}$ is equal to the slope of the $k$-th segment $S^{c,k}$ and can be obtained from vector $\alpha^a$.

\subsubsection{Implementation details}
\label{sec_roif_3}
We present in this Section the details to vectorize the filtering condition we developed in equation~\ref{eq_roi_cond_2}. We define matrix $\mathbbm{D}_h$ to contain $\forall (k,l)$ the quantities $\mathbbm{d}^{k,l}_{h}$. We illustrate $\mathbbm{D}_h$ in equation~\ref{eq_DY_h_elements} and compute it using equation~\ref{eq_DY_h_formula}:
\begin{equation}\label{eq_DY_h_elements}
\mathbbm{D}_h = \begin{bmatrix}
\mathbbm{d}^{1,1}_{h}     & \mathbbm{d}^{1,2}_{h}   & \dots  & \mathbbm{d}^{1,N^d}_{h}   \\
\mathbbm{d}^{2,1}_{h}     & \mathbbm{d}^{2,2}_{h}   & \dots  & \mathbbm{d}^{2,N^d}_{h}   \\
\vdots & \vdots & \ddots  & \vdots\\
\mathbbm{d}^{N^c,1}_{h}   & \mathbbm{d}^{N^c,2}_{h} & \dots  & \mathbbm{d}^{N^c,N^d}_{h}
\end{bmatrix}
= (\mathbbm{d}^{k,l})\in \mathbb{R}^{N^c \times N^d}
\end{equation}

\begin{equation}\label{eq_DY_h_formula}
\mathbbm{D}_h = |\alpha^c \times (x^d_h)^T + B^c - Y^d_h|\oslash ((1^c + (\alpha^c)^{\odot 2})^{\odot \nicefrac{1}{2}} \times (1^d)^T)
\end{equation}
where $(.)^T$ denotes vector transpose, $B^c\in \mathbb{R}^{N^c \times N^d}$ contains intercept values of all functions of set $g^c = \{g^{c,k}\}_{k = 1:N^c}$, $Y^d_h \in \mathbb{R}^{N^c \times N^d}$ contains $y$ Cartesian coordinates of all endpoints in set $p^{d}_h = \{p^{d,l}_h\}_{l = 1:N^d}$, and $1^u = [1, 1,\dots, 1]^T \in \mathbb{R}^{N^u \times 1}$ is an all-ones vector. We provide in what follows the details for getting $\alpha^c$, $B^c$, and $Y^d_h$. We establish $\alpha^c$ by indexing vector $\alpha^a$ using vector $I_{b,a}$, which produces vector $\alpha^b$. We index the latter using vector $I_{c,b}$, producing thus desired slopes $\alpha^c$. We already established both vectors $I_{b,a}$ and $I_{c,b}$ in Section~\ref{sec_lsf_2}. We show elements of $B^c$ in equation~\ref{eq_B_c_elements}, and compute them by calculating vector $\beta^c$ according to equation~\ref{eq_beta_c} and broadcasting the result to $N^d$ columns, as shown in equation~\ref{eq_B_c}. We show elements of matrix $Y^d_h$ in equation~\ref{eq_Y_d_h_elements} and compute them using equation~\ref{eq_Y_d_h_formula}.
\begin{equation}\label{eq_B_c_elements}
B^c =
\begin{bmatrix}
\beta^{c,1}   & \beta^{c,1}   & \dots  & \beta^{c,1}\\
\beta^{c,2}   & \beta^{c,2}   & \dots  & \beta^{c,2}\\
\vdots        & \vdots        & \ddots & \vdots     \\
\beta^{c,N^c} & \beta^{c,N^c} & \dots  & \beta^{c,N^c}
\end{bmatrix}
\in \mathbb{R}^{N^c \times N^d}
\end{equation}

\begin{equation}\label{eq_beta_c}
\begin{split}
\beta^c =& [\beta^{c, 1}, \beta^{c, 2},\dots, \beta^{c, N^c}]^T \in \mathbb{R}^{N^c \times 1}\\
=& y^c_h - (\alpha^c \odot x^c_h)
\end{split}
\end{equation}

\begin{equation}\label{eq_B_c}
B^c = \beta^c \times (1^d)^T
\end{equation}

\begin{equation}\label{eq_Y_d_h_elements}
Y^d_h = \begin{bmatrix}
y^{d,1}_h   & y^{d,2}_h   & \dots  & y^{d,N^d}_h\\
y^{d,1}_h   & y^{d,2}_h   & \dots  & y^{d,N^d}_h\\
\vdots        & \vdots        & \ddots & \vdots \\
y^{d,1}_h   & y^{d,2}_h   & \dots  & y^{d,N^d}_h
\end{bmatrix}
\in \mathbb{R}^{N^c \times N^d}
\end{equation}

\begin{equation}\label{eq_Y_d_h_formula}
Y^d_h = 1^c \times (y^d_h)^T \in \mathbb{R}^{N^c \times N^d}
\end{equation}

To exploit the distances in matrix $\mathbbm{D}_h$, we define $Q_h$ as in equation~\ref{eq_Q_h_elements} and compute it using equation~\ref{eq_Q_h_formula}:

\begin{equation}\label{eq_Q_h_elements}
Q_h =
\begin{bmatrix}
q^{1,1}_h     & q^{1,2}_h     & \dots  & q^{1,N^d}_h\\
q^{2,l}_h     & q^{2,2}_h     & \dots  & q^{2,N^d}_h\\
\vdots        & \vdots        & \ddots & \vdots \\
q^{N^c,l}_h   & q^{N^c,2}_h   & \dots  & q^{N^c,N^d}_h
\end{bmatrix}
= (q^{k,l}_h) \in \mathbb{R}^{N^c \times N^d}
\end{equation}

\begin{equation}\label{eq_Q_h_formula}
Q_h = \mathbbm{D}_h \leq T_{roi} \in \mathbb{R}^{N^c \times N^d};~\forall h \in ~\{s, e\}
\end{equation}
where $q^{k,l}_h$ is a Boolean scalar whose truth indicates that endpoint $p_h^{d,l}\in \mathcal{R}^{c,k}$ and $T_{roi} \in \mathbb{R}^{N^c \times N^d}$ is a matrix whose all elements are equal to $t_{roi}$ (see equation~\ref{eq_roi_cond_2}). The comparison performed in equation~\ref{eq_Q_h_formula} is equivalent to that of equation~\ref{eq_roi_cond_2}. We further process matrix $Q_h$ using equation~\ref{eq_q_formula}:
\begin{equation}\label{eq_q_formula}
\begin{split}
q^d &= \overset{\downarrow}{\lor}(Q_s \tilde{\land} Q_e) \in \mathbb{R}^{1 \times N^d}\\
&= [q^{d,1}, q^{d,2}, \dots q^{d,N^d}]
\end{split}
\end{equation}
where $\overset{\downarrow}{\lor}(.)$ operator computes the logical \textit{or} along the vertical axis of matrix $Q_h$\footnote{In other words, $\overset{\downarrow}{\lor}$ performs the logical \textit{or} operation on elements of each column}, $\tilde{\land}$ computes the element-wise logical \textit{and}, and $q^{d,l}$ is a logical Boolean whose truth indicates that segment $S^{d,l}$ is encompassed by at least one region of interest $\mathcal{R}^{c,k}$. We note that the \textit{and} operator used in equation~\ref{eq_q_formula} reflects the \textit{and} operator in equations~\ref{eq_roi_cond_1} and~\ref{eq_roi_cond_2}, whereas the \textit{vertical or} operator in equation~\ref{eq_q_formula} reflects the $\exists$ operator in equations~\ref{eq_roi_cond_1} and~\ref{eq_roi_cond_2}. Thus, equation~\ref{eq_roi_cond_2}, which is equivalent to equation~\ref{eq_roi_cond_1}, becomes equivalent to equation~\ref{eq_roi_cond_3}. Subsequently, we create vector $I_{e,d}$ to contain indices where elements of vector $q^d$ (see equation~\ref{eq_q_formula}) are true Booleans. Finally, we establish set $S^e$ by indexing $S^d$ using $I_{e,d}$. The final filtered set of segments output by the ROIF (stage 3 in Fig.~\ref{flow}) is computed as in equation~\ref{eq_S^f}. Fig.~\ref{fig_filtering_steps_of_segments_SA} and~\ref{fig_filtering_steps_of_segments_SF} show an example of the original set of segments $S^a$ and the corresponding filtered set $S^f$, respectively.

\begin{equation}\label{eq_roi_cond_3}
S^{d,l} \in S^e\text{ if } q^{d,l} = \mathsf{1}
\end{equation}

\begin{equation}\label{eq_S^f}
S^f = S^c \cup S^e
\end{equation}

\subsection{Segments to Edge Pixels (STEP)}
\label{sec_step}
To get filtered edge points, we must find coordinates of all pixels along each segment in set $S^f$. The number of pixels to locate for each segment $S^{f,n}$ is equal to its length $L^{f,n}$. Thus, we represent coordinates of pixels along segment $S^{f,n}$ as in equation~\ref{eq_z^n_element} and compute them using equation~\ref{eq_z^n_formula}:

\begin{equation}\label{eq_z^n_element}
\begin{split}
x^n &= [x^{n,1}, x^{n,2},\dots x^{n,L^{f,n}}]^T \in \mathbb{R}^{L^{f,n}\times 1}\\
y^n &= [y^{n,1}, y^{n,2},\dots y^{n,L^{f,n}}]^T \in \mathbb{R}^{L^{f,n}\times 1}
\end{split}
\end{equation}

\begin{equation}\label{eq_z^n_formula}
\begin{split}
x^n &= \frac{x_e^{f,n}-x_s^{f,n}}{L^{f,n} - 1}\odot z^n + x_s^{f,n}\\
y^n &= \frac{y_e^{f,n}-y_s^{f,n}}{L^{f,n} - 1}\odot z^n + y_s^{f,n}
\end{split}
\end{equation}
where $z^n = [0, 1, 2, \dots, L^{f, n} - 1]^T \in \mathbb{R}^{L^{f,n}\times 1}$. Because detected segments do not have the same length, we cannot further vectorize equation~\ref{eq_z^n_formula} as we did with previous equations. Therefore, we iterate over all segments of $S^f$. At each iteration, we compute equation~\ref{eq_z^n_formula} and append the result as expressed in equation~\ref{eq_xy^out}:
\begin{equation}\label{eq_xy^out}
\begin{split}
x^{out} &:= x^{out~\frown}~x^n;~\forall n \in \{1, 2, \dots, N^f - 1, N^f\}\\
y^{out} &:= y^{out~\frown}~y^n;~\forall n \in \{1, 2, \dots, N^f - 1, N^f\}
\end{split}
\end{equation}
where $:=$ represents value assignment, $^{\frown}$ denotes vector appending, and $x^{out}$ and $y^{out}$ are vectors that would contain Cartesian coordinates of pixels along all segments in set $S^f$. Both $x^{out}$ and $y^{out}$ are initialized to an empty vector at the first iteration. Although this process is not fully vectorized, the effect on real-time performance is not significant as the portion of survived segments in $S^f$ is tiny compared to other sets of segments. We will see in Section~\ref{sec_experimental_results} that the real-time performance achieved is satisfying.

Next, we use vectors $x^{out}$ and $y^{out}$ to establish an edge map image $E^{'}(x,y)$ (see stage 4 in Fig.~\ref{flow}) corresponding to the downsampled image $I'_{red}$. Fitting the horizon line on $E^{'}$ works well but establishing horizon parameters corresponding to the original frame size requires scaling up the horizon position $Y^{'}$ corresponding to $E^{'}$: $Y = Y^{'}\times\frac{1}{\kappa}$. This affects the detected line because we amplify the error corresponding to $Y^{'}$ by $\frac{1}{\kappa}$. Scaling up parameters of the horizon is necessary not only for comparison with Ground Truth (GT) position $Y^{GT}$ but for subsequent applications as well. For instance, the computation of transformation matrices involved in stabilizing original video frames is directly related to horizon parameters $\{Y, \phi\}$~\cite{cai, fgsl, smtj2015}. Therefore, we infer the horizon line by relocating edge points of $E'$ on an edge map $E$ (see stage 4 in Fig.~\ref{flow}) with the original frame size. Thus, we significantly mitigate the amplified error and leverage image downsizing. We compute $E$ by upsampling $E^{'}$ to the original frame size using a bilinear interpolation method, which outputs grayscale image $E^{''}$ with thicker edges. We shrink the thickness of the latter to one pixel using the nonmaximum suppression method, followed by single-thresholding of $E^{''}$ with $E_{th} = 254$. This outputs image $E$, which contains relocated edge points. Figure~\ref{fig_edge_output} shows the edge map $E$ corresponding to challenging conditions, such as weak horizon edges, coastal boundaries, and highly textured regions induced by sea waves, glints, and clouds; the large number of segments $S^a$ in Fig.~\ref{fig_edge_output_segments} indicates the high amount of noise; Fig.~\ref{fig_edge_output_edges} demonstrates that the high image clutter is significantly suppressed while favoring horizon edges to form the most prominent line. Edges of the horizon line corresponding to the right image of Fig.~\ref{fig_edge_output_org} are dramatically affected by the smooth color transition between the sea and sky. Yet, the corresponding edge map in Fig.~\ref{fig_edge_output_edges} demonstrates that the filtering stages managed to keep enough horizon edges to allow an accurate detection. This will be further demonstrated qualitatively and quantitatively in Section~\ref{sec_experimental_results}.
\begin{figure*}[!h]
	\centering
	\subfigure[]{
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_1_org.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_2_org.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_3_org.png}
		\label{fig_edge_output_org}}

	\subfigure[]{
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_1_segs.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_2_segs.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_3_segs.png}
		\label{fig_edge_output_segments}}

	\subfigure[]{
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_1_edges.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_2_edges.png}
		\includegraphics[width = \figwidthedge, keepaspectratio]{fig_filtered_3_edges.png}
		\label{fig_edge_output_edges}}
	\caption{\footnotesize{Illustration of the output edge map $E$; (a) original images; (b) all detected line segments; (c) the output edge map $E$. For visualization convenience, the binary edges in (c) are dilated and inverted, and the horizon edges are colored in green.}}
	\label{fig_edge_output}
\end{figure*}

\subsection{Horizon line inference}
\label{sec_horizon_inference}
We infer the horizon line in three major steps (see stage 5 in Fig.~\ref{flow}). The Outlier Handler Module (OHM) verifies that the global maximum of the accumulated Hough space, denoted as $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\}$, is not an outlier. We define outlier lines as false-positive horizons persisting over a very brief period. This issue occurs on other algorithms we tested in this paper as well. In our case, outlier lines appear primarily due to noisy segments of the sea region that occasionally line up in a way to survive the ROIF. The OHM considers $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\}$ as an outlier line if the condition in equation~\ref{eq_outlier_condition} is satisfied:
\begin{equation}\label{eq_outlier_condition}
(|Y^{\mathcal{H}} - Y^{prv}| > \Delta Y^{th}) \lor (|\phi^{\mathcal{H}} - \phi^{prv}| > \Delta \phi^{th}) = \mathsf{1}
\end{equation}
where $\lor$ is the logical \textit{or}, $\{Y^{prv}, \phi^{prv}\}$ are parameters of the most recent horizon line, and $\{\Delta Y^{th}, \Delta \phi^{th}\}$ are scalar thresholds. If $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\}$ is not an outlier, we consider it as a coarse estimation of the horizon line: $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\} = \{Y^{crs}, \phi^{crs}\}$. Otherwise, if $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\}$ is an outlier, the OHM will consider that the coarse horizon line $\{Y^{crs}, \phi^{crs}\}$ is one of the longest $M$ Hough lines satisfying equation~\ref{eq_outlier_condition}. If multiple lines satisfy such a condition, we select the line corresponding to the minimum value $|Y^{\mathcal{H}} - Y^{prv}|$. To refine the coarse line $\{Y^{crs}, \phi^{crs}\}$, we follow Jeong et al.~\cite{jeong(roi)} by applying the least-squares fitting algorithm on \textit{inlier edge points}, i.e., edge points that voted on $\{Y^{crs}, \phi^{crs}\}$.

Ettinger et al.~\cite{ettinger} pointed out that using previous detections of the horizon line to infer the horizon may get the algorithm into a failure state; concretely, if the algorithm detects a faulty line, it will incorrectly update the parameters of the most recent horizon line $\{Y^{prv}, \phi^{prv}\}$. Therefore, the outlier condition in equation~\ref{eq_outlier_condition} becomes useless if the difference between that faulty line and the true horizon exceeds a certain threshold. When this issue occurs, the faulty horizon is unlikely to be persistently detected on subsequent frames, either due to the changing sea noise or the linear property of the horizon that would compete with that faulty line. Thus, the outlier condition triggers on a significant number of consecutive frames. The OHM counts this number, denoted as $N_{outs}$, and compares it to a threshold $N^{th}_{outs}$. If $N_{outs} > N^{th}_{outs}$, the OHM considers that the algorithm is in a failure state. To get the algorithm out of this state, the OHM avoids finding a substitute line and directly refines the longest Hough line $\{Y^{\mathcal{H}}, \phi^{\mathcal{H}}\}$ using the least-squares fitting. This process quickly converges to correctly updating the parameters $\{Y^{prv}, \phi^{prv}\}$.

\section{Experiments}
\label{sec_exp}
\subsection{Implementation of benchmarked methods}
\label{sec_impl_of_bench}
We compare our algorithm with the best state-of-the-art algorithms~\cite{li2021sea, liang, jeong(roi), gershikov(1)}. The only source code we could obtain is a Matlab implementation of Liang and Liang's method~\cite{liang}. We implemented all other algorithms~\cite{li2021sea, jeong(roi), gershikov(1)}, including ours, using Python 3.8. The main packages we used are \textit{Numpy} and \textit{OpenCV}. We run all algorithms, one at a time, on a computer with 8 GB of RAM and an Intel\textregistered~Core\texttrademark~i5-6300U CPU @ 2.40GHz. For a fair comparison, we made sure to stop all tasks that may use the computer's CPU during the execution of each algorithm.

\subsection{Evaluation protocol}
\label{sec_eval_protocol}
We parameterize the detected horizon line using the position $Y$ and tilt $\phi$ we previously exposed in Fig.~\ref{fig_horizon_rep}. We quantify the detection error by computing $Y^{\epsilon} = |Y - Y^{GT}|$ and $\phi^{\epsilon} = |\phi - \phi^{GT}|$, where the pair $\{Y^{GT}, \phi^{GT}\}$ corresponds to parameters of the GT (annotated) horizon. Eventually, we use errors computed to extract the following six statistical metrics: $\mu$, $\sigma$, Q25, Q50, Q75, and Q95, where $\mu$ and $\sigma$ correspond to the mean and standard deviation, respectively, and Q$P$ is the $P$-th percentile. To quantify the real-time performance, we measure the mean processing time of all algorithms.

\subsection{Experimental results and discussion}
\label{sec_experimental_results}
Fig.~\ref{fig_visual_results} shows that our method is highly effective against edge-degraded horizons (see Figs.~\ref{fig_visual_results_1}, ~\ref{fig_visual_results_2}, and ~\ref{fig_visual_results_3}) while handling other maritime challenges as well, such as glints, waves, shore edges, ship occlusion, and textured clouds. Tables~\ref{tab_result_smd_onboard} and~\ref{tab_result_smd_onshore} show that our algorithm and Jeong et al.'s method~\cite{jeong(roi)} perform significantly better on the SMD. Table~\ref{tab_result_smd_onboard} shows that our method performed the best in terms of 9 metrics (out of 12), while Jeong et al.'s method~\cite{jeong(roi)} scored the best on the remaining three metrics. Such a good performance of~\cite{jeong(roi)} is expected because the employed multi-scale median filter is highly effective against a variety of sea clutter when the horizon is long with strong edges, which is the case in most SMD-onboard videos. Table~\ref{tab_result_smd_onshore} shows that the performance of all algorithms decreased. We explain this by the fact that the horizon in SMD-onshore videos is relatively shorter due to ship occlusion and mixed with additional types of clutter such as linear wakes and contours of ships. Nevertheless, Table~\ref{tab_result_smd_onshore} shows that our method performs remarkably better on 10 (out of 12) metrics while scoring satisfying values on the remaining two metrics. Table~\ref{tab_execution_time} shows the execution time of all methods on each dataset. Our method is not the fastest, but its speed-accuracy trade-off makes it suitable for real-time target tracking.
\begin{figure}[!b]
	\centering

	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_1_org.png}
		\label{fig_visual_results_org_1}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_2_org.png}
		\label{fig_visual_results_org_2}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_3_org.png}
		\label{fig_visual_results_org_3}}

	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_4_org.png}
		\label{fig_visual_results_org_4}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_5_org.png}
		\label{fig_visual_results_org_5}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_6_org.png}
		\label{fig_visual_results_org_6}}

	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_1_zardoua.png}
		\label{fig_visual_results_1}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_2_zardoua.png}
		\label{fig_visual_results_2}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_3_zardoua.png}
		\label{fig_visual_results_3}}

	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_4_zardoua.png}
		\label{fig_visual_results_4}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_5_zardoua.png}
		\label{fig_visual_results_5}}
	\subfigure[]{
		\includegraphics[width = \figwidthc, keepaspectratio]{fig_vis_comp_6_zardoua.png}
		\label{fig_visual_results_6}}

	\caption{\footnotesize{Visual results of the horizon detected by our method: (a) to (f) are the original images; (g) to (l) are the same images drawn with the detected horizon line}}
	\label{fig_visual_results}
\end{figure}

\begin{table*}[!h]
	\centering
	\caption{\label{tab_result_smd_onboard}\footnotesize{Quantitative results on \textit{onboard} videos of the Singapore Maritime Dataset}}\vspace{-2mm}
	\setlength{\extrarowheight}{5pt}
	{\footnotesize
\begin{tabular}{P{0.11\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.01\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}
		        P{0.045\linewidth}}
	        \cline{2-13}
			&\multicolumn{6}{c}{$Y^{\epsilon} = |Y - Y^{GT}|$ in pixels}&\multicolumn{6}{c}{$\phi^{\epsilon} = |\phi - \phi^{GT}|$ in degrees}\\
			\cline{2-7}
			\cline{9-13}
			&$\mu$&$\sigma$&Q25&Q50&Q75&Q95&~~~&$\mu$&$\sigma$&Q25&Q50&Q75&Q95\\
			\hline

			Ours
			&1.95&1.59&0.84&1.71&2.71&4.57&~~~&0.16&0.13&0.060&0.13&0.23&0.41

			\\
			\hline

			Li~\cite{li2021sea}
			&34.64&73.96&14.35&22.72&30.84&71.49&~~~&0.85&1.46&0.37&0.64&1&1.70

			\\
			\hline

			Liang~\cite{liang}
			&17.77&46.93&0.90&2.13&7.65&93.43&~~~&4.71&3.38&2.53&4.27&6.28&10.09

			\\
			\hline

			Jeong~\cite{jeong(roi)}
			&3.84&33.09&0.68&1.43&2.44&4.86&~~~&0.26&1.61&0.065&0.14&0.25&0.50

			\\
			\hline

			Gershikov~\cite{gershikov(1)}
			&301.44&145.02&186.67&334.50&410.28&510.67&~~~&6.55&6.85&1.60&3.84&9.17&22.10

			\\
			\hline
		\end{tabular}
	}
\end{table*}

\begin{table*}[!h]
	\centering
	\caption{\label{tab_result_smd_onshore}\footnotesize{Quantitative results on \textit{onshore} videos of the Singapore Maritime Dataset}}\vspace{-2mm}
	\setlength{\extrarowheight}{5pt}
	{\footnotesize
		\begin{tabular}{P{0.11\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.01\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}
				P{0.045\linewidth}}
			\cline{2-13}
			&\multicolumn{6}{c}{$Y^{\epsilon} = |Y - Y^{GT}|$ in pixels}&\multicolumn{6}{c}{$\phi^{\epsilon} = |\phi - \phi^{GT}|$ in degrees}\\
			\cline{2-7}
			\cline{9-13}
			&$\mu$&$\sigma$&Q25&Q50&Q75&Q95&~~~&$\mu$&$\sigma$&Q25&Q50&Q75&Q95\\
			\hline

			Ours
			&5.87&13.54&1.20&2.56&5.15&22.72&~~~&0.22&0.25&0.05&0.14&0.33&0.66

			\\
			\hline

			Li~\cite{li2021sea}
			&50.49&97.13&7.24&17.73&48.14&246.26&~~~&1.07&2.23&0.17&0.39&1.13&4.26

			\\
			\hline

			Liang~\cite{liang}
			&31.50&49.19&2.50&17.65&33.13&151.97&~~~&3.40&4.13&0.48&1.72&4.85&12.81

			\\
			\hline

			Jeong~\cite{jeong(roi)}
			&7.86&12.56&1.68&3.61&8.60&29.62&~~~&0.28&0.22&0.12&0.24&0.37&0.72

			\\
			\hline

			Gershikov~\cite{gershikov(1)}
			&119.03&93.29&30.62&99.83&185.74&284.41&~~~&5.09&5.37&1.19&2.89&7.20&17.20

			\\
			\hline
		\end{tabular}
	}
\end{table*}

\begin{table}[!h]
	\centering
	\caption{\label{tab_execution_time}\footnotesize{Mean computational time per frame on the SMD (Onboard and Onshore) (in milliseconds)}}\vspace{-2mm}
	\setlength{\extrarowheight}{5pt}
	{\footnotesize
		\begin{tabular}{
				P{0.24\linewidth}
				P{0.31\linewidth}
				P{0.31\linewidth}}
			\hline
			&Onboard&Onshore\\
			\hline

			Ours
			&107.6&101.4

			\\
			\hline

			Li~\cite{li2021sea}
			&486.5&426.0

			\\
			\hline

			Liang~\cite{liang}
			&73.7&79.4

			\\
			\hline

			Jeong~\cite{jeong(roi)}
			&90.6&95.7

			\\
			\hline

			Gershikov~\cite{gershikov(1)}
			&50.6&53.0

			\\
			\hline
		\end{tabular}
	}
\end{table}

\section{Conclusion}
\label{conclusion}
The sea horizon line plays a fundamental role in the overall maritime tracking system. Our focus was to detect in real-time the horizon line under different sea clutter without failing in the case where the horizon edge response is weak. The results we reported indicate the our method can detect the horizon lines with degraded edges, which is due to our filtering approach that does not affect horizon edges and leverages relevant properties of grown line segments extracted with a low edge magnitude threshold. While incorporating temporal information is useful in avoiding outlier detections, we believe that spatial information related to the relative positioning of EO sensors onboard the marine platform would increase the algorithm performance.
\printbibliography
\end{document}


text:
\section{Contribution}

The proposed approach introduced by \citet{zhang2018}, \citet{Siddiqui2019} and \citet{zhao2019} is to apply Deep Learning models to extract patterns not only from single point data but also from sky images. In particular, 2D inputs such as images can be treated by a specific type of neural network termed a CNN. The learning of such models is based on the use of filters which are trained to find task relevant patterns in an image. In parallel with other past meteorological data, such information can be given to an Artificial Neural Network (ANN) to forecast irradiance. The global network made of the CNN and the ANN can be trained simultaneously using supervised learning. The original contributions of this study are the following:
\vspace{-0.5\baselineskip}

\begin{itemize}
    \item \citet{Siddiqui2019} uses CNNs to extract patterns on sky images from the past 4 to 6 hours as general information about the cloud cover of that day to adapt forecasts up to 6h ahead. We, however, set up our protocol to learn how the clouds visible from the lens at time $t$ will induce the variability of the solar resource in the short-term upcoming future. In other words, we limit the forecasting time window to a narrower 20-min window from auxiliary data (past irradiance measurements, solar position of the sun and its sine / cosine transformations) and images taken at time $t$ and $(t-2)$min.

    \vspace{0.4\baselineskip}

    \item Instead of training and validating the model on distinct days, which does not let the model learn from same day past samples, we tried a new setting to generate both sets, which improves forecasting performances by around 10\%. With this approach, the model is validated on afternoon samples only, morning samples being in the training set. This way, the protocol is closer to real life applications where past samples of the same day are available for model fine-tuning.

    \vspace{0.4\baselineskip}

    \item To assess the model's performance we implemented the MSE forecast skill based on the smart persistence model, which is using the last in-situ pyranometric measurements only. A clear-sky index is computed with a concomitant model of clear-sky (here the ESRA model). This clear-sky index is then multiplied by the corresponding clear-sky prediction of the model to obtain the smart persistence forecast. This gives a reliable performance score to compare the proposed Deep Learning model with other approaches.

    \vspace{0.4\baselineskip}

    \item To increase the receptive field of the network without using dilated convolution \citep{Dilated_convolutions,Siddiqui2019}, we implemented a ResNet architecture \citep{Resnet}, which enables solvers to handle a deeper model.
\end{itemize}

\section{Approach}

\subsection{Dataset}
\label{dataset}

The chosen approach is to provide the Deep Learning model with sky images taken by hemispherical cameras on the ground (see Figure~\ref{fig:img_short_long}) and a range of metadata such as past irradiance measurements or the angular position of the sun.

\vspace{0.5\baselineskip}
The dataset used in this study originated from the SIRTA laboratory \cite{sirta}. Samples were collected over a period of seven months from March 2018 to September 2018. Images of the sky have RGB resolution of 768 by 1024 pixels. Samples are taken every two minutes and composed of two images with different exposures (referred to as long and short exposures, see Figure~\ref{fig:img_short_long}).

\begin{figure}[H]
\centering
\begin{subfigure}{0.23\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/20180509151800_01.jpg}

\end{subfigure}
\begin{subfigure}{0.23\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/20180509151800_03.jpg}

\end{subfigure}%
\caption{Images of the sky taken by a hemispherical camera with short and long exposures (source:~\cite{sirta})}
\label{fig:img_short_long}
\end{figure}

\vspace{-0.7\baselineskip}

The short exposition time provides more details on the region close to the sun, whereas the long exposure focuses on the rest of the sky, and in particular on distant clouds. In addition to this, the dataset has a range of metadata, in particular the Global Horizontal Irradiance (GHI) measurements which are averaged over a minute. Also, the angular position of the sun is available through the Zenith and Azimuthal angles.

\subsection{Network Architecture}
\label{network_architecture}

The proposed model follows the configuration presented in Figure~\ref{fig:network_architecture_2}. It is composed of two distinct networks merged into one which output the irradiance estimate. On the one side, a CNN made of ResNet units is used to extract features from sky images and on the other side, an ANN treats available auxiliary data (past irradiance measurements, angular position of the sun, etc). Both outputs are fed into another ANN, which integrates them to give its prediction.

\begin{figure}[H
\centering
\includegraphics[width=0.5\textwidth]{Figs/forecasting_CNN_model_4img.png}
\caption{Schematic network architecture}
\label{fig:network_architecture_2}
\end{figure}

\vspace{-0.5\baselineskip}
Hyperparameters  of the model such as the number of filters / nodes per layer or the number of layers are mostly hand-tuned given the model's performance. Another approach presented in Section~\ref{models_leaning} was used to select the number of filters per layer. In addition, the number of layers, and in particular the number of convolutional layers with a stride of 2 (the filter convolves around the input volume by shifting two units at a time) was set such that the output of the last convolutional layer was computed from each pixel of the input images.

\vspace{0.5\baselineskip}
The CNN architecture is composed of a set of convolutional layers with 32 filters in each and strides of 1 or 2. ResNet units using residual connections are implemented to increase the depth of the network. This way, if a set of convolutional layers does not improve the overall performance, the network is able to bypass the given 'residual connection' through a parallel identity connection (see Figure~\ref{fig:network_architecture_1}). The vectorisation layer reshapes the 2D output of the last convolutional layer into a 1D layer whose nodes are connected to dense layers with 512 and 64 nodes.

\vspace{0.5\baselineskip}
The parallel ANN is composed of two layers of 16 nodes plus a residual connection, which merges with the other network through a concatenation layer. Following this, another set of 2 densely connected layers integrate the incoming information to output the irradiance estimate through the final node.

\begin{figure}[H
\centering
\includegraphics[width=0.34\textwidth]{Figs/Forecasting_CNN_architecture_32.png}
\caption{Detailed network architecture}
\label{fig:network_architecture_1}
\end{figure}

\subsection{Training and Validation}
\label{training_and_validation}

Each sample given to the network is composed of four images in total, the short and the long exposure taken at time $t$ and $(t-2)$min. Each 2D input was originally downsized to a grey scale 150x150px resolution. In parallel, metadata are made of the past irradiance measurements at time $t$ and $(t-2)$min, the angular position of the sun at time $t$ defined by its Zenith and Azimuthal angles, plus the cosine and sine of those angles.

\vspace{0.5\baselineskip}
The assessment of the model was performed on a setting of the training and the validation sets composed of 16,000 and 4000 samples respectively. Samples collected from 8am to 7pm over the seven months period were randomly allocated to the validation and the training sets from two groups comprising data from distinct days. This way we prevent the model from using samples of the same day seen in the training set to adapt its forecasts on samples of the validation set.

\vspace{0.5\baselineskip}
In a second setting however, we tried to shift morning samples from the validation set to the training set to see how training the model on past samples of the same day improves its predictions. This configuration could be seen as closer to real life applications as data are streamed continuously and could be used to fine-tune the model given previous measurements.

\vspace{0.5\baselineskip}
The loss function used by the model as a reference to assess its own performance is the unregularised MSE. The learning rate associated with the minimisation of the loss was set to $10^{-3}$. The performance of the model was assessed using the forecast skill metrics based on the persistence of the clear-sky index on both training and validation set settings.

\vspace{0.5\baselineskip}
Hyperparameter tuning was performed to achieve the best forecasting performance on the 10-min ahead forecast. The same network architecture has been used to train models for the 2-min to 20-min ahead forecasts. The only difference is that the number of filters per convolutional layer was reduced to 16 for the 2-min and 4-min ahead forecasts due to a tendency to overfitting with 32 filters.

\section{Results}
\label{results}

Figure~\ref{fig:fs_different_time_windows} shows the performance of the model obtained with the same architecture for different time windows from 2 to 20-min. Contrary to~\cite{zhang2018}, the longer the forecast window, the better the forecast skill, until it reaches a plateau from the 12-min forecast window, which tempers the conclusion of~\cite{zhang2018} about a decrease of the skill score for longer time horizon forecasts.

\vspace{-0.5\baselineskip}

\begin{figure}[H
\centering
\includegraphics[width=0.38\textwidth]{Figs/fs_different_time_windows.png}
\caption{Forecasting skill on the validation set using the MSE based on the smart persistence model for different time windows, the training and validation sets being generated from different days}
\label{fig:fs_different_time_windows}
\end{figure}

Also, it is interesting to note that showing past data of the same day during the training strongly increases the model's performance with a skill score on the 10-min ahead forecast increasing by around 10\% from 0.4 to 0.44 for instance (see Table~\ref{table:results}). These results highlight the need to take same day samples of the past into account even for short term forecasting. This could be achieved through active learning or forecasting from a longer sequence of samples from the past (see~\cite{Siddiqui2019}).

\renewcommand{\arraystretch}{1.1}
\begin{table}[H]
\begin{center}
\begin{small}
\begin{tabular}{|c|ccc|}
\hline
Settings & 2-min & 10-min & 20-min\\
\hline
Distinct days & 0.20 & 0.4 & 0.43 \\
Validation on afternoon samples & \textbf{0.26} & \textbf{0.44} & \textbf{0.45} \\
\hline
\end{tabular}
\end{small}
\end{center}
\caption{Forecast Skill (MSE) on the validation set based on the smart persistence model}
\label{table:results}
\vskip -0.1in
\end{table}

\section{Model's learning assessment}
\label{models_leaning}

Often considered as black boxes, Deep Learning models foster research focusing on their interpretability. In that respect, different methods have been proposed to visualise models' learning, from intermediate activations and filter visualisation \citep{dl_with_python} to feature maps \citep{Zeiler2013}.

\vspace{-0.7\baselineskip}

\subsection{Intermediate activation visualisation}

\vspace{-0.8\baselineskip}

\begin{figure}[h!]
    \begin{minipage}[b]{0.55\linewidth}
        Visualising the layer activation of a trained model is useful for two reasons. Firstly, one can see what trained filters focus on for a given input. Figure~\ref{activation_visualisation} shows the transformation of figure~\ref{activation_visualisation_image}, once passed through 60 filters of the first convolutional layer. As we can see, some filters focus on specific patterns of the image such as the sun or the distant sky. However, some filters seem to look at useless areas of the image such as the black background or even simply return a blank output.

    \end{minipage}\hfill
    \begin{minipage}[b]{0.41\linewidth}
        \centering
        \includegraphics[width=0.9\textwidth]{Figs/response_filters_image.png}
        \caption{Example of an image of the sky, which is passed through the first layers to visualise intermediate activations (see Figure~\ref{activation_visualisation})}
        \label{activation_visualisation_image}
    \end{minipage}
\end{figure}

\vspace{-0.5\baselineskip}
This leads to the second use of this analysis technique, which is selecting the number of filters per layer in the network. Figure~\ref{activation_visualisation} indicates that the model did not make use of all 60 filters during the training. This gives a hint that decreasing the number of filters in this layer would not penalise the learning process, and would reduce the complexity of the model, i.e shorten the training time while preventing the model from overfitting. In the final configuration, the number of filters per convolutional layer was reduced to 32.

\begin{figure}[H
\centering
\includegraphics[width=0.46\textwidth]{Figs/response_filters.png}
\caption{2D output resulting from the convolution stage by filters of the first layer}
\label{activation_visualisation}
\end{figure}

\vspace{-2\baselineskip}

\subsection{Filter Visualisation}

As explained in \cite{dl_with_python}, the different patterns learnt by a model can be visualised by generating, for each filter of the network, the visual pattern that it responds to. To quantify the response of a filter to an image, one can define a loss function returning the average value of the output resulting from the convolution by the given filter. Starting from a noisy image, the visual pattern is then formed using gradient ascent on the pixel of the input to maximise the loss function measuring the response of the filter to that input.

\vspace{0.5\baselineskip}
Some images generated by this method for filters of the first, third and fifth convolutional layers are presented in Figure~\ref{fig:filters_first_layer}. One can notice that patterns spotted by filters of the first layers are rather abstract and difficult to interpret. However, as we go deeper into the network, corresponding filters respond to patterns of clouds. And the higher the convolutional layer, the more complex the corresponding images: a repetition of small patterns for images of the first and third layers, but complex shapes of cloud covers for the fifth layer. In particular, in Figure~\ref{fig:filters_first_layer}, some filters seem to focus on sparse cloud covers and others at more dense covering.

\vspace{-0.5\baselineskip}

\begin{figure}[H]
\centering
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/1layer_1filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/1layer_2filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/1layer_8filter.png}

\end{subfigure}%

\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/3layer_2filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/3layer_3filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/3layer_5filter.png}

\end{subfigure}%

\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/5layer_2filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/5layer_3filter.png}

\end{subfigure}
\begin{subfigure}{0.15\textwidth}
    \centering
    \includegraphics[width=1\linewidth]{Figs/Filters/5layer_5filter.png}

\end{subfigure}%

\caption{Tuned images leading to strong responses by filters of the first, third and fifth convolutional layers (first, second and third row respectively)}
\label{fig:filters_first_layer}
\end{figure}

\vspace{-2\baselineskip}

\section{Conclusions}

This study shows that the Deep Learning framework is a promising approach to irradiance forecasting with CNNs being able to extract relevant features from sky images. Also, integrating same day historical data in model predictions proves to be a sensible aspect of short-term forecasting.

\begin{acknowledgements}

The authors would like to acknowledge SIRTA for providing the sky images and irradiance measurements used in this study. We also would like to thank Prof. Philippe Blanc for his valuable advice. This research was supported by Engie, EPSRC and the University of Cambridge.

\end{acknowledgements}

\vspace{-2\baselineskip}

redpajama_set_name: RedPajamaArXiv