\section{Introduction} \vspace{-0.2cm}

Recent advances in 3D hardware sensors (e.g., 3D scanners, LiDARs and RGB-D cameras) have improved the accessibility of point clouds, which can intuitively reflect the spatial geometric information of objects in 3D space.
Recently, 3D semantic segmentation of point clouds, which aims to label each point with a corresponding category information, has drawn increasing attention due to its broad applicability to the fields of general environment perception.
It enables various higher-level potential applications, including autonomous driving, human-computer interaction, virtual reality and robotics.

With the significant progress of deep learning on grid-shaped images/videos~\cite{NIPS2012_imagenet,resnet}, recent studies began to explore how to apply convolutional neural networks (CNN) on irregular-structured point clouds~\cite{qi2017pointnet,qi2017pointnet++}.
For example, some previous methods converted point clouds data to regular data representation by employing multi-view images~\cite{su2015multi,lawin2017deep} or voxels~\cite{maturana2015voxnet,wu20153d,qi2016volumetric}, and then CNN-like operations (e.g., 2D CNN and 3D CNN) can be performed on unstructured point clouds.
Although these methods have made some progress, the projection and voxelization steps inevitably lead to information loss of original 3D data.
Furthermore, as a pioneering work, PointNet~\cite{qi2017pointnet} was proposed to directly process point clouds  without intermediate data conversion. It employed the shared multi-layer perceptron (MLP) to learn features of each point and then a symmetric function is used to aggregate global features. This network can achieve permutation invariance for unordered 3D point clouds, but cannot well consider the local geometric information between points (including direction, shape and topological structure).
Subsequently, several dedicated neural modules have been proposed to aggregate local geometric information of point clouds for further boosting the semantic segmentation performance, such as Pointnet++~\cite{qi2017pointnet++}, RS-CNN~\cite{liu2019relation},
Geo-CNN~\cite{yang2019modeling}, PointWeb~\cite{zhao2019pointweb} and ShellNet~\cite{zhang2019shellnet}.
For example, PointNet++ adopted a hierarchical neural network to process a set of points sampled in a metric space~\cite{qi2017pointnet++};  RS-CNN defined the low-level relation as a compact feature vector with 10 channels, then contextual shape-aware representation was learned for all points~\cite{liu2019relation}.

However, most of these approaches always aggregate/mix these geometric information together to explore the relationship between a point and its neighbor points. Such aggregated methods essentially belong to the isotropic filtering on all operated points therein, which tends to lose the information of 3D geometric structures.

In this paper, we propose a novel spatial transformer point convolution (STPC) method to achieve anisotropic convolution filtering on point clouds, and finally boost the performance of 3D semantic segmentation task.
Given a 3D point clouds sample, our proposed STPC framework can predict the point-wise categorization in an end-to-end fashion.
Inspired by the classic bag-of-words model~\cite{bagoftheword1,bagoftheword}, we specifically build a spatial direction dictionary to represent these latent spatial directions/geometric components of local point clouds data.
By virtue of the dictionary, these unordered and unstructured neighbor points can be coordinated into the canonical dictionary space by using sparse direction encoding.
Subsequently, in the transformed dictionary space, a typical image-like convolution operation can be applied to perform the anisotropic filtering process for robustly capturing these finer geometric structure information of point clouds.
Both dictionary learning and encoding processes are encapsuled into a network module, and the entire spatial transformer point convolution is also jointed to optimize in an end-to-end neural network.
The proposed STPC framework can better capture these subtle geometric structure information, especially for those finer variances of local regions in 3D point clouds.

We summarize the main contributions as three folds:
\begin{itemize}
	\vspace{-0.15cm}\item We propose a novel spatial transformer point convolution (STPC) framework to deal with semantic segmentation on point cloud data, which successfully performs anisotropic convolution filtering on unstructured data like the standard convolution on images.
	\vspace{-0.15cm}\item We propose a direction dictionary induced spatial direction encoding method, which transforms unordered neighbor points into a latent atom-coordinated system and further well encode those subtle structure variances of local neighbor points.
	\vspace{-0.15cm}\item Comprehensive evaluations on three point cloud datasets (including S3DIS, Semantic3D, SemanticKITTI)  demonstrate the superiority of our proposed STPC when compared with other state-of-the-art methods in the point cloud semantic segmentation problem.
\end{itemize}\vspace{-0.3cm}
\section{Related Work}\vspace{-0.2cm}
In this section, we focus on point clouds data analysis and briefly review three types of existing deep learning methods according to their underlying technologies.

\textbf{Multiview-based convolution methods:} Some existing methods \cite{su2015multi,lawin2017deep,feng2018gvcnn} project point clouds into a set of renderings from different viewpoints and then the standard convolutional neural network can be applied. MVCNN\cite{su2015multi} projected a 3D object into multiple views and extracted the corresponding view-wise features, then simply max-pools multi-view features into a global descriptor for accurate object recognition. Felix et al.\cite{lawin2017deep} projected a 3D point clouds onto 2D planes from multiple virtual camera views. Then, a multi-stream FCN was used to predict pixel-wise scores on synthetic images. The final semantic label of each point was obtained by fusing the reprojected scores over different views. However, these view-based convolution methods lost spatial information during projection, and can not fully exploit the underlying geometric and structural information. Also, the performance of multi-view semantic segmentation methods is sensitive to viewpoint selection and occlusions.

\textbf{Voxelization-based convolution methods:} Volumetric representation naturally preserves the neighborhood structure of 3D point clouds. Its regular data format allows direct application of standard 3D convolutions~\cite{maturana2015voxnet,wu20153d,tchapmi2017segcloud,huang2016point,le2018pointgrid}. VoxNet~\cite{maturana2015voxnet} introduced a volumetric occupancy network to achieve robust 3D object recognition. 3D ShapeNet~\cite{wu20153d} was proposed to represent 3D shape by a probability distribution of binary variables on voxel grids. To alleviate the problem of computation and memory footprint growing cubically with the resolution, Kd-trees~\cite{klokov2017escape} and octrees~\cite{riegler2017octnet,wang2017o-cnn} were introduced to efficiently model point clouds. However, the voxelization step inherently causes discretization artifacts and information loss. It is non-trivial to select an appropriate grid resolution in practice.

\textbf{Point-based methods.} The point-based method directly takes raw point clouds as input without converting them to other formats. As a pioneering work, PointNet~\cite{qi2017pointnet} used shared multi-layer perceptron to learn features of each point and then a symmetric function was used to aggregate global features. This network can achieve permutation invariance for unordered 3D point clouds but the local geometric information between points can not be captured. Subsequently,
neighbouring feature pooling methods~\cite{qi2017pointnet++,zhao2019pointweb,zhang2019shellnet}, graph message passing methods~\cite{wang2018local,wang2019graph,wang2019dynamic}, kernel-based convolution methods~\cite{atzmon2018point,groh2018flex-convolution,hua2018pointwise,thomas2019kpconv}, and attention-based aggregation methods~\cite{xie2018attentional,yang2019modeling} were proposed to learn per-point local features. Although these methods have shown promising results, they did not implement the anisotropic filtering on point clouds.
To generate anisotropic filtering on point clouds, PointCNN~\cite{li2018pointcnn} learnt a $\mathcal{X}$-transformation from the input points to weight the input features associated with the points and permute the points into a latent and potential canonical order. TangentConv~\cite{tatarchenko2018tangent} was proposed by projecting local neighbor points to local tangent planes and processing them with 2D convolutions. A-CNN~\cite{komarichev2019cnn} can order neighboring in a
clockwise/counterclockwise manner on a tangent plane and then apply annular convolutions. FPConv~\cite{lin2020fpconv} performed a local flattening by learning a weight map to softly project surrounding points onto a 2D grid.
Our STPC is different from the above anisotropic methods in three folds: i) introduces a flexible spatial direction dictionary as the latent coordinate system (each atom may be viewed as a coordinate direction), ii) uses the inspirit of sparse coding and transformer to sparsely encode neighbor points into the dictionary system, and iii) filter kernels are defined on latent coordinate directions for anisotropic convolution.

\vspace{-0.2cm}
\section{Our Approach}\vspace{-0.2cm}

\subsection{Overview}

The standard convolution on images/videos is conditioned to grid shape region, and performs anisotropic filtering in order to encode subtle texture appearances. The point clouds data are often embedded in an unstructured and non-regular 3D space, and thus the standard convolution kernels are difficult to be defined on point clouds data due to its irregular property.
Most existing convolution methods on point clouds ~\cite{qi2017pointnet,qi2017pointnet++} take the sum-/max-aggregation on local neighbor points without considering the local space topology of points therein.
As a contrast, we attempt to construct an \textit{anisotropic filtering} on 3D point clouds to better encode these finer characteristics of points in each local 3D spatial region. To this end, we introduce spatial direction dictionaries to coordinate these nearest points into latent normal spaces, where each atom of the dictionary implicitly defines a spatial direction.
The previous filter operated on all nearest points can be untied to adapt to different atoms, which naturally results to anisotropic convolution filtering analogous to the standard convolution. This anisotropic filtering can be named as \textit{spatial transformer point convolution (STPC)}.

Given one reference point $p_i$, we can sample the $K$ nearest neighbor points as a set $\mcN(p_i)=\{p_k|k=1,\cdots, K\}$. Obviously, each point has 3D coordinate value, denoted as $\p_k=(x_k,y_k,z_k)^\intercal} \def\st{\text{s.t.~}$\footnote{Here the lowercase $p$ denotes one general point, while the bold-type $\p$ denotes its spatial x-y-z coordinates. }. Point clouds often carry some additional attributes, e.g., RGB color or intermediate learnt features. We denote the attributes of point $p_{i}$ as a vector $\x_i$. Next, we perform an anisotropic
filtering on the point set $\mcN(p_i)$ to produce a robust response of the reference point $p_i$. Below we define one-time spatial transformer point convolution process, which can be encapsulated into an STPC module as one layer. Formally,
\begin{align}
\x_i^{l+1}& \longleftarrow \psi_\text{conv}([\tbx_1^l,\tbx_2^l,\cdots,\tbx_M^l]),\label{eqn:x_ilp1}\\
\tbx_m^l &\longleftarrow \mcT_{\text{f}} (\{(\alpha_{km}^l,\mcF_\text{f}^l(\p_k,\x_k^l))| p_k\in\mcN(p_i)\}),\quad m=1,\cdots M,\label{eqn:tbx_ml}\\
\bm{\alpha}_k^l:&[\alpha^l_{k1},\alpha^l_{k2},\cdots, \alpha^l_{kM}]\longleftarrow\mcT_{\text{s}}(\A^l,\mcF_\text{s}(\p_k)),\quad p_k\in\mcN(p_i), \label{eqn:alpha_jl}\\
\A^l&:[\a_1^{l},\a_2^{l},\cdots,\a_M^{l}]\longleftarrow \mcF_\text{a}(\a_1^{l-1},\a_2^{l-1},\cdots,\a_M^{l-1}),\label{eqn:A_l}
\end{align}
where $l$ is the layer number of neural network, $M$ is the number of atoms $\{\a_1,\a_2,\cdots,\a_M\}$ in the spatial direction dictionary $\A$, $\bm{\alpha}_k$ is the encoding coefficient w.r.t the point $p_k$, $\tbx_m$ is the accumulated component at the $m$-th atom for all neighbor points of $p_i$, $\x_i^{l+1}$ is the anisotropic convolution response of point $p_i$. In the above formulas,
\begin{itemize

	\item Direction dictionary leaning $\mcF_\text{a}$. Due to spatial encoding, this dictionary should emphasize spatial information (source from points positions), rather than attribute information. Meantime, the dictionary should be dynamic and depended on prior states by accompanying with stacked modules.
	\item Spatial direction encoding $\mcT_\text{s}$. The point $p_k$ is represented with the atoms of the direction dictionary by using sparse coding or nearest neighbor reconstruction. The representation coefficients describe the latent coordinate information of this point in the dictionary space. The spatial transformation $\mcF_\text{s}$ is used to extract spatial features from 3D coordinates.
	\item Spatial transformer $\mcT_\text{f}$. The neighbor points set $\mcN(p_i)$ is projected onto atom points, and forms a regular coordinate system defined by $M$ atoms. The feature transformation $\mcF_\text{f}$ is used to learn point features/attributes.
	\item Anisotropic convolution (i.e., STPC) $\psi_\text{conv}$. As the dictionary implicitly defines a regular coordinate system, the spatial transformed features $\{\tbx_m|m=1,\cdots, M\}$ just corresponds to the atoms therein. Hence, we could perform different filtering for each of them, which really works like the standard convolution on different relative spatial positions.
\end{itemize}

\begin{figure}[!t]
	\centering
	\includegraphics[width=1.0\linewidth]{framework6-7.png}
	\caption{The overview of spatial transformer point convolution (STPC) module. Here the direction dictionary learning, spatial direction encoding, spatial anisotropic filtering can be encapsulated into the STPC module, and jointly learnt in an end-to-end fashion. The STPC module can finally achieve an anisotropic convolution filtering on point clouds for capturing those finer variances of local region.

	} \vspace{-0.6cm}
	\label{conv}
\end{figure}
The above convolution process can be framed into a STPC network module as shown in Fig.~\ref{conv}, whose implementation details are introduced in section~\ref{aaa}, ~\ref{bbb} and ~\ref{ccc}.

\vspace{-0.2cm}
\subsection{Direction Dictionary Learning}\label{aaa} \vspace{-0.2cm}
In order to produce consistent responses under different permutation layouts, most existing methods mix these geometric information together to explore the relationship between a point and its neighbor points.
Intuitively, these geometric/direction information of point clouds lie in a dense space with several latent components/directions.
Inspired by the idea of bag-of-words, here we learn a spatial direction dictionary in the $c$-dimensional space, which can be then used to sparsely quantize local neighborhood structure  of point clouds.
Specifically, we build a spatial direction dictionary $\A=[\a_{1}, \a_{2}, \dots, \a_{M}]$ with $M$ atoms through random initialization, where each atom $\a_{m} \in \mathbb{R}^c $ refer to a spatial direction information.
Here the initialized dictionary can be noted as $\A^0=[\a_{1}^0, \a_{2}^0, \dots, \a_{M}^0]$.
Due to the sparsity of point clouds, a certain number of atoms are sufficient to encode each neighborhood region, and the redundancy and expensive computing costs can be also avoided in practice.
With the prior states of direction dictionary in previous layer, we can dynamically optimize/update it while satisfying the constrain of atom decorrelation.
\begin{equation}
\begin{aligned}
&\min \sum_{p=1}^{M}\sum_{q\neq{p}}^{M}\mathcal{S}_{(\a_p^l,\a_q^l)}, \qquad \text{s.t.} \ \ \A^l=\mcF_\text{a}( \A^{l-1}),
\label{eq:eps4}
\end{aligned}
\end{equation}
where $\A^{l}: [\a_1^{l},\a_2^{l},\cdots,\a_M^{l}] $ and $\A^{l-1} : [\a_1^{l-1},\a_2^{l-1},\cdots,\a_M^{l-1}]$ denote the states of direction dictionary in the $l$-th and $l-1$-th layers, the dictionary learning function $\mcF_\text{a}$ is implemented by a shared multi-layer perceptron.
In the direction dictionary learning process, we use the cosine similarity to measure the similarity between each of them.
\begin{equation}\label{eq:eps5}
\mathcal{S}_{(\a_p^l,\a_q^l)}=\frac{\a_p^l\cdot\a_q^l}{\|\a_p^l\|\|\a_q^l\|}, \qquad p,q=1,2,\dots,M,
\end{equation}
where the greater the cosine similarity $\mathcal{S}_{(\a_p^l,\a_q^l)}$ is,  the more similar the spatial directions between $\a_p^l$ and $\a_q^l$ are. The constraint is that the atoms of the spatial direction dictionary are discrete and not similar between them.
We can iteratively learn the spatial direction dictionary by minimizing the above constraint, and then the atoms can discretely represent different spatial  directions in the high-dimensional dictionary space.

\vspace{-0.2cm}
\subsection{Spatial Direction Encoding}\label{bbb}\vspace{-0.2cm}
For a single point, its x-y-z coordinates can not fully explore the local geometric information with its neighbor points, which is important for the semantic segmentation task. Several recent works\cite{yang2019modeling,liu2019relation,hu2019randla} have explicitly embedded the low-level elements of neighbor points to encode local features.
These above methods focus on expressing local structure relationships by employing multiple low-level elements, but the distinguishable spatial direction information of each point would be drowned out.
Therefore, we introduce a spatial direction encoding method to capture and represent implicit geometric structures and distinguishable direction information.
For the $i$-th point, its spatial direction information can be represented as $[\d_1, \d_2,\dots, \d_K]\in\mathbb{R}^{K\times c}$, and each of them can be encoded as:
\begin{equation}\label{eq:eps3}
\begin{aligned}
\d_{k} = \mcF_\text{s}(\p_i-\p_{k}), \ \ \ \  p_k \in \mcN(p_{i}),
\end{aligned}
\end{equation}
where the point set $\{p_{k}|k=1, 2, \dots, K\}$ refers to the $K$ nearest neighbor points of the point $p_{i}$,
$c$ is the feature dimension of $\d_{k}$,
and the spatial transformation $\mcF_\text{s}$ is performed by the shared multi-layer perceptron for extracting the spatial direction information.

In the $l$-th convolution layer, we can encode the local spatial information of point $p_i$ by employing the direction dictionary $\A^{l}$ and spatial direction information $[\d_k | k=1, 2, \dots, K]$.
For one nearest point $p_k$ of the point $p_i$, we can encode it with the atoms of the  direction dictionary by using sparse coding.
Analogous to the self-attention mechanism~\cite{vaswani2017attention} in the NLP field, $\d_{k}$ and $\a_m^{l}$ refer to a query and one dictionary atom, respectively.
For describing the latent coordinate information of the point $p_k$, we can get the representation coefficient by computing the correlation between spatial direction information $\d_{k}$ and the atom $\a_m^{l}$.
Specifically, we use the cosine similarity as Eqn. (\ref{eq:eps5}) to measure the correlation between them and an exponential function is used to make the representation coefficients more sparse.
The  spatial encoding process $\mathcal{T}_{\text{s}}$ can be expressed as computing the direction-dependent relationship between $p_k$ and each atom $\a_m^{l} \in \A^l$:
\begin{equation}\label{eq:eps8}
\alpha_{km}^l=\frac{\exp\left(\mathcal{S}_{(\d_{k},\a_m^l)}\right)}{\sum_{m=1}^{M}\exp\left(\mathcal{S}_{(\d_{k},\a_m^l)}\right)},
\end{equation}
where $\alpha_{km}^l$ may be viewed as the attention score/coefficient of the point $p_k$ w.r.t the atom $\a_m$. Thus the spatial information of the point $p_k$ in the dictionary space is encoded as $\bm{\alpha}_k^l=[\alpha^l_{k1},\alpha^l_{k2},\cdots, \alpha^l_{kM}]$.
Besides, the study on sparse coding~\cite{yang2009linear} demonstrates that the sparsity of representation coefficients is important to overcomplete dictionaries. Hence, we discard those trivial coefficients in $\bm{\alpha}$ by setting a small threshold value (i.e., $\tau$ = 0.01), which results into sparse representation in practice.

\vspace{-0.2cm}
\subsection{Spatial Anisotropic Filtering}\label{ccc}\vspace{-0.2cm}

To encode the unordered local information of the point $p_i$, we introduce a spatial transformer $\mathcal{T}_{\text{f}}$ for projecting the neighbor points set $\mathcal{N}(p_i)$ into a canonical order dictionary space. After obtaining the attention score, the feature of all neighbor points $p_k\in\mcN(p_i)$ can be transformed into the canonical ordered dictionary space.
\begin{equation}\label{eq:eps9}
\begin{aligned}
\tbx_m^l &= \mcT_{\text{f}} (\{(\alpha_{jm}^l,\mcF_\text{f}^l(\p_k,\x_k^l))| p_k\in\mcN(p_i)\}) \\
 &= \sum_{k=1}^{K} \alpha_{km}^l
\mcF_\text{f}^l(\p_k,\x_k^l), \quad m=1,\cdots M,
\end{aligned}
\end{equation}
where the feature transformation $\mcF_\text{f}^l(\p_k,\x_k^l)$ refers to the feature extracting process by employing the shared multi-layer perceptron, the x-y-z coordinates $\p_k$ and the geometry-aware features $\x_k^l$ of the point $p_{k}$ are used as inputs.
Thus, any disordered neighborhood can be transformed into the canonical ordered dictionary space, and then we can perform an anisotropic convolution filtering $\psi_\text{conv}$ to produce a robust response of the point $p_i$.
\begin{equation}
\begin{aligned}
\x_i^{l+1} &= \psi_\text{conv}([\tbx_1^l,\tbx_2^l,\cdots,\tbx_M^l]) \\

           & =\w_1*\tbx_1^l + \w_2*\tbx_2^l + \dots + \w_M*\tbx_M^l, \\
\end{aligned}
\end{equation}
where $[\tbx_1^l,\tbx_2^l,\cdots,\tbx_M^l]$ denotes the spatial transformed features with the ordered structure in the dictionary space, $[\w_1, \w_2, \dots, \w_M]$ refers to the learnt weights for performing anisotropic convolution.

\vspace{-0.2cm}
\section{Experiments}\vspace{-0.2cm}

\subsection{Experiment Setting}

To demonstrate the efficiency of our proposed STPC method, we conduct semantic segmentation experiments on three large-scale point clouds datasets, including S3DIS~\cite{armeni2017joint}, Semantic3D~\cite{hackel2017semantic3d} and SemanticKITTI~\cite{behley2019semantickitti}.
Given an point clouds as input, the widely-used encoder-decoder architecture is used to predict the point-wise label in an end-to-end way.
In the encoding stage, we stack five STPC modules to extract per-point features, each of which is followed by a random sampling operation~\cite{hu2019randla} to reduce the size of the point clouds.
In the decoding stage, five upsampling layers are stacked to propagate features by employing the nearest point interpolation technique.
Meanwhile, skip connection is applied to concatenate the upsampled features with the intermediate features produced by encoding layers.
To predict the categories of input point clouds, we stack several fully-connected layers on this responses, and use the cross-entropy loss during training.
Follow the same protocols as in~\cite{hu2019randla}, we split piont cloud data into training set and testing set, and adopt the same original features as inputs (3D coordinates and color information for S3DIS and Semantic3D, and only 3D coordinates for SemanticKITTI). We adopt the mean IoU (mIoU), mean class Accuracy (mAcc) and Overall Accuracy (OA) over the total classes as the standard metrics.
The learning rate is initialized at 0.01 and decreases by 5\% after each epoch. The number of nearest points K is set as 16.
As in RandLA-Net~\cite{hu2019randla}, we use the Adam optimizer with default parameters and the same data processing methods. Our learnt spatial direction dictionary contains 25 atoms (i.e., $M=25$), the dimension of each atom is 16 (i.e., $c=16$), each of which represents a different spatial direction.
For smooth updates, the learning rate of direction dictionary learning is set to 0.01 times the global network.
All models in our experiment are trained and tested based on an NVIDIA RTX2080Ti GPU.

\begin{table}[!t]
	\caption{Comparison of semantic segmentation performance on the S3DIS dataset (Area 5)}
	\label{s3dis_Area_5}
	\centering
	\resizebox{\textwidth}{!}{
		\begin{tabular}{rcccccccccccccccc}
			\toprule
			Method&OA(\%)&mAcc(\%)&mIoU(\%)&ceil.&floor&wall&beam&col.&wind.&door&table&chair&sofa&book.&board&clut. \\
			\midrule
			PointNet~\cite{qi2017pointnet}&-&49.0&41.1&88.8&97.3&69.8&0.1&3.9&46.3&10.8&58.9&52.6&5.9&40.3&26.4&33.2\\
			SegCloud~\cite{tchapmi2017segcloud}&-&48.9&57.4&90.1&96.1&69.9&0.0&18.4&38.4&23.1&70.4&75.9&40.9&58.4&13.0&41.6\\
			TangentConv~\cite{tatarchenko2018tangent}&-&62.2&52.6&90.5&97.7&74.0&0.0&20.7&39.0&31.3&77.5&69.4&57.3&38.5&48.8&39.8\\
			SPGraph~\cite{landrieu2018large}&86.4&66.5&58.0&89.4&96.9&78.1&0.0&\textbf{42.8}&48.9&61.6&84.7&75.4&69.8&52.6&2.1&52.2\\
			PointCNN~\cite{li2018pointcnn}&85.9&63.9&57.2&92.3&98.2&79.4&0.0&17.6&22.8&62.1&74.4&80.6&31.7&66.7&62.1&56.8\\
			PointConv~\cite{wu2019pointconv}&85.4&64.7&58.3&92.8&96.3&77.0&0.0&18.2&47.7&54.3&\textbf{87.9}&72.8&61.6&65.9&33.9&49.3\\
			PAT~\cite{yang2019modeling}&-&70.8&60.1&\textbf{93.0}&\textbf{98.5}&72.3&\textbf{1.0}&41.5&\textbf{85.1}&38.2&57.7&83.6&48.1&67.0&61.3&33.6\\
			PointWeb~\cite{zhao2019pointweb}&87.0&66.7&60.3&92.0&\textbf{98.5}&79.4&0.0&21.1&59.7&34.8&76.3&88.3&46.9&69.3&64.9&52.5\\
			KPConv rigid~\cite{thomas2019kpconv}&-&70.9&65.4&92.6&97.3&81.4&0.0&16.5&54.5&\textbf{69.5}&80.2&90.1&66.4&74.6&63.7&58.1\\
			KPConv deform~\cite{thomas2019kpconv}&-&72.8&\textbf{67.1}&92.8&97.3&\textbf{82.4}&0.0&23.9&58.0&69.0&81.5&\textbf{91.0}&75.4&\textbf{75.3}&\textbf{66.7}&\textbf{58.9}\\

			Our STPC&\textbf{88.5}&\textbf{75.3}&66.6 &91.7 &96.7 &82.1 &0.0 &37.3 &64.7 &52.6 &79.7 &89.2 &\textbf{76.3} &72.5 &66.5 &56.1 \\
			\bottomrule
	\end{tabular}}
\end{table}
\begin{table}[!t] \vspace{-0.4cm}
	\caption{Comparison of semantic segmentation performance on the Semantic3D dataset (reduced-8)}
	\label{Semantic3d}
	\centering
	\resizebox{\textwidth}{!}{
		\begin{tabular}{rcccccccccc}
			\toprule
			Method&mIoU(\%)&OA(\%)&man-made&natural&high veg.&low veg.&buildings&hard scape&scanning art.&cars\\
			\midrule
			SegCloud~\cite{tchapmi2017segcloud} &61.3 &88.1 &83.9 &66.0 &86.0 &40.5 &91.1 &30.9 &27.5 &64.3\\
			ShellNet~\cite{zhang2019shellnet} &69.3 &93.2 &96.3 &90.4 &83.9 &41.0 &94.2 &34.7 &43.9 &70.2\\
			GACNet~\cite{wang2019graph} &70.8 &91.9 &86.4 &77.7 &\textbf{88.5} &\textbf{60.6} &94.2 &37.3 &43.5 &77.8\\
			SPGraph~\cite{landrieu2018large} &73.2 &94.0 &97.4 &92.6 &87.9 &44.0 &83.2 &31.0 &63.5 &76.2\\
			KPConv rigid~\cite{thomas2019kpconv} &74.6 &92.9 &90.9 &82.2 &84.2 &47.9 &94.9 &40.0 &\textbf{77.3} &\textbf{79.7}\\
			RGNet~\cite{truong2019fast} &74.7 &94.5 &97.5 &\textbf{93.0} &88.1	&48.1 &94.6	&36.2 &72.0 &68.0\\
			RandLA-Net~\cite{hu2019randla} &\textbf{77.4} &\textbf{94.8} &95.6 &91.4 &86.6 &51.5 &\textbf{95.7} &\textbf{51.5} &69.8 &76.8\\
			Our STPC &76.2 &\textbf{94.8} &\textbf{97.9} &\textbf{93.0} &88.1 &49.3 &95.4 &47.3 &59.0 &79.6\\
			\bottomrule
	\end{tabular}} \vspace{-0.2cm}
\end{table}
\vspace{-0.2cm}
\subsection{Results and Comparisons}\vspace{-0.2cm}

\textbf{S3DIS dataset}\cite{armeni2017joint}: Table~\ref{s3dis_Area_5} shows the performance of our STPC method and comparisons with several state-of-the-arts on OA, mAcc and mIoU metrics. our STPC method can significantly outperform four baselines: 2.1\% over SPGraph~\cite{landrieu2018large}, 2.6\% over PointCNN~\cite{li2018pointcnn}, 3.1\% over PointConv~\cite{wu2019pointconv} and 1.5\% over PointWeb~\cite{zhao2019pointweb} in term of OA score.
When compared with these existing point clouds semantic segmentation approaches, our method in term of mAcc score also achieves the best result, representing improvement of 26.3\%, 13.1\%, 11.4\%, 4.5\%, and 2.5\% over PointNet~\cite{qi2017pointnet},  TangentConv~\cite{tatarchenko2018tangent},
 PointCNN~\cite{li2018pointcnn}, PAT~\cite{yang2019modeling}, KPConv deform~\cite{thomas2019kpconv}.
The proposed STPC method can also obtain a better performance than these exiting methods on mIoU score, except for KPConv deform method~\cite{thomas2019kpconv}.
Furthermore, we also report the mIoU scores for each class in Table~\ref{s3dis_Area_5}.
Generally, when compared with these baselines, our STPC method shows the comparable performance, e.g, 76.3\% \textit{vs} 75.4\%~\cite{thomas2019kpconv} for sofa, 89.2\% \textit{vs} 91.0~\cite{thomas2019kpconv} for chair, 66.5\% \textit{vs} 64.9~\cite{zhao2019pointweb} for board, and 82.1\% \textit{vs} 79.4\%~\cite{li2018pointcnn} for wall.
It demonstrates that our STPC method performs very well on predicting the point-wise labels by considering the anisotropic convolution filtering technology on 3D point clouds data.

\textbf{Semantic3D dataset}~\cite{hackel2017semantic3d}:
We report the results of comparisons between the baselines and our STPC method
on the Semantic3D dataset (reduced-8). As shown in Table~\ref{Semantic3d}, comparisons with previous
methods~\cite{tchapmi2017segcloud,zhang2019shellnet,thomas2019kpconv,truong2019fast} demonstrate that our proposed
STPC obtains the best performance, achieving improvements
of 0.3\% over RGNet~\cite{truong2019fast}, 1.9\% over  KPConv rigid~\cite{thomas2019kpconv}, 1.6\% over ShellNet~\cite{zhang2019shellnet} and 6.7\% over SegCloud~\cite{tchapmi2017segcloud} in term of OA score.
The STPC achieves 76.2\% in term of OA score, which can outperform most existing methods~\cite{tchapmi2017segcloud,zhang2019shellnet,thomas2019kpconv,truong2019fast} , except the RandLA-Net method~\cite{hu2019randla}.
For example, when compared with the classic SPGraph~\cite{landrieu2018large}, we can boost the semantic segmentation performance by 0.8 and 3.0\% in terms of OA and mIoU scores, respectively.
Furthermore, when predicting the labels for each category, the mIoU scores for man-made and natural can achieve the best performance with 97.9\% and 93.0\% segmentation results.
This indicates that our SPTC can better capture these finer local features of point clouds and then improve the discriminative capability of semantic segmentation network.

\textbf{SemanticKITTI dataset}~\cite{behley2019semantickitti}:
We finally compare the proposed STPC method with two families of recent point clouds semantic segmentation approaches, including point based methods~\cite{qi2017pointnet,landrieu2018large,qi2017pointnet++,tatarchenko2018tangent,hu2019randla} and projection based approaches~\cite{wu2018squeezeseg,wu2019squeezesegv2,behley2019semantickitti,milioto2019rangenet++}.
The quantitative comparison of segmentation results is reported in Table~\ref{SemanticKITTI}.
When compared with all point based methods and projection based approaches, the STPC can achieve the best performance in term of mIoU score, and also obtain the better results in most categories (11 out of 19).
Especially in the motorcyclist category, the mIoU results of all these existing methods are quite low (e.g., 0.0\% with PointNet~\cite{qi2017pointnet}, SPLATNet~\cite{su2018splatnet} and PointNet++~\cite{qi2017pointnet++}, 0.9\% with SqueezeSeg~\cite{wu2018squeezeseg}, 7.2\% with RandLA-Net~\cite{hu2019randla}), but we can achieve 15.3\% segmentation result, which is much higher than the second-ranked TangentConv method~\cite{tatarchenko2018tangent} by improving 7.2\%.
Compared with the best projection based approach
(i.e., RangeNet53++~\cite{milioto2019rangenet++}), we outperform it by 2.4\% in term of  mIoU score, and obtain the better results in 14 categories.  Compared with the best point based method (i.e., RandLA-Net~\cite{hu2019randla}), we outperform it by 0.7\% in term of  mIoU score, and achieve the better results in 16 categories.
Some visualization results between our STPC method and the compared RandLA-Net~\cite{hu2019randla} can be found in Fig.~\ref{semantickitti_compare}.
These above experimental results indicate that our proposed STPC method can better predict the point-wise category information of point clouds data by improving the predictive capability of the network.

\begin{table} \vspace{-0.4cm}
	\caption{Comparison of semantic segmentation performance on the SemanticKITTI dataset}
	\label{SemanticKITTI}
	\centering
	\resizebox{\textwidth}{!}{
		\begin{tabular}{rcccccccccccccccccccccc}
			\toprule
			Methods &Size &\rotatebox{90}{mIoU(\%)} &\rotatebox{90}{road} &\rotatebox{90}{sidewalk} &\rotatebox{90}{parking} &\rotatebox{90}{other-ground} &\rotatebox{90}{building} &\rotatebox{90}{car} &\rotatebox{90}{truck} &\rotatebox{90}{bicycle} &\rotatebox{90}{motorcycle} &\rotatebox{90}{other-vehicle} &\rotatebox{90}{vegetation} &\rotatebox{90}{trunk} &\rotatebox{90}{terrain} &\rotatebox{90}{person} &\rotatebox{90}{bicyclist} &\rotatebox{90}{motorcyclist} &\rotatebox{90}{fence} &\rotatebox{90}{pole} &\rotatebox{90}{traffic-sign} \\
			\midrule
			SqueezeSeg~\cite{wu2018squeezeseg} &\multirow{5}{1cm}{64*2048 \centering{pixels}} &29.5 &85.4 &54.3 &26.9 &4.5 &57.4 &68.8 &3.3 &16.0 &4.1 &3.6 &60.0 &24.3 &53.7 &12.9 &13.1 &0.9 &29.0 &17.5 &24.5\\
			SqueezeSegV2~\cite{wu2019squeezesegv2} & &39.7 &88.6 &67.6 &45.8 &17.7 &73.7 &81.8 &13.4 &18.5 &17.9 &14.0 &71.8 &35.8 &60.2 &20.1 &25.1 &3.9 &41.1 &20.2 &36.3\\
			DarkNet21Seg~\cite{behley2019semantickitti} & &47.4 &91.4 &74.0 &57.0 &26.4 &81.9 &85.4 &18.6 &26.2 &26.5 &15.6 &77.6 &48.4 &63.6 &31.8 &33.6 &4.0 &52.3 &36.0 &50.0\\
			DarkNet23Seg~\cite{behley2019semantickitti} & &49.9 &\textbf{91.8} &74.6 &64.8 &\textbf{27.9} &84.1 &86.4 &25.5 &24.5 &32.7 &22.6 &78.3 &50.1 &64.0 &36.2 &33.6 &4.7 &55.0 &38.9 &52.2\\
			RangeNet53++~\cite{milioto2019rangenet++} & &52.2 &\textbf{91.8} &\textbf{75.2} &\textbf{65.0} &27.8 &87.4 &91.4 &25.7 &25.7 &34.4 &23.0 &80.5 &55.1 &64.6 &38.3 &38.8 &4.8 &58.6 &47.9 &\textbf{55.9}\\
			\midrule
			PointNet~\cite{qi2017pointnet} &\multirow{6}*{50K pts} &14.6 &61.6 &35.7 &15.8 &1.4 &41.4 &46.3 &0.1 &1.3 &0.3 &0.8 &31.0 &4.6 &17.6 &0.2 &0.2 &0.0 &12.9 &2.4 &3.7\\
			SPGraph~\cite{landrieu2018large} & &17.4 &45.0 &28.5 &0.6 &0.6 &64.3 &49.3 &0.1 &0.2 &0.2 &0.8 &48.9 &27.2 &24.6 &0.3 &2.7 &0.1 &20.8 &15.9 &0.8\\
			SPLATNet~\cite{su2018splatnet} & &18.4 &64.6 &39.1 &0.4 &0.0 &58.3 &58.2 &0.0 &0.0 &0.0 &0.0 &71.1 &9.9 &19.3 &0.0 &0.0 &0.0 &23.1 &5.6 &0.0\\
			PointNet++~\cite{qi2017pointnet++} & &20.1 &72.0 &41.8 &18.7 &5.6 &62.3 &53.7 &0.9 &1.9 &0.2 &0.2 &46.5 &13.8 &30.0 &0.9 &1.0 &0.0 &16.9 &6.0 &8.9\\
			TangentConv~\cite{tatarchenko2018tangent} & &40.9 &83.9 &63.9 &33.4 &15.4 &83.4 &90.8 &15.2 &2.7 &16.5 &12.1 &79.5 &49.3 &58.1 &23.0 &28.4 &8.1 &49.0 &35.8 &28.5\\
			RandLA-Net~\cite{hu2019randla} & &53.9 &90.7 &73.7 &60.3 &20.4 &86.9 &94.2 &\textbf{40.1} &26.0 &25.8 &\textbf{38.9} &81.4 &61.3 &66.8 &49.2 &48.2 &7.2 &56.3 &49.2 &47.7\\
			\midrule
			Our STPC &50K pts &\textbf{54.6} &90.8 &74.1 &63.6 &5.3 &\textbf{90.7} &\textbf{94.7} &34.4 &\textbf{48.9} &\textbf{39.7} &24.5 &\textbf{82.7} &\textbf{62.1} &\textbf{67.5} &\textbf{51.1} &\textbf{48.9} &\textbf{15.3} &\textbf{61.5} &\textbf{51.4} &47.9\\
			\bottomrule
	\end{tabular}} \vspace{-0.4cm}
\end{table}
\begin{figure*}[!t]
\begin{minipage}[t]{0.5\linewidth}
  \centering
  \includegraphics[width=6.5cm]{keys.pdf}\\
  \caption{Segmentation results with different number of dictionary atoms on the S3DIS dataset.}
  \label{change atoms number}
\end{minipage} \hspace{0.2cm}
\begin{minipage}[t]{0.5\linewidth}
  \centering
  \includegraphics[width=6.5cm]{attention0.pdf}
  \caption{The direction-dependent relationship between one neighbor point $p_k$ and each atom $\a_m$.}
  \label{attention0} \vspace{-0.2cm}
\end{minipage}
\end{figure*}
\begin{figure*}
	\begin{minipage}[b]{0.6\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{semantickitti_compare1.png}\\
		\caption{Semantic segmentation results of RandLA-Net~\cite{hu2019randla} and ours on the SemanticKITTI dataset.
		}
		\label{semantickitti_compare}
	\end{minipage}
	\begin{minipage}[b]{0.4\linewidth}
		\centering
		\tiny
		\resizebox{1\textwidth}{!}{
			\begin{tabular}{lc}
				\toprule
				methods & mIoU(\%)\\
				\midrule
				Conv. on unordered points &64.15\\
				Isotropic Conv. after max-pooling &65.31\\
				Isotropic Conv. after mean-pooling &64.17\\
				Isotropic Conv. after sum-pooling &64.87\\
				Anisotropic Conv. (ours) &66.56\\
				\bottomrule
		\end{tabular}}
				\captionof{table}{Performance with several network variants on the S3DIS dataset.}
				\label{variants}
	\end{minipage}
\end{figure*}
\vspace{-0.2cm}
\subsection{Ablation Study}\vspace{-0.2cm}
In this section, we conduct the following ablation studies for our spatial transformer point convolution module. All ablated networks are trained on area 1$\sim$4 and 6, and tested on area 5 of S3DIS dataset.

As illustrated in Figure.\ref{change atoms number}, we compare the semantic segmentation performance with different number of atoms (i.e., $M=1, 4, 9, 16, 25, 36$) in the spatial direction dictionary.
When the number of atoms increases from 1 to 25, the performance of point clouds semantic segmentation is also improving from 64.87\% to 66.56\% in term of mIoU score.
If the number of atoms continues to increase, the segmentation performance would be slightly decline.
Our STPC method achieves the best result when the number of dictionary atoms is 25 ($M=25$).
If the number of atoms is too small, the local geometric information cannot be captured finely in the anisotropic convolution filtering process.
When the number of atoms equals to 1, our method would degenerate to the sum operation, while
if the number of atoms is too large, the learnt weights tend to average on each atom, and these finer features would be lost.
Therefore, a suitable number of dictionary atoms can make the proposed SPTC method to capture the feature varieties of point clouds.

For better understanding the direction dictionary learning process,
we visualize the learned representation coefficients in  Eqn.~\ref{eq:eps8}, which can reflect the direction-dependent relationship between one neighbor point $p_k$ and each atom $\a_m \in \A$. As shown in Fig.\ref{attention0}, the vertical ordinates refer to the index of the K nearest neighbor points (i.e., $p_k \in \mathcal{N}(p_i)$, $k=1, 2, \dots, 16$) and the horizontal  ordinates denote the index of dictionary atoms (i.e., $\a_m \in \A$,  $m=1, 2, \dots, 25$).
The local neighbor points around the point $p_i$ can be projected onto different dictionary atoms.
It indicates that the feature-dense local point clouds can be sparsely encoded by sparse discrete direction dictionary.

For exploring the effectiveness of the anisotropic filtering/convolution in our network,  we evaluate the segmentation performance with several network variants on the S3DIS dataset.
Specifically, we directly perform convolution on unordered point clouds (``Conv. on unordered points"), and conduct isotropic convolution by aggregating local spatial points after max-, mean-, and sum-pooling operations, named ``Isotropic Conv. after max-pooling", ``Isotropic Conv. after mean-pooling", and ``Isotropic Conv. after sum-pooling", respectively.
As can be seen in Table~\ref{variants}, when compared with ``Conv. on unordered points", our proposed anisotropic filtering can improve the segmentation perform from 64.15\% to 66.56\%.
When compared with three different isotropic convolutions,
the performance of our proposed anisotropic convolution is better than
three isotropic convolution methods, representing improvements of 1.25\%, 2.39\% and 1.69\% over ``Isotropic Conv. after max-pooling", ``Isotropic Conv. after mean-pooling", and ``Isotropic Conv. after sum-pooling".
It proves that our anisotropic filtering method (i.e., STPC) can capture and represent implicit geometric structures, which is more robust to express those finer variances of local regions.

\section{Conclusion}
In this work, we track the point clouds semantic segmentation problem with the spatial transformer point convolution method, which can better predict the point-wise category information of 3D point clouds.
By employing the spatial direction information of point clouds, we learn the spatial direction dictionary to represent those latent geometric components.
By projecting these unordered neighbor points into the canonical dictionary space, we introduce the spatial transformer point convolution to perform the anisotropic filtering process.
The direction dictionary learning, spatial direction encoding, spatial anisotropic filtering processes can be integrated into an unified network and jointly optimized in an end-to-end fashion.
Extensive experimental results clearly demonstrate the effectiveness of the proposed spatial transformer point convolution method.

\clearpage
\textbf{\large Broader Impact}

\textbf{Future Societal Consequences.}
Recent advances in 3D hardware sensors (e.g., 3D scanners, LiDARs and RGB-D cameras) has improved the accessibility of point clouds, which can intuitively reflect the spatial geometric information of objects in 3D space.
3D semantic segmentation of point clouds, which aims to label each point with a corresponding category information, has drawn increasing attention due to its broad applicability to the fields of general environment perception.
It enables various higher-level potential applications, including autonomous driving, human-computer interaction, virtual reality and robotics.
Meanwhile, with the development of artificial intelligence, deep learning has made significantly progress in the field of computer vision.

This work "Spatial Transformer Point Convolution" focuses on the anisotropic convolution filtering method on point clouds.
Different from these previous methods, our proposed STPC can perform the anisotropic convolution filtering on point clouds by expressing those finer variances of local regions, and then further improve the performance of point clouds semantic segmentation. Extensive experiments have proved its discriminative capability.

This work will be of great significance for the applications in the 3D semantic segmentation domain and will further promote the development of artificial intelligence.

\textbf{Ethical Consideration.}
This work mainly performs scene understating in 3D space, and facilitate the
development of some high-level vision applications. For example, it can help understand our living environment to facilitate the human-computer interaction.  For the negative outcome, it will depend on the specific task and the criteria for assessing positive and negative. Besides, in this paper, all used datasets are publically available as academic research, and the evaluation metrics are also standard.

\clearpage

\small
\bibliographystyle{plain}


\section{Introduction}

Since the first discovery of a planet orbiting a main sequence star \citep{mayor_queloz_1995}, the study of exoplanets has grown tremendously. We have now begun to observe these planetary bodies in reflected light \citep[e.g.][]{chauvin2004giant}, and future large space telescopes currently under consideration are expected to make the first high-contrast reflected-light observations of rocky Earth-mass objects. In this paper, we use the NASA Planetary Spectrum Generator \citep[PSG, \url{https://psg.gsfc.nasa.gov};][]{villanueva2018planetary} to produce a grid of reflected-light spectra of rocky Earth-mass planets in Earth-like orbits around Sun-like stars with a range of simplified atmospheric abundance and aerosol profiles, with the goal of exploring and understanding the morphology of such spectra. Understanding the range of possible spectral morphologies will help future instrument designers and mission planners to efficiently plan exoplanet surveys with an eye towards characterization.

The field of exoplanetology has advanced rapidly in the past 25 years, from teasing out radial velocity detections \citep{mayor_queloz_1995} to atmospheric characterization \citep{charbonneau2002detection} and the success of sensitive transit detection surveys \citep{borucki2010kepler,ricker2014transiting,gillon2017seven}. Now, instruments such as the Very Large Telescope's Spectro-Polarimetric High-Contrast Exoplanet Research \citep[VLT-SPHERE;][]{beuzit2019sphere}, the Gemini South Telescope's Gemini Planet Imager \citep[GPI;][]{macintosh2015discovery}, the Palomar Observatory Project 1640 \citep{hinkley2011new} and the Subaru Telescope's Subaru Coronagraphic Extreme Adaptive Optics \citep[SCExAO;][]{jovanovic2015subaru} have successfully achieved direct imaging of dozens of brown dwarfs and other massive substellar objects. Future telescopes such as the Nancy Grace Roman Space Telescope \citep[RST, formerly WFIRST;][]{spergel2013wide} aim to enlarge this inventory. In addition, several large space telescope mission studies, such as the Large Ultra Violet / Optical / InfraRed Surveyor \citep[LUVOIR;][]{luvoir2019luvoir} and the Habitable Exoplanet Observatory \citep[HabEx;][]{gaudi2020habitable}, include as major mission objectives the direct imaging of Earth-sized exoplanets in Earth-like orbits around Sun-like stars in reflected (visible) light. The reflected-light spectrum of an exoplanet is a sensitive measure of the planet's atmospheric composition through photometric and spectroscopic measurements. Therefore, in advance of observations by RST and subsequent next-generation exoplanet imaging missions, recent studies \citep[such as][]{lupu2016developing,nayak2017atmospheric,batalha2018color,feng2018characterizing,smith2020detecting} have used a variety of sophisticated models to explore their potential scientific yield and develop tools to help characterize observed planets.

The use of state-of-the-art modeling techniques to produce large databases of model observations for comparison against real observations is common both within the exoplanet literature and across many fields of astronomy. For example, \citet{goyal2018library} published a model databases exploring giant planets in transmission spectroscopy, while \citet{batalha2018color} and \citet{macdonald2018exploring} computed publicly available grids of giant planet reflection spectra. \citet{allard2012models} and \citet{allard1996model} presented grids of atmospheres for low-mass stars, brown dwarfs and young giant planets. Many other databases of simulated observations exist, some going much further afield.

Future space missions will likely incorporate spatially-resolved spectroscopy (i.e. an integral field spectrograph, or IFS) to gather spectra of every source in the image plane sampled by a high-contrast exoplanet imaging instrument, and our model grid is produced at moderate resolution to accommodate future investigations of the science yield from spectroscopic analysis.  However, direct-imaging surveys of exoplanets will likely continue to utilize broadband photometric analysis methods even in the era of IFSs.  Photometric searches for Earth-like planets will most likely require an initial search without spectroscopy in order to detect background sources and efficiently characterize the system. These observations will either be photon-limited (assuming noiseless detectors) or detector-noise-limited (assuming more standard detectors), both of which benefit from photometric binning. The combination of multiple broadband photometric data points into ``colors'' has been the workhorse of astronomical classification from the early days of stellar physics up until today. Recently, some thought has been given towards applying similar methods to the characterization and study of rocky exoplanets like those in our own Solar System \citep{krissansen2016pale}, and there have been efforts to understand the efficiency of such a survey for giant planets with RST \citep{batalha2018color}.

In this study, we have applied these methods of database computation and color-color analysis towards the understanding of the science products of future exoplanet observation missions which aim to make direct observations of small, rocky exoplanets in reflected light. Acknowledging that exoplanets may vary wildly in composition, and that these variations can have dramatic consequences for observations, we first sought to parameterize and simplify the description of a rocky exoplanet atmosphere in order to reduce the dimensionality of our model grid. We then derived a physically meaningful range for each dimension, and explored the resulting parameter space with the NASA Planetary Spectrum Generator (PSG). Here, we present the resulting 141,600-model database, as well as some observations of the morphology of spectra contained therein and a simulated color-color analysis of the database.

In Section~\ref{sec:methods}, we discuss the tools and methods used in this project, including PSG. We also discuss the design of the parameter space for the grid. In Section~\ref{sec:database}, we present and validate our model database. In Section~\ref{sec:observation} we perform some basic photometric analysis using our database as a test case. Finally, in Section~\ref{sec:discussion}, we review and discuss the database and the findings from our photometric analysis, and suggest some future studies that could be performed by using this spectral database and avenues for its expansion.

\section{Methods}
\label{sec:methods}

In this section, we will discuss the techniques and choices employed in the development and calculation of our database of rocky exoplanet albedo models (hereafter ``the database''). In Subsection~\ref{subsec:design} we will describe the choices of parameters, and the chosen ranges of those parameters, that describe the parameter space explored by the database. In Subsection~\ref{subsec:PSG}, detail will be given of our choice of atmospheric forward model and its implementation in the computation of the database.

\subsection{Grid Design}
\label{subsec:design}

\begin{deluxetable*}{llll}[!]
    \tablecaption{Model Parameterization}
    \tablewidth{0pt}
    \tablehead{
    \colhead{Parameter Symbol} & \colhead{Description} & \colhead{Values Computed} & \colhead{Category}
    }
    \startdata
    C$_{abs}$ & Dominant absorbing chemical & [ H$_2$O, CH$_4$, CO$_2$, SO$_2$, O$_2$ ] & Atmosphere \\
    log(M$_{abs}$) & Dominant absorber mixing ratio & [ -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0, -0.01 ] & Atmosphere \\
    H$_2$/(H$_2$ + N$_2$) & Background Gas Hydrogen-Nitrogen Ratio & [ 0, 0.02, 0.05, 0.15, 0.5 ] & Atmosphere \\
    log(P$_0$) [bar] & Surface Atmospheric Pressure & [ -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0 ] & Atmosphere \\
    A$_s$ & Surface Albedo & [ 0.05, 0.2, 0.35, 0.5 ] & Bulk \\
    log($l$) [kg/m$^2$] & Column Cloud Mass & [ -3.0, -2.5, -2.0, -1.5, -1.0 ] & Cloud \\
    log(P$_t$) [bar] & Cloud top pressure & [ -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0] & Cloud \\
    \enddata
    \tablecomments{The parameters, and ranges of values, used to generate the primary database of 136,000 model object geometric albedo spectra presented in this paper. The full parameter space described by these parameters is explored, except for models where P$_t$ $\geq$ P$_0$. In addition, an auxiliary grid of 5,600 geometric albedo spectra were generated for model atmospheres without clouds, for a total of 141,600 models.}
    \label{table:params}
\end{deluxetable*}

The problem of describing a planet's atmosphere and predicting its observable signatures can be quite complex, and to make it tractable at this stage we have sought to find and limit ourselves to the most important planetary parameters and determine a realistic range of values for each of these parameters. For a rocky planet, we might naturally divide these parameters into two groups: atmospheric parameters, describing the gaseous envelope surrounding the planet, and bulk parameters, describing the solid rocky surface and interior of the planet. Another division is useful when describing the atmosphere, as any aerosols, clouds, dust and haze suspended in the atmosphere suggest many specialized ``cloud'' parameters. As a first simplification, we chose to model our planets in one dimension, as additional dimensions increase computation time by orders of magnitude. Further, as is common in the literature \citep[i.e.][]{feng2018characterizing,smith2020detecting, ross2019simulated,lingam2019photosynthesis}, we reduced the surface characteristics to a simple isotropic grey reflecting surface\footnote{Recent publications have studied the repercussions of this choice. See Section~\ref{sec:discussion} for more discussion on this point.}, with a given numerical albedo value A$_s$ such that 0 $\leq$ A$_s$ $\leq$ 1.

It is also convenient to express mass in terms of $g$, or surface gravity, as this term (a function of mass and surface radius) is used in calculations of the atmospheric scale height. Direct-imaging measurements are much less sensitive to scale heights than transmission spectroscopy; as such, we chose to fix this gravity parameter at Earth-standard 9.8 m/s$^2$ in order to eliminate the parameter and allow us to explore other, more impactful parameters. Further, in our  model, the planet’s radius and orbital distance become only scalar coefficients on the overall planet’s insolation; therefore, we also fix the planet's radius at 1 R$_\Earth$. Surface temperature is set equal to the atmospheric temperature at the surface (described below), while internal temperature is deemed irrelevant for reasons explained above. Thus, we have reduced the bulk parameters to one, A$_s$.

We next consider atmospheric parameters. A full and complete description of a planetary atmosphere would require a host of parameters. Atmospheric temperature varies with altitude in all known planets, as do the mixing ratios of the atmosphere’s many constituent chemical species. Describing this variation in a tractable grid would require many parameters, yet using fixed vertical profiles typically captures the main molecular morphologies dominating a spectrum. Therefore we have constructed our atmosphere to be isothermal and vertically well-mixed, such that temperature and mixing ratios do not vary with altitude. In addition, we do not consider variations in the planet’s surface or atmospheric temperature, fixing both at 300 K. Tests which varied this parameter from 200K to 500K indicate it has no appreciable effect on continuum albedo, and only a slight effect on absorption features.

This leaves a simplified atmosphere described by its chemical mixing ratios and pressure profile. With our fixed planetary surface gravity, the pressure profile can be described completely with a single parameter: the surface atmospheric pressure, P$_0$. The atmospheric pressure at any altitude can then be computed from the scale height, $H = kT/mg$. Lastly, we chose to simplify the composition of our atmosphere from potentially hundreds of constituent chemical species to only three: two of which, H$_2$ and N$_2$, describe an inert background gas; and a third which represents the dominant absorbing molecule in the atmosphere. The composition of the background gas can be described by the ratio of H$_2$ to N$_2$, which we parameterize into the fraction H$_2$/(H$_2$+N$_2$) which lies between 0 and 1. The remainder of the atmosphere requires two parameters: the choice of dominant absorbing molecule, and the mixing ratio M$_{abs}$ of that molecule in the atmosphere. Thus we have described our atmosphere with four parameters: the non-continuous choice of dominant absorber C$_{abs}$, plus continuous variables of M$_{abs}$, H$_2$/(H$_2$+N$_2$), and P$_0$.

Finally, we consider the addition of aerosols. Like the atmosphere, the addition of aerosols opens us to many free parameters. These include the aerosols’ composition, altitude, vertical extent, and particle size and density. There may be multiple layers of clouds and/or hazes at different altitudes, and they may even have different compositions. In order to maintain computational feasibility, we were forced to discard many of these possibilities. We chose a single aerosol type - a simple representative ``white'' cloud - and assigned it an effective particle radius of 10 $\mu$m. These two choices nearly eliminate any variation in the aerosol’s properties as a function of wavelength. While this assumption breaks down at longer wavelengths, white clouds provide a reasonable approximation to both water and sulfuric acid clouds at near-UV, optical, and most near-IR wavelengths. In addition, with these choices, the cloud's optical thickness, $\tau$, becomes approximately proportional to the cloud deck’s density.

\begin{figure}
    \centering
    \includegraphics[width=3.3in]{Figures/AtmoProfile.png}
    \caption{An example of our atmospheric profile, showing a .1\% water vapor atmosphere with no H$_2$ and Earth-like surface pressure. Clouds profiles with top pressures at 10$^{-0.5}$ bars and 10$^{-2.0}$ bars are shown, each with different abundances (column masses of 10$^{-2.5}$ and 10$^{-1.5}$ kg/m$^2$ respectively). The N$_2$ and H$_2$O abundances are volume mixing ratios, while the cloud abundance is measured in mass mixing ratio.}
    \label{fig:atmoprof}
\end{figure}

As a further simplifying choice, we constructed a standard cloud profile to be placed in our atmospheres. This profile, whose vertical location is defined by the pressure level P$_t$ at the top of the cloud, consists primarily of several pressure layers of ``thick'' clouds. Above the cloud top, the cloud particle mass mixing ratio drops rapidly to 10$^{-10}$ over 3.33 scale heights. Beyond this point, no cloud particles are included in the atmosphere. Below the thick cloud ``core'', the mass mixing ratio likewise drops off precipitously to 10$^{-9}$ over 3 scale heights, where it remains down to the planet's surface; see Figure~\ref{fig:atmoprof} for an example atmosphere profile. The mass mixing ratio at any given level is designed such that the sum of cloud mass across all pressure levels is equal to a chosen cloud column mass (hereafter $l$) in units of kg/m$^2$, and yet the cloud is localized at a specific pressure level. Thus, we have reduced our cloud description to two parameters: P$_t$ and $l$.

Through the method described above, we have reduced our rocky exoplanet models from potentially dozens of dimensions to seven. In the process, however, we have discarded some potentially valuable parameters. We plan to gradually restore some of these lost dimensions in future work. Additionally, we recognize that some of the parameters we chose not to explore (such as planet mass) may be determined more easily through other means (such as radial velocity or astrometry measurements).

We expect that the seven parameters (described in Table~\ref{table:params}) adequately describe a rocky Earth-like exoplanet to a first approximation. The combination of these parameters, and the values (shown in Table~\ref{table:params}) chosen for them, gave us a parameter space consisting of 136,000 models. To supplement this, we computed an auxiliary aerosol-free grid, in which all aerosol-related parameters were ignored; thus, noting that we have two ``cloud'' parameters, this auxiliary grid is five-dimensional, containing 5,600 models. The final database thus contains 141,600 model geometric albedo spectra.

\subsection{PSG}
\label{subsec:PSG}

The individual model albedo spectra of our grid were generated using the NASA Planetary Spectrum Generator \citep[PSG;][]{villanueva2018planetary}, an online browser-accessible radiative transfer suite. PSG, developed at the NASA Goddard Space Flight Center, is a comprehensive all-in-one modeling package capable of accurately simulating observations of planetary and sub-planetary-mass objects. It is capable of employing a variety of different molecular line lists including collision-induced absorption cross-sections, as well as pre-computed aerosol scattering models and surface material reflectances, and employs data from sources such as MERRA-2 and MCD to further improve models of Solar System objects. In addition to modeling spectra, PSG also includes software for simulating observations, including observatory throughput and detector noise, and calculates observational uncertainties.

The online, publicly available browser-accessible PSG application is frequently updated with improved functionality. For this project, we used the installable Dockerized version of the PSG application program interface (API); for reference, the installation was performed on 12/13/19. For detailed information about PSG, we recommend \citet{villanueva2018planetary}, or the documentation available on the PSG website (psg.gsfc.nasa.gov). We installed and ran multiple instances of PSG on the Goddard Private Cloud (GPC) at NASA Goddard, and utilized the GPC Job Engine\footnote{https://gitlab.com/mdmoore25404/gpc-je} to distribute calculations across $\sim100$ virtual machines dedicated to the project.

When performing scattering calculations, as required by this work, convergence of the radiative transfer calculations in situations of high optical depth $\tau$ can be challenging. As our model database includes an opacity source that is in some cases extremely dense (clouds), this issue required some special care. For scattering calculations, PSG internally restricts the opacity per layer to be moderate levels ($\tau<1$) in order to ensure convergence of the scattering calculations. Therefore, to avoid this from becoming an issue, we pre-computed the optical depth at each layer within the core of the cloud deck. Any individual layer which we found to have a $\tau \geq1$ was divided into two layers, thus reducing the opacity of each individual layer. Because PSG performs computations including the spherical geometry of the planet's atmosphere, the initial optical depth evaluation was performed near the limb of the planet where the path length of a photon through each layer is greatest.

\section{Results}
\label{sec:database}

\begin{figure*}[hbtp]
    \centering

    \includegraphics[width=7.0in]{Figures/DatabaseSurvey.png}
    \caption{A survey of the variable parameters in our database. The \textbf{bold} line represents the same model, our fiducial model with an H$_2$O/N$_2$ atmosphere, in each sub-panel. The sub-panels show explorations of the dominant absorbing chemical C$_{abs}$ (offset for clarity), its mixing ratio log(M$_{abs})$, H$_2$ content in the background gas, surface pressure log(P$_0$), surface albedo A$_S$, cloud density parameter log($l$), and cloud top altitude P$_t$. A normalized stellar spectrum from a G2 star and the albedo spectrum of a 1-bar N$_2$ atmosphere are included in the bottom right for reference.  See Section~\ref{sec:database} for more details.}
    \label{fig:Database}
\end{figure*}

The full database of geometric albedo spectra for 136,000 cloudy and 5,600 cloud-free model rocky exoplanets is available for download\footnote{ http://doi.org/10.5281/zenodo.3743500}. In addition, the database can used to produce models across the full parameter range of the grid using the REPAST online tool at the Exoplanet Modeling and Analysis Center (EMAC) website\footnote{http://emac.gsfc.nasa.gov/repast}.

In this Section we will explore some features of the albedo spectra contained within the database. To better highlight the differences between models and to aid in comparisons, we adopt a set of ``fiducial'' parameters. The albedo spectra associated with our fiducial model is present in all figures in this section. This fiducial model has an H$_2$O mixing ratio of $10^{-3}$, and a background gas of pure N$_2$ (such that the parameter H$_2$/(N$_2$+H$_2$) = 0). The surface pressure is set at 1 bar, and the surface albedo is 0.05. We choose a column cloud mass of $10^{-2.5}$ g/kg$^2$, with the cloud top at $10^{-0.5}$ bars.

In Figure~\ref{fig:Database} we plot a sample of the model spectra when varying each of our model parameters.  At the top left, we compare the full spectra for our various absorbing species. It is clear that, at this mixing ratio, CO$_2$ and, beyond 0.3 $\mu$m, O$_2$ have negligible impacts on the spectrum, while many of the H$_2$O and CH$_4$ features in the near infrared region overlap. SO$_2$ has a well-known very strong  absorption feature to the blue of 0.4 $\mu$m. For clarity, these spectra have been plotted to include an artificial offset from each other.

When plotting the absorber mixing ratio for H$_2$O, we see an expected progression from no absorption for extremely low mixing ratios, to very strong absorption for extremely high mixing ratios. Of interest are the increasing number of detectable features in the optical at high mixing ratios, as well as the impact that the changing mean molecular weight has on the effective slope of the spectrum at short wavelengths; this region is dominated by Rayleigh scattering at low absorber concentrations, but drops faster as more short-wavelength flux is absorbed.  Similar distortions of the overall spectral shape are present for the other absorbers.

We find that surface pressure has a significant impact on spectral shapes, as the Rayleigh slope changes dramatically with surface pressure. In addition, high atmosphere pressure has a similar effect on molecular spectral features as increased mixing ratios. The background gas composition has a relatively minor impact at low surface pressure; even when the background gas is 50\% H$_2$, the variation in the spectra is almost undetectable. However, at higher surface pressures than that chosen for our fiducial model, we do note that the effects of H$_2$-H$_2$ collision-induced absorption come into play, creating features beyond 1 $\mu$m.  The three parameters affecting the surface and the aerosol properties (A$_S$, P$_t$ and $l$) all affect the overall continuum shape of the spectra.  Raising the bulk reflectivity of the surface of the planet clearly raises continuum beyond the Rayleigh slope, while having significantly lesser effects on the shape of the Rayleigh slope itself or on the bottom levels of molecular absorption features. The impact of cloud density on the continuum albedo is very similar, but there is a greater effect on the Rayleigh slope and a slightly lesser effect on molecular absorption features in comparison. At high surface pressures and high altitude clouds, the tendency of denser clouds to suppress Rayleigh effects by preventing access to the lower portions of the atmospheres competes with their tendency to raise the overall albedo at all wavelengths.  The impact of cloud top pressure/altitude on the continuum is subtle, with a slight loss of the Rayleigh slope for high-altitude clouds, but the greater impact is on molecular features. High-altitude clouds mask molecular absorption features, as a smaller amount of the atmosphere is accessible. This effect is most noticeable in high-pressure atmospheres.

In comparison to the albedo spectra, we also include a stellar spectrum for a G2 star in Figure~\ref{fig:Database}, normalized to the peak of the stellar output, in order to provide context for the impact of the stellar spectral energy distribution on the planet's observed spectrum (star $*$ albedo). Due to the increased flux at the peak of the stellar spectrum, planetary flux and absorption features in this region will be much easier to detect than features in lower-flux regions. For example, the high absorption of SO$_2$ lies in a region of very low solar output; thus this absorption feature may be difficult to detect due to large observational uncertainties. We also plot a representative albedo spectrum, depicting a 1-bar N$_2$ atmosphere over a 0.5 A$_s$ surface with no clouds or absorbing molecules.

\subsection{Comparison with Solar System Rocky Bodies}
\label{subsec:solarsystem}

\begin{figure*}[hbtp]
    \centering

    \includegraphics[width=5.7in]{Figures/Fourplot.png}
    \caption{This figure compares each of the atmosphere-rich rocky bodies of the Solar System to one or more analog models from our database for each planet, highlighting the differences between the original, vetted PSG simulation of the planet and our simplified model. For Earth, Mars, and Titan, we also include a transitional model which includes the most significant simplifications but retains much of the complexity of the original model. The parameters for the database analogs are listed in Table~\ref{table:models}.}
    \label{fig:fourplot}
\end{figure*}

\begin{deluxetable*}{cclcccl}[!]
    \tablecaption{Database Analogs to Solar System Planets}
    \tablewidth{0pt}
    \tablehead{
    \colhead{Analog} & \colhead{Associated} & \colhead{} & \colhead{} & \colhead{} & \colhead{} & \colhead{} \\
    \colhead{Name} & \colhead{S.S. Planet} & \colhead{C$_{abs}$} & \colhead{M$_{abs}$} & \colhead{P$_0$} & \colhead{A$_s$} & \colhead{Cloud Param.}
    }
    \startdata
    M1 & Mars & CO$_2$ & -0.01 & -1.0 & 0.2 & None \\
    E1 & Earth & H$_2$O & -2.0 & 0.0 &  0.35 & None \\
    V1 & Venus & CO$_2$ & -0.01 & 2.0 &  0.2 & log($l$)=-1.5, log(P$_t$)=0.0  \\
    V2 & Venus & SO$_2$ & -5.0 & 2.0 &  0.2 & log($l$)=-1.5, log(P$_t$)=0.0 \\
    T1 & Titan & CH$_4$ & -2.0 & 0.0 &  0.2 & log($l$)=-2.0, log(P$_t$)=-1.0 \\
    \enddata
    \centering
    \tablecomments{All database analogs have no H$_2$.}
    \label{table:models}
\end{deluxetable*}

Since our spectral database examines a range of atmospheric properties for an Earth-mass planet, it is a natural question to ask how these model predictions compare with the albedo spectra of the terrestrial-mass bodies with significant atmospheres in our own Solar System - Venus, Earth, Mars and Titan

However, since our model grid is highly simplified in terms of the atmospheric composition, cloud composition and structure, and surface albedo, any attempt to match the spectra of Solar System planets to the outputs for a specific region of parameter space in the current database is not justified.  In particular, the Solar System terrestrial bodies have a mix of dominant absorbing species: the Venusian atmosphere includes both SO$_2$ and CO$_2$, Titan's spectrum includes absorption from many hydrocarbon species, and Earth's atmosphere obviously includes a number of atmospheric constituents with significant absorption signatures.  Similarly, the cloud properties of Venus and Titan are complex, with combinations of photochemical hazes contributing significantly to the optical spectrum.  Finally, the surface albedo spectrum of Mars has a significant wavelength-dependent reflectivity across the visible spectral region that produces the dominant color contribution to the reflected-light spectrum.  We expressly choose to avoid the added complexities of examining a myriad of specific cloud and haze species and composition-specific surface albedo reflectivities in this study, in order to focus on a small number of variables that are agnostic to the particulars of atmosphere and surface chemistry; we leave an examination of these additional factors for future work, as discussed in Section~\ref{sec:discussion}.

Instead, in order to identify where the atmospheric abundances, surface pressure, cloud properties and average surface albedos for the four atmosphere-rich Solar System rocky bodies would fall within our grid, we have defined analogs within our spectral database, choosing values for these four parameters that most closely approximate the Solar System values. Figure~\ref{fig:fourplot} depicts the spectral differences between the real Solar System bodies and our grid analogs.  It is clear that for Mars and Titan, the dominant change is the drop in flux at short wavelengths; for Mars this is due to the surface reflectivity, while for Titan this is due to the reflectivity of hydrocarbon hazes.  For Earth, our database analog primarily diverges in the near-UV due to absorption from O$_3$; since our grid only has a single absorber in each spectrum, only spectral features of H$_2$O are present.

Finally, for Venus we define two database analogs, since both SO$_2$ and CO$_2$ produce significant absorption across the Venus spectrum. Interestingly, Figure~\ref{fig:fourplot} still shows a difference in the morphology of the SO$_2$ band between the realistic Venus spectrum and our database analog; this is due to the impact of scattering on incoming short-wavelength radiation and the low altitude of SO$_2$ in Venus' atmosphere, which diminishes the absorption.

In Table~\ref{table:models} we list the properties for each of the grid analogs to the Solar System bodies, and the associated label names that we include in subsequent figures. In Section~\ref{sss:analogSS} we will explore what impact of the shift from detailed Solar System planet model to our database analog on our ability to differentiate between them using photometric colors.

\section{Photometric Color Analysis}
\label{sec:observation}

In this Section we use our spectral database to explore the potential to differentiate the properties of rocky planet atmospheres from a multi-band photometric survey. This strategy is motivated by the work of \citet{krissansen2016pale}, hereafter KT16, who showed that producing integrated bandpass filter measurements from optical spectra of Solar System bodies could distinguish potentially Earth-like exoplanets from planets similar to the other (non-habitable) rocky planets in our Solar System. They determined that bandpass ranges of 431-531 nm (``blue''), 569-693 nm (``green''), and 770-894 nm (``red'') could be combined to create sufficient separation in color-color space between Earth and many Solar System bodies.

In Section~\ref{subsec:bandpass} we examine two different sets of photometric bandpasses based on the current expectations for the LUVOIR mission concept, as well as the bandpasses from KT16. We use color-color plots as case studies in Section~\ref{subsec:plots} to consider the efficiency of a purely photometric survey for distinguishing planetary atmospheric characteristics. Section~\ref{subsec:statistics} will review the statistics employed in this analysis.

We note that all of our photometric colors are actually ``albedo colors''- we do not include the impact of the stellar spectrum in the color-color plots.  This is to allow a clear comparison to the example spectra shown in Figure~\ref{fig:Database}. Multiplying a stellar spectrum to the albedo spectra would only apply a multiplication factor to every value, but would not change the statistical analysis of filter positions since they are determined through correlations between colors and various parameters.  The impact of the stellar spectrum is taken into account in Section~\ref{sec:observation} when calculating the actual observable flux for each filter.

\subsection{LUVOIR Photometric Bandpass Case Studies}
\label{subsec:bandpass}

The LUVOIR space observatory concept is one of the major exoplanet direct imaging missions under study for future development; due to the large aperture size, LUVOIR would search hundreds of stars for Earth-like planets as part of an initial reconnaissance survey \citep{luvoir2019luvoir}. We therefore chose to use assumptions about the expected LUVOIR bandpasses for two of our three case studies for a photometric survey. While many features of such a survey are yet to be decided, it is expected that the first step would be the discovery of planets by imaging at a wavelength near the peak photon flux of the planet. As such, we examined an Earth spectrum (generated by PSG using MERRA-2 data) to find this peak flux in a filter bandpass. We found this to be at approximately 0.58 $\mu$m. As the LUVOIR concept would be able to observe in ultraviolet, optical, and a red or near-infrared channel simultaneously, we paired this with a second photometric point at 0.4 $\mu$m. We justify this wavelength choice in Section~\ref{subsec:statistics}. Finally, we paired these two with a third filter choice at an indeterminate point red of 0.58 $\mu$m. For Case 1, we assumed the photometric bandpasses of the filters would be 10\%, based on the expected achievable bandpass width for an on-axis telescope design baselined for the 15-meter LUVOIR-A concept \citep{luvoir2019luvoir}.

The second case study is identical to the first, except that the width of the filter bandpass has been changed to 20\%. This bandpass would be achievable with an off-axis design, as baselined for the smaller LUVOIR-B concept.

\subsection{Statistical Methodology}
\label{subsec:statistics}

For the LUVOIR bandpasses, we employed the Pearson correlation coefficient and Grubbs's test for outliers to guide our decisions of bandpasses to explore in color-color space.

The Pearson correlation coefficient r is a measure of how close the data points lie to a linear correlation; $-1\leq$ r $\leq 1$ where r $=\pm 1$ indicates a perfectly (positive or negative) linear correlation and r$=0$ indicates no correlation. The Pearson correlation coefficient r$_{xy}$ between variables x and y can be calculated for a sample by

\begin{equation}
    r_{xy} = \frac{\sum_{i=1}^{n}{(x_i-\bar{x})(y_i-\bar{y})}}{\sqrt{\sum_{i=1}^{n}{(x_i-\bar{x})^2}}\sqrt{\sum_{i=1}^{n}(y_i-\bar{y})^2}}
\end{equation}

where n is the sample size, x$_i$, y$_i$ are individual sample points, and $\bar x$, $\bar y$ are the sample means for x and y respectively. Interpretation of the Pearson coefficients requires some caution as they do not capture such information as the `slope' of the correlation. However, Pearson correlation coefficients are useful for gaining some intuition for the value of each filter.

Grubbs's test for outliers, as the name suggests, is useful for locating unusual data points in a sample. The one-sided Grubbs test statistic G is defined as

\begin{equation}
    G=\frac{y_{max}-\bar y}{s}
\end{equation}
\nopagebreak
where s is the sample standard deviation, y$_{max}$ is the largest value in the sample, and $\bar y$ is defined as above. Then, if

\begin{equation}
    G>\frac{N-1}{\sqrt{N}} \sqrt{\frac{t_{\alpha/N,N-2}^2}{N-2+t_{\alpha/N,N-2}^2}}
\end{equation}

where N is the sample size and $t_{\alpha/N,N-2}^2$ is the upper critical value of a t-distribution, we reject the hypothesis that y$_{max}$ is not an outlier with a significance level $\alpha$. We then remove $y_{max}$ from the sample and repeat the test until no outliers are found.

To determine the best filter positions to distinguish the impact of each specific parameter, we calculated the Pearson correlation coefficients for the color of different filters combined with our 0.58 $\mu$m central filter, and then applied Grubb's test in order to help find ``abnormally strong'' correlations. Like the Pearson correlation coefficient, the Grubbs test should only be regarded as a guide. It was used only to highlight potentially useful combinations for further study.

\subsection{Survey Results}
\label{subsec:plots}

In the following section, we will review the results of the three cases laid out in Section~\ref{subsec:bandpass}. In each of the figures below, we highlight five models in our grid chosen to approximate the atmosphere-rich terrestrial bodies within our Solar System. Details of these five models are given in Table~\ref{table:models} and Section~\ref{sss:analogSS}.

\begin{figure}[ht!]

    \includegraphics[width=3.3in]{Figures/PCCPlotFinal033020.png}
    \caption{The absolute values of Pearson correlation coefficients calculated for individual 10\% filters in the UV band with a variety of central wavelengths. These measure how well a change in the parameter value correlates with a change in the filter response. Values close to 0 indicate very little or no correlation, while values close to $\pm 1$ indicate very strong, almost linear correlations. Notice the shift in correlation emphasis from SO$_2$ concentration to atmospheric and cloud parameters around 0.4 $\mu$m. We chose this location for the center of our UV bandpass as mentioned in Section~\ref{subsec:bandpass}.}
    \label{fig:PearsonUV}
\end{figure}

\begin{figure}[t!]

    \includegraphics[width=3.3in]{Figures/LUVOIRchems10.png}
    \caption{Our model database, and several Solar System planets, plotted on our chosen axes. Each axis is the logarithm of the ratio of two 10\% filters, centered at the indicated wavelengths. Here we have highlighted the dominant absorbing chemical in each model on the plot. We note the four distinct branches separating SO$_2$, O$_2$, H$_2$O, and CH$_4$. On these axes, CO$_2$-dominated atmospheres fail to manifest a distinct branch due to a lack of spectral features in any of the bands.}
    \label{fig:case1}
\end{figure}

\subsubsection{LUVOIR Bandpasses: Filter Position Selection}
\label{sss:LUVOIRfilterchoices}

For our UV filter options, the only parameters that showed significant Pearson correlation coefficients were SO$_2$ (due to the strong absorption feature) and surface pressure (due to Rayleigh scattering); see Figure~\ref{fig:PearsonUV}) for a visual representation of the Pearson coefficients. In order to be sensitive to the impact of both parameters, we decided to choose a filter position at 0.4 $\mu$m for both our 10\% and 20\% cases, since it covers a region of reasonable sensitivity to both parameters.

For the 10\% filter case, we found that a third filter located at 0.80 $\mu$m, in combination with the 0.4 and 0.58 $\mu$m filters, produced an outlier-level correlation for four of the five chosen molecular species concentrations. We note, however, that this correlation was still quite low (approximately 0.2-0.25). For our 20\% filter scenario we found that a similarly located filter to that chosen in the 10\% scenario, at 0.78$\mu$m, yielded similar - though diluted - results.

In addition to finding outliers in our 0.8 and 0.78 $\mu$m filters, which we ultimately chose as our final selections, we also noted several other filters which correlated well with several parameters. For example, blue and green optical filters (at 0.58, 0.66, and 0.68 $\mu$m) correlated strongly with cloud density when combined with our pre-chosen filters, while some filters in the near-infrared showed correlation with O$_2$ as well as CH$_4$ (at 1.34 and 1.36 $\mu$m) and H$_2$O (at several locations between 1.38 and 1.46 $\mu$m). However, we did not choose these as they were only able to distinguish two or, in some cases, three of the five molecules. We did not find any strong or outlier correlations for CO$_2$ in our study. In each case, we plotted the filter combinations in color space in order to see how distinct (or degenerate) the correlations between different parameters was. Our final choice of the 0.8 and 0.78 $\mu$m filters for 10\% and 20\% respectively was made based on achieving the maximum color-color correlation for the most parameters.

\subsubsection{LUVOIR Bandpasses: Color-Color Plots}
\label{sss:cases1n2}

\begin{figure}[!]

    \includegraphics[width=3.3in]{Figures/LUVOIRmixes10.png}
    \caption{As Figure~\ref{fig:case1}, however in this plot we have highlighted the mixing ratio of the primary absorber. Notice that the H$_2$O and O$_2$ mixing ratio must exceed 0.1, and CH$_4$ mixing ratio must exceed 0.01, before the models move away from the central grouping of models. This suggests that photometric characterization will only be effective for identifying atmospheres with a high concentration of molecular absorbers.}
    \label{fig:mixrats}
\end{figure}

\begin{figure}[t!]
    \includegraphics[width=3.3in]{Figures/LUVOIRchems20.png}
    \caption{As Figure~\ref{fig:case1}, however in this plot we use 20\% bandpass filters. The results are very similar to Figure~\ref{fig:case1}. However, note that H$_2$O and CH$_4$ are more mixed and harder to distinguish in this case. The difficulties observed in Figure~\ref{fig:mixrats} are also present here.}
    \label{fig:case2}
\end{figure}

\begin{figure}
    \includegraphics[width=3.3in]{Figures/transition.png}
    \caption{Here, we highlight the shift in color-color space from each of the terrestrial solar system planets to the analogs referenced in the paper. As shown, the differences are extremely minor for the relatively simple Earth and Venus models. Titan and Mars shift in a more significant way, as shown in Figure~\ref{fig:fourplot} and discussed in Section~\ref{sec:discussion}. The axes here are the same as those used in the 10\% bandpass case. The numbers adjacent to each point refer to the numbered models in Figure~\ref{fig:fourplot}. }
    \label{fig:transit}
\end{figure}

Figure~\ref{fig:case1} shows a color-color plot of our grid using our three chosen 10\% filters, as well as our Solar System database analogs. Note that each of H$_2$O, CH$_4$, CO$_2$, and SO$_2$ appear to follow their own branch on this plot. However, the separation from the central ``group'' of data points does not manifest until a significant concentration of the chemical is present in the atmosphere, as indicated in Figure~\ref{fig:mixrats}. In addition, we observe that our Solar System analog objects all lie in the central region where the branches overlap; the only analog S.S. planet with any appreciable separation is E1, which is due to the absorption from H$_2$O.

Figure~\ref{fig:case2} shows a color-color plot for our chosen 20\% filters. This set shows similar trends as those observed in Figure~\ref{fig:case1}; most notably, there are four branches corresponding to H$_2$O, CH$_4$, CO$_2$, and SO$_2$. However, there is a greater degree of degeneracy between H$_2$O and CH$_4$ on these axes, and the same clustering of solar system analog objects.

\subsubsection{Exploration of Solar System Analogs in Color-Color Space for 10\% Bandpasses}
\label{sss:analogSS}

In Figures \ref{fig:case1}, \ref{fig:mixrats}, and \ref{fig:case2}, it is clear that our Solar System database analogs are tightly clustered and difficult to differentiate relative to other models covering a variety of atmospheric properties. It is instructive to explore this phenomenon and how it compares against the positions of real (non-analog) Solar System objects on these axes.

Figure~\ref{fig:transit} traces the movement in color-color space for each Solar System planet as the model parameters are adjusted to match the chosen analog models from our database; the numbers proceed from the realistic models of Earth, Mars, Titan, and Venus to our E1, M1, T1, and V1 and V2 analog objects, and are matched to the numbers of the spectra presented in Figure~\ref{fig:fourplot}. It is clear that Earth and Venus have their positions shifted only slightly by the transition, since both planets' spectra are dominated by clouds that are featureless across the visible region of the spectrum, and the key wavelength-dependent variation is driven by atmospheric absorption.

In contrast, the positions of Mars and Titan are changed significantly compared with our database analogs.  For Mars, this is due to the removal of the detailed surface reflectivity present in the realistic Mars spectral model, while the presence of tholin hazes significantly affect the spectral color of Titan.  As stated earlier, we chose not to include either of these factors in the current database to preserve simplicity and a minimal parameter set; we leave this work to future studies.

\subsubsection{Bandpasses from Krissansen-Totton et al.}
\label{sss:case3}

\begin{figure}[t!]
    \includegraphics[width=3.3in]{Figures/KT16chems.png}
    \caption{Our database is plotted on the same axes as used in Figure 2 in from \citet{krissansen2016pale}. Note the strong clustering of our Solar System analog models on these axes, as well as the lack of a distinct SO$_2$ branch.}
    \label{fig:KT16}
\end{figure}

Figure~\ref{fig:KT16} shows our database plotted on the red-blue-green axes chosen by KT16 as described in Section~\ref{subsec:solarsystem}. It is notable that our Solar System analogs models are very tightly clustered on these axes, even more-so than in Figure~\ref{fig:case1} and Figure~\ref{fig:case2}, except now the Mars analog is slightly separated.  We also note that the KT16 bandpasses reduce the overlap between the methane and water ``branches'' at lower concentrations. However, the SO$_2$ branch is completely eliminated, and the O$_2$ branch is shorter in both axes.

\subsubsection{Estimating Observing Time for Discriminating Atmospheric Properties of Rocky Exoplanets}
\label{sec:obstime}

\begin{figure*}[hbtp]
    \centering
    \includegraphics[width=6.0in]{Figures/ChemsWithSolar10.png}
    \caption{In this plot, a G2 stellar spectrum has been applied to the raw albedo spectra prior to calculating spectroscopic data from 10\% filters, in order to approximate actual planetary fluxes for a planet around a solar-type star. Further, the X- and Y-axes have been changed to show the absolute difference of the filter response ratio from that of a cloud-free, pure N2 atmosphere at 1 bar P$_0$. Included are the 3-$\sigma$ noise floors for a planet at 10 parsec for 1, 3, 10, and 100 hour observations with a 15m telescope, calculated with PSG (assuming zero detector noise), which indicate how difficult it would be to distinguish between a given model and a simple world with no absorbers. Letters indicate the position of the Solar System planet analogs described in Section~\ref{subsec:solarsystem}.}
    \label{fig:solarplot}
\end{figure*}

In order to assess the effectiveness with which different planetary atmosphere compositions can be distinguished for a specific observatory and instrument suite, it is necessary to actually calculate the measurement uncertainty of observations in the chosen photometric bands, and compare that to the expected difference in photometric color for different models with different atmospheric compositions.

In Figure~\ref{fig:solarplot}, we have multiplied a G2 stellar spectrum with our raw model albedo spectra to determine real photometric colors in the LUVOIR 10\% bandpasses, and then calculated the difference in color of that model compared with a cloud-free, pure N2 atmosphere with 1-bar surface pressure; this allows us to evaluate how well a certain model could be distinguished from an absorber-free atmosphere.  We also calculated the photon-limited uncertainty in the determination of each color for an observation of a solar-type star at 10 parsec with a 15-meter telescope at a variety of integration times using PSG and the telescope and instrument parameters provided in the LUVOIR study report \citep{luvoir2019luvoir}. We chose to use a ``perfect'' detector for this calculation, in order to maintain applicability to multiple future observing projects; thus, no detector noise, read noise, dark current, or other sources of uncertainty are considered here. As such, the noise floors depicted are likely optimistic.

Under these conditions, we found that a 3-hour minimum observation would be required to begin distinguishing the H$_2$O and CH$_4$ models from a pure-N$_2$ atmosphere. With a 10-hour observation, atmospheres with high fractions of O$_2$ and SO$_2$ could be identified.  Finally, with a 100-hour observation, a planet with Earth's H$_2$O abundance could be positively identified.

\section{Discussion}
\label{sec:discussion}

We have computed a database of 141,600 geometric albedo spectra for rocky Earth-mass exoplanets. Individual planet spectra can be derived from the database using the online calculator hosted at the EMAC website\footnote{http://emac.gsfc.nasa.gov/repast}, and the whole database is publicly available\footnote{http://doi.org/10.5281/zenodo.3743500}.

Using this database, we simulated a photometric color-color analysis using a variety of filter combinations. We examined photometric filters with 10\% and 20\% bandpasses in the ultraviolet, visible, and near-infrared. We found that, given a plausible discovery bandpass at 0.58 $\mu$m and a secondary ultraviolet bandpass at 0.4 $\mu$m, adding a third filter bandpass at 0.8 $\mu$m (for 10\% bandpasses) or 0.78 $\mu$m (for 20\% bandpasses) allowed discrimination between four of the five molecular absorbers simulated in this study. We note however that exceptional concentrations of the chemical species we examined are required to distinguish planets using shorter integration times; in addition, we did not include all sources of instrument noise, and therefore the results could be even less promising than described here.

As mentioned in Section~\ref{subsec:design}, we made a number of other simplifying assumptions in the interest of computational feasibility. We intend to address several of these assumptions in future expansions on this database, beginning with the choice of a gray surface albedo. While this simplifying assumption is common in the literature and is reasonable as a first approximation, several studies \citep{madden2020surfaces, o2018vegetation, schwieterman2015nonphotosynthetic, kaltenegger2007spectral} have shown that the surface composition of rocky exoplanets can have a significant effect on reflected light spectra. The impact of using a grey surface albedo is also quite notable in comparison to objects in our Solar System; see for example Figures \ref{fig:fourplot} and \ref{fig:transit}. However, the addition of appropriate surfaces would correct this discrepancy, and allow a model database to more appropriately cover the regimes that Solar System objects inhabit. Because of this, and other comparisons (e.g. to Earth), expanding our database to cover a variety of surfaces is a very high priority for the further development of our models.

Our list of dominant absorbers did not include molecular species that are produced primarily by photochemistry in cool rocky planet atmospheres. In particular, O$_3$ is a byproduct of O$_2$ photochemistry, and is responsible for a significant spectral absorption feature in Earth's spectrum.  Similarly, we did not include CO, or various chemical and photochemical byproducts of CH$_4$, which are present in the atmospheres of Venus and Titan. We decided not to include these spectral absorbers since their abundance is highly dependent on a planet's bulk atmospheric composition and the spectral energy distribution of a planet's host star, and therefore the parameter space where these absorbers would be important will be restricted. For this study, we were focused on bulk constituents that are expected to be present in the majority of cool rocky planet atmospheres. We intend to explore ways to include more realistic chemistry into future spectral databases.

We included molecular hydrogen H$_2$ as a potential background gas component in this database, since it has been suggested \citep[][e.g. in]{ramirez2014warming} that Mars retained an atmosphere with significant amounts of H2 briefly after formation, with climate implications \citep[see also][for implications for close-orbiting planets.]{koll2019hot}. However, the detection of H$_2$ in terrestrial planet atmospheres primarily relies on H$_2$-H$_2$ collision induced absorption, and our results show that H$_2$-H$_2$ CIA absorption does not become apparent until concentrations and/or pressures become fairly extreme. Thus, this parameter may be restricted in future database releases in favor of other, more impactful dimensions.

Additional areas of future work include variations in aerosols; our study was limited to a featureless cloud by computational necessity. However, the addition of high-altitude photochemical hazes are also worthy directions to explore, as these are common in rocky planet atmospheres within the Solar System and have an impact comparable to that of planetary surfaces. An exploration of other particulate matter (such as dust) suspended in the atmosphere is also under consideration. Also, potential follow-ups to the photometric analysis presented in this paper include the application of machine learning tools \citep[as in][]{batalha2018color} to make the filter choices more agnostic.

\acknowledgements
    The authors would like to acknowledge Y. K. Feng and Dr. Joshua Krissansen-Totton for providing albedo spectra for validation; Dr. Natasha Batalha for guiding us towards useful analysis techniques; and Dr. Chris Stark and Dr. Matthew Bolcar for descriptions of expected LUVOIR spectroscopic approaches. The material is based upon work supported by NASA under award number 80GSFC17M0002.

\software
    SciPy \citep{virtanen2019scipy}, NumPy \citep{van2011numpy}, matplotlib \citep{hunter2007matplotlib}, pandas \citep{mckinney2011pandas}, astropy \citep{astropy:2013,astropy:2018}, PSG \citep{villanueva2018planetary}