# Performance comparison of 3D correspondence grouping algorithm for 3D plant
point clouds

Shiva Azimi
Department of Electrical Engineering
Indian institute of Technology Delhi
New Delhi, India 110016
<EMAIL_ADDRESS>
Tapan K. Gandhi
Department of Electrical Engineering
Indian institute of Technology Delhi
New Delhi, India 110016
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

Plant Phenomics can be used to monitor the health and the growth of plants.
Computer vision applications like stereo reconstruction, image retrieval,
object tracking, and object recognition play an important role in imaging
based plant phenotyping. This paper offers a comparative evaluation of some
popular 3D correspondence grouping algorithms, motivated by the important role
that they can play in tasks such as model creation, plant recognition and
identifying plant parts. Another contribution of this paper is the extension
of 2D maximum likelihood matching to 3D Maximum Likelihood Estimation Sample
Consensus (MLEASAC). MLESAC is efficient and is computationally less intense
than 3D random sample consensus (RANSAC). We test these algorithms on 3D point
clouds of plants along with two standard benchmarks addressing shape retrieval
and point cloud registration scenarios. The performance is evaluated in terms
of precision and recall.

_K_ eywords Computer Vision $\cdot$ Correspondence Grouping $\cdot$ Plant
Phenotyping $\cdot$ 3D Point Cloud

##  1 Introduction

Providing food to increasing world population has become a challenging task
due to loss of arable land, climatic changes and many man made factors. This
requirement for increasing the crop yield has necessitated the development of
modern tools and techniques for increasing crop yield [1]. Plant phenotyping
deals with the measurement of changes in phenomes in reaction to genetic and
environmental changes and is an important plant science technique that can be
used to monitor plant health and growth.

Recently, image-based high-throughput plant phenotyping has become an
important computer vision research problem [2]. Measurements of visual plant
features can be done using either two dimensional (2D) [3], [4] or three
dimensional (3D) models [5] of the plant. However, since 3D models can capture
much more information regarding the plant physiology as compared to 2D models,
3D models are preferred for plant phenotyping.

Finding correct correspondences between 3D shapes is a cornerstone in 3D
computer vision tasks like 3D object recognition [6], point cloud registration
[7], and 3D object categorization [8]. The first stage in local feature-based
matching is detecting the keypoints on the surface and generating descriptors
for every keypoints. Generating the vector feature descriptors surrounding
each keypoint is computationally demanding, and thus, detecting keypoints is
required to reduce the amount of computation required. In the second stage,
crude initial matches are generated for marking the similarities between two
3D shapes. However, the initial matches have a high number of false positives
due to reasons such as the residual errors from the proceeding steps like
errors in keypoint localization, noise, point density variations, occlusion,
overlaps, etc. To ensure accurate transformation estimation and hypothesis
formulation, outliers must be filtered out from the initial matching, making
correspondence grouping (CG) important [9].

CG can play an important role in 3D model based plant phenotyping and
automation of the farming process. Certain applications that require
correspondence matching include creating 3D models, plant species recognition,
recognition of plant parts, unsupervised template learning for fine-grained
plant structures , etc. In autonomous farming also, techniques based on
correspondence matching can be used in many ways. For instance, a general
model of a harvesting machine can be produced in large numbers and then any
farmer having no robotics skills can teach it on-site to identify and pick the
required kind of fruit. In the same way, an autonomous weeding machine can be
mass produced, then taught to distinguish weeds from the crops at different
kinds of farms of widely, improving farming efficiency.

Although, there are studies in the existing literature comparing the
performance of 3D CG algorithms [10], to the best of our knowledge, there is
no work in the existing literature that deals specifically with the
performance evaluation of the CG algorithms as applied to 3D point cloud data
for plant health classification and plant growth analysis. Working with 3D
point cloud data in the case of plants offers special challenges:

  * •

Plants are made up of very fine structures that makes it very difficult to
make a perfect scan. Thus the keypoint detectors local feature descriptors
have to be robust to noise and holes present in the point cloud data.

  * •

The lighting conditions in farming keep changing, offering different
illumination conditions. This also changes the color and the texture of the
plants, thus, making it very difficult to use the color information for
decision making.

In summary, in this paper, we consider four most common 3D CG algorithms, viz,
Nearest Neighbor Similarity Ratio (NNSR) [11], [12], Random Sample Consensus
(RANSAC) [13], [9], Geometric Consistency (GC) [14], [15] and 3D Hough Voting
(3DHV) [16], [17]. We have also implemented 3D MLEASAC, which combines Maximum
Likelihood Estimation (MLE) and 3D RANSAC and is motivated by the approach
taken in [18], [19] for 2D images. Fig.1 is represented the pipeline of our
work.

The remainder of the paper is structured as follows. Section II briefly
overviews the related work. We give a description of the 3DCG algorithms that
we have considered for evaluation in Section III. performance evaluation is
provided in Section IV. Finally, in Section V, we provide some perspectives
and conclude the paper.

##  2 Related work

3DCG is still an active area in computer vision. 3DCG works based on local
shape feature matching and is commonly used in other applications such as
object recognition, pose estimation, point cloud registration, etc. Image
based correspondences features have been explored in depth in the existing
literature. Torresani [20] presented a method for finding correspondences
between parse image features by graph matching. They solve the matching
problem by minimizing a image matching objective function which depend on the
appearance, geometric and spatial features. This method for real-world image
matching is able to achieve global optimal by an unknown non-rigid mapping.
Lowe [21] is described how to use SIFT features for object recognition in
image. In this paper, they investigated Nearest Neighborhood algorithm for
matching individual feature to a database of feature base on Euclidean
distance of their feature vector. They found good matches in a large number of
features with using distinctive descriptors. This method is extended to 3D in
[9]. In [22] is modified RANSAC method for 3D object recognition and pose
estimation. They obtained optimum number of iteration by using the combination
algorithms of descriptor, hash table and RANSAC. In [14], Chen et al. used
geometric constraints for 3DCG algorithm in 3D object recognition. They
improved the version of Johnson et al. [23] algorithm for clustering the
inliers into disjointed clusters. Tombari et al. [16] introduced 3D Hough
Voting for 3D object recognition by using Hough Transform [24]. In the
proposed approach each correspondence can cast a vote for presence of the
unique reference point of object in the 3D Hough space simultaneously. Anders
et al. [9] introduced two stage local and global method for casting a vote for
finding 3DCG. On the local stage, they caste a vote for finding
correspondences by geometric consistency and in the global stage, they cast a
vote by covariance constraint. At the final stage, they compared their
algorithm with existing ones.

##  3 3D correspondence grouping algorithms

In this section, we provide a brief overview the 3DCG algorithms. Given a
model $\mathcal{S}$ and a scene $\mathcal{S^{\prime}}$, where an raw
correspondence set $\mathcal{C}$ is obtained after comparing the feature sets
$\mathcal{F}$ and $\mathcal{F^{\prime}}$ belongng to $\mathcal{S}$ and
$\mathcal{S^{\prime}}$, respectively, the goal is to obtain a subset
$\mathcal{C}_{inlier}$ of $\mathcal{C}$ that identifies the correct
correspondences between $\mathcal{S}$ and $\mathcal{S^{\prime}}$, i.e., the
inliers. An element in $\mathcal{C}$ is given as:
$c=\\{p,p^{\prime},s_{\mathcal{F}}\left(f,f^{\prime}\right)\\}$, where $p$ and
$p^{\prime}$ are keypoints belonging to the model (source) and the scene
(target), respectively, $f$ and $f^{\prime}$ are the descriptors for these
keypoints, and $s_{\mathcal{F}}\left(f,f^{\prime}\right)$ is the similarity
score assigned to $c$. We now briefly explain the 3DCG algorithms considered
in this paper.

###  3.1 Nearest neighbor similarity ratio

NNSR is a basic method for determining the initial correspondences [11]. This
method is based on Lowe’s ratio [21]. It judges match features base on the
ratio of two closest neighbors and obtained by Euclidean ($L2$) distance
between an invariant feature vector $f\in\mathcal{F}$ with the nearest
$f^{\prime}_{1}\in\mathcal{F^{\prime}}$ and second nearest matching feature
$f^{\prime}_{2}\in\mathcal{F^{\prime}}$. A correspondence will be as an inlier
if the Euclidean distance between $f$, $f^{\prime}_{1}$ and $f$,
$f^{\prime}_{2}$ satisfies:

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}\frac{\left\|{f - f'_{1}}
\right\|_{L2}}{\left\|{f - f'_{2}} \right\|_{L2}} \leq{t_{nnsr}},
{}$\hfil&{\vrule width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(1)}\cr}$

where, $t_{nnsr}$ is the threshold. The NNSR will be used as a baseline for
all comparison purposes.

Figure 1: Performance comparison of 3D correspondence grouping algorithms

###  3.2 RANSAC

RANSAC is an iterative method for judging the correctness of a correspondence
set by calculating the number of inliers returned in each iteration [13]. The
RANSAC method uses random samples from the initial correspondences set
$\mathcal{C}$ in each iteration and finds the transformation $T$ corresponding
to those samples. The steps are repeated $N_{ransac}$ iteration and the
transformation is compared with a threshold $t_{ransac}$ . The transformation
returning the maximum number of inliers is judged the optimal $T^{*}$. The
correspondences in $\mathcal{C}$ that do not agree with $T^{*}$ are classified
as outliers.

###  3.3 Geometric consistency

CG method is based on geometric constrains for finding inliers and is not
dependent on the local feature space [14]. The algorithm looks for consistency
and compatibility in terms of consistent cluster. After finding initial
matching set $\mathcal{C}$, we find the distance between keypoints from the
seed in the model and scene and then we find $L2$-distances between the
pairwise correspondences to compare with the geometric constraint threshold
$t_{gc}$. Finally, we choose the largest cluster as optimum one and classify
the corresponding correspondences as inliers. The compatibility score for the
given correspondences pair $c_{1}$ and $c_{2}$ is given as:

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}d\left( {c_{1},c_{2}} \right)
= \left| {d\left( {p_{1},p_{2}} \right) - d\left( {p_{1}',p_{2}'} \right)}
\right| \leq{t_{gc}}, {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(2)}\cr}$

where ${d\left(p_{1},p_{2}\right)=\left\|{(p_{1}-p_{2}}|\right|_{L2}}$.

###  3.4 3D Hough voting

3DHV uses features by casting a vote in 3D Hough space for finding inliers
[16]. In this method, after finding initial matching set $\mathcal{C}$, for
each keypoint $p_{i}$ in the source shape (model) that is in $\mathcal{C}$ ,
the vector between source shape centeroid $C_{S}$ and $p_{i}$ is computed in
the global reference frame (GRF):

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}{\bf{V}}_{i,G}^{S} = {C_{S}}
-p_{i}, {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(3)}\cr}$

After that ${\bf{V}}_{i,G}^{S}$ is transformed to local reference frame (LRF),
as we want our method to be translation and rotation invariant.

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}{\bf{V}}_{i,L}^{S} =
{\bf{R}}_{i}^{S} \cdot{\bf{V}}_{i,G}^{S}, {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(4)}\cr}$

Where ${\bf{R}}_{i}^{S}$ is rotation matrix, and each line in
${\bf{R}}_{i}^{S}$ being a unit vector of the LRF of $p_{i}$. Similarly, we
can obtain ${\bf{V}}_{i,G}^{S^{\prime}}$ corresponding to $p_{i}^{\prime}$.
Then, ${{\bf{V}}_{i,L}^{S}=\bf{V}}_{i,L}^{S^{\prime}}$, because of the
invariance in LRF. The vector ${\bf{V}}_{i,G}^{S^{\prime}}$ is computed by
transforming ${\bf{V}}_{i,L}^{S^{\prime}}$ into GRF.

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}{\bf{V}}_{i,G}^{S'} =
{\bf{R}}_{i}^{S'}\cdot{\bf{V}}_{i,L}^{S'} + {p_{i}}', {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(5)}\cr}$

Eventually the feature ${f_{i}}^{\prime}$ could cast a vote in 3D Hough space
through the vector ${\bf{V}}_{i,G}^{S^{\prime}}$. The inliers correspond to
the peak of Hough space.

###  3.5 3D MLESAC

MLEASAC combines 3D RANSAC with MLE for finding inliers. The novelty of our
algorithm is generalization of RANSAC in 3D. This algorithm modifies RANSAC in
a way such that we look for the transformation that maximizes the likelihood
function rather than maximizing the number of inliers. Our implementation of
3D MLESAC consists of the following basic steps.

Similar to 3D RANSAC, given $N_{mlesac}$ iterations, in each iteration, five
random components are sampled from the initial matching $\mathcal{C}$. We use
these samples to calculate the affine transformation $T$ and we transform all
model keypoint that are in $\mathcal{C}$ according to $\mathcal{C}$. We then
define the likelihood function for each $i=1,...,N$, with $N$ being the number
of column entries in $\mathcal{C}$ as:

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}\Pr\left( d_{i} \right) =
\gamma\frac{1}{{\sqrt{2\pi{\sigma^{2}}} }}\exp\left( { -
\frac{{{d_{i}^{2}}}}{{2{\sigma^{2}}}}} \right) + \left( {1 - \gamma}
\right)\frac{1}{\nu}, {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(6)}\cr}$

where, $d_{i}$ is the Euclidean distance between transformed model keypoints
$p_{i}$ with matched points $p_{i}^{\prime}$ in scene, $\nu$ is a constant,
$\sigma$ is the standard deviation, and $\gamma$ is the mixing parameter,
which has been calculated according to the method described in [18].

Since, distances $d_{i}$ between $p_{i}^{S}$ and $p_{i^{\prime}}^{S^{\prime}}$
are independent of each other, the joint likelihood function for all $d_{i}s$
given as:

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}\Pr\left( d \right) =
\prod\limits_{i = 1}^{N} {\left( {\gamma\frac{1}{{\sqrt{2\pi{\sigma^{2}}}
}}\exp\left( { - \frac{{{d_{i}}^{2}}}{{2{\sigma^{2}}}}} \right) + \left( {1 -
\gamma} \right)\frac{1}{\nu}} \right)} .\hphantom{(7)} {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(7)}\cr}$

To simplify the calculations, we take logarithm on both sides of the above,
giving the negative log-likelihood function as;

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREc\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTc\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle- L = - \sum\limits_{i = 1}^{N}
{\log\left( {\gamma\frac{1}{{\sqrt{2\pi{\sigma^{2}}} }}\exp\left( { -
\frac{{{d_{i}}^{2}}}{{2{\sigma^{2}}}}} \right) + \left( {1 - \gamma}
\right)\frac{1}{\nu}} \right)} .\hphantom{(8)} $\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(8)}\cr}$

Finally, the transformation that results in the minimum of negative log-
likelihood function is compared with $t_{mlesac}$ and it is considered the
optimum transformation $T^{*}$. The correspondences in $\mathcal{C}$ that do
not agree with $T^{*}$ are classified as outliers.

Note that the first term on the RHS of (3.5), refers to the correspondences
that are inliers while the second term corresponds to the outliers, and is
assumed to have a uniform distribution.

(a) Plant

(b) B3R

(c) U3M

Figure 2: Visualization (snapshot) of sample model (right) and scene (left)
from each dataset

##  4 Experimental evaluation and discussion

Figure 3: Performance of the evaluated CG algorithms in terms of precision and
recall for the plant dataset

Figure 4: Performance of the evaluated CG algorithms in terms of precision and
recall for the B3R dataset

In this section, we evaluate and provide the implementation details of all the
algorithms discussed in the previous section on two chosen benchmarks and a
plant dataset. Fig.2 shows the samples of models and scenes in different
datasets. These datasets have scenes with different noise levels, point
densities, and partial overlaps. We determine the correctness of a
correspondence using the ground truth affine transformation $\textbf{T}_{GT}$
between the model and the scene. A given correspondence $c$ between points
$p$, $p^{\prime}$ is declared correct if:

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}\left| {{{\bf{T}}_{GT}}\left(
p \right) - p'} \right| \leq t, {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(9)}\cr}$

Where $t$ is the threshold. The calculated inliers by all algorithms are
measured in terms of the precision and recall criteria which are defined in
terms of the inlier set ${{C_{in}}}$, the correctly judged set ${{C_{cor}}}$,
and the ground truth set ${{C_{GT}}}$ as

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}{\rm{Precision}} =
\frac{{\left| {{C_{cor}}} \right|}}{{\left| {{C_{in}}} \right|}},
{}$\hfil&{\vrule width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(10)}\cr}$

and

$\halign to=0.0pt{\relax\global\advance\@IEEEeqncolcnt by
1\relax\begingroup\csname
@IEEEeqnarraycolPREC\endcsname#\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\csname
@IEEEeqnarraycolPOSTC\endcsname\relax\relax\relax\relax\relax\relax\relax\relax\relax\relax\endgroup&\global\advance\@IEEEeqncolcnt
by 1\relax\bgroup#\egroup&\global\advance\@IEEEeqncolcnt by 1\relax\hbox
to\z@\bgroup\hss#\egroup\cr\hfil$\displaystyle{}{\rm{Recall}} = \frac{{\left|
{{C_{cor}}} \right|}}{{\left| {{C_{GT}}} \right|}}. {}$\hfil&{\vrule
width=0.0pt,height=0.0pt,depth=0.0pt}&0.0pt{\hss(11)}\cr}$

###  4.1 Dataset and experimental setup

####  4.1.1 Plant dataset

The models for this dataset have been obtained from the dataset containing 3D
scans of plant shoot architectures provided by Salk Institute of Biological
Science [25]. This dataset contains $559$ 3D plant shoot architectures from 4
species (Arabidopsis, tomato, tobacco, and sorghum) scanned under multiple
conditions (ambient light (control), heat, high-light, shade, drought). Each
plant was scanned every 1 or 2 days through roughly 20 days of development.
In our work we have used 26 model of sorghum, tobacco and tomato plants in
three different conditions of control, shade and heat in different growth
stage and different variety because this data has the same categories and the
same number of replicates which makes it easier to compare the results. The
scenes are created by rotating the models and by using different levels of
noise and down-sampling on these models, giving a total of $80$ scenes for
each model under consideration.

####  4.1.2 B3R dataset

The Bologna 3D retrieval (B3R) dataset [26], with 6 models, is used to
evaluation the algorithms robustness with respect to noise and point
densities. The models are taken from the Stanford Repository, however the
scenes here are created by rotating the models, and then adding $10$ levels of
noise and $8$ levels of downsamplig, giving a total of $80$ scenes for each
model.

Table 1:  Implementation platforms of the used techniques Technique |  Platform
---|---
ISS Keypoints Detector |  PCL
3D SIFT Descriptor |  MATLAB
NNSR |  MATLAB
RANSAC |  MATLAB
GC |  MATLAB
3D HV |  MATLAB
3D MLESAC |  MATLAB
Table 2: Parameters of algorithms NNSR |  $t_{nnsr}$ $20$
---|---
RANSAC |  $t_{ransac}$ 0.01 $N_{ran}$ 1000
GC |  $t_{GC}$ 0.01
3D MLESAC |  $t_{mlesac}$ 0.01 $N_{mlesac}$ 1000

Figure 5: Performance of the evaluated CG algorithms in terms of precision and
recall for the U3M dataset

(a) NNSR- B3R

(b) NNSR- Plant

(c) NNSR- U3M

(d) RANSAC- B3R

(e) RANSAC- Plant

(f) RANSAC- U3M

(g) MLESAC- B3R

(h) MLESAC- Plant

(i) MLESAC- U3M

(j) 3DHV- B3R

(k) 3DHV- Plant

(l) 3DHV- U3M

(m) GC- B3R

(n) GC- Plant

(o) GC- U3M

Figure 6: Visualization of the evaluated CG algorithms for the three datasets

####  4.1.3 U3M dataset

The UWA 3D modeling (U3M) [27] dataset belongs to the point cloud (2.5D view)
registration scenario. There are $22$, $16$, $16$ and $21$ $2.5D$ views
respectively captured from the Chef, Chicken, T-rex, and Parasaurolophus
models. There scenes represent different levels of overlaps. In our
experiments, we have taken first $16$ views of each model.

####  4.1.4 Implementations details

In this paper, we have used 3D Intrinsic shape signature (ISS) for detecting
keypoints and 3DSIFT for describing the local features. This choice is
motivated by the results obtained in our previous work, where this detector-
descriptor pair outperformed other combinations of popular detectors and
descriptors in case of plant 3D point clouds. In the keypoint detection stage,
around $3\%$ of points in the point cloud are chosen as interest points, with
the non maximum radius set to round $4$ point cloud resolution ($pr$). For the
3DSIFT descriptor, the search radius is set to $8pr$ and the size of feature
descriptor is 128. NNSR is used for finding initial matching. RANSAC, GC,
3DHV, and 3D MLESAC are used as 3DCG algorithms. Table.1 gives the
implementation platforms of the used algorithms and the parameters used in all
algorithms are listed in Table.2.

###  4.2 Results and discussions

####  4.2.1 Performance on the plant dataset

The performance of the 3DCG algorithms for the plant dataset are shown in
Fig.3. Noise is expected to degrade on the discriminative capabilities of
feature descriptor, thus creating a certain amount of false matches. This is
obvious in the plots obtained for precision and recall for all the 3DCG
algorithms. The precision decreases with increasing noise variance. Similar
trend is observed for the recall values, however, the deterioration is
smaller, as compared to the precision. Down sampling the scenes also has
negative impact on the performance, as seen in the precision and recall plots.
This can be attributed to the false matches that become more common as the
number of points become less. As indicated by the plots, RANSAC and MLESAC
appear to be the best two ones among all evaluated algorithms, considering
their overall precision and recall performance, however, MLESAC generally
outperforms RANSAC, while taking less time. This indicates that for larger
point clouds of objects having highly complex shapes, MLESAC is a better
choice as a 3DCG algorithm.

####  4.2.2 Performance on the B3R dataset

The performance of the 3DCG algorithms for the B3R dataset are shown in Fig.4.
It follows similar trends, as in the case of the plant as seen in Fig.3. As
indicated by the plots here also, RANSAC and MLESAC appear to be the best two
ones among all evaluated algorithms, considering their overall precision and
recall performance, with RANSAC performing slightly better. However, in our
experiments, MLESAC is quite faster as compared to RANSAC.

####  4.2.3 Performance on the U3M dataset

A common trend to all algorithms is that their performance generally degrades
as the degree of overlap drops as depicted in Fig.5. This is owing to the fact
that the ratio of outliers in the initial correspondence set is closely
correlated to the ratio of overlapping regions. MLESAC outperform in all
algorithm in term of precious and recall.

####  4.2.4 Visualization

Fig.6 shows some visual results of 3DCG algorithms. These figures display the
strength of each method. For example in NNSR method, The number of outliers
are more than other methods, because of occlusion and overlap in the feature
matching. We can also see the location and number of correspondences are also
different.

##  5 Conclusions

In this paper, we compared the performance of different 3D correspondence
grouping algorithms in context of 3D plant point clouds and also provided the
corresponding performance evolution of some standard publically available
datasets for comparison. The performance was evaluated in terms of precision
and recall, the results were presented in terms of the corresponding plots.
Experimental results show that the RANSAC and MLESAC perform quite closely and
are better than 3D Hough Voting and Geometric Consistency. Also, in general,
MLESAC is much faster compared to RANSAC, making it a better algorithm for
applications involving 3D point clouds of plants as they are usually very
large in size. In future works, we intend to prepare our own dataset of 3D
point clouds of crops such as rice and wheat and will try to extend the
analysis further to perform plant identification for plant phenomics
applications.

## References

  * [1] Andrew J Challinor, J Watson, David B Lobell, SM Howden, DR Smith, and Netra Chhetri.  A meta-analysis of crop yield under climate change and adaptation.  Nature Climate Change, 4(4):287, 2014.
  * [2] Noah Fahlgren, Malia A Gehan, and Ivan Baxter.  Lights, camera, action: high-throughput plant phenotyping is ready for a close-up.  Current opinion in plant biology, 24:93–99, 2015.
  * [3] Massimo Minervini, Hanno Scharr, and Sotirios A Tsaftaris.  Image analysis: the new bottleneck in plant phenotyping [applications corner].  IEEE signal processing magazine, 32(4):126–131, 2015.
  * [4] Rachana Panwar, Kusha Goyal, Nilay Pandey, and Nitin Khanna.  Imaging system for classification of local flora of uttarakhand region.  In 2014 International Conference on Power, Control and Embedded Systems (ICPCES), pages 1–6. IEEE, 2014.
  * [5] Suxing Liu, Lucia Acosta-Gamboa, Xiuzhen Huang, and Argelia Lorence.  Novel low cost 3d surface model reconstruction system for plant phenotyping.  Journal of Imaging, 3(3):39, 2017.
  * [6] Anders G Buch, Henrik G Petersen, and Norbert Krüger.  Local shape feature fusion for improved matching, pose estimation and 3d object recognition.  SpringerPlus, 5(1):297, 2016.
  * [7] Huan Lei, Guang Jiang, and Long Quan.  Fast descriptors and correspondence propagation for robust global point cloud registration.  IEEE Transactions on Image Processing, 26(8):3614–3623, 2017.
  * [8] Yu Xiang, Wongun Choi, Yuanqing Lin, and Silvio Savarese.  Data-driven 3d voxel patterns for object category recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1903–1911, 2015.
  * [9] Anders Glent Buch, Yang Yang, Norbert Kruger, and Henrik Gordon Petersen.  In search of inliers: 3d correspondence by local and global voting.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2067–2074, 2014.
  * [10] Jiaqi Yang, Ke Xian, Yang Xiao, and Zhiguo Cao.  Performance evaluation of 3d correspondence grouping algorithms.  In 3D Vision (3DV), 2017 International Conference on, pages 467–476. IEEE, 2017.
  * [11] Marius Muja and David G Lowe.  Fast matching of binary features.  In 2012 Ninth conference on computer and robot vision, pages 404–410. IEEE, 2012.
  * [12] Antti Hietanen, Jussi Halme, Anders Glent Buch, Jyrki Latokartano, and J-K Kämäräinen.  Robustifying correspondence based 6d object pose estimation.  In 2017 IEEE International Conference on Robotics and Automation (ICRA), pages 739–745. IEEE, 2017.
  * [13] Martin A Fischler and Robert C Bolles.  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM, 24(6):381–395, 1981.
  * [14] Hui Chen and Bir Bhanu.  3d free-form object recognition in range images using local surface patches.  Pattern Recognition Letters, 28(10):1252–1262, 2007.
  * [15] Jiaqi Yang, Yang Xiao, Zhiguo Cao, and Weidong Yang.  Ranking 3d feature correspondences via consistency voting.  Pattern Recognition Letters, 117:1–8, 2019.
  * [16] Federico Tombari and Luigi Di Stefano.  Object recognition in 3d scenes with occlusions and clutter by hough voting.  In 2010 Fourth Pacific-Rim Symposium on Image and Video Technology, pages 349–355. IEEE, 2010.
  * [17] Anders Glent Buch, Lilita Kiforenko, and Dirk Kraft.  Rotational subgroup voting and pose clustering for robust 3d object recognition.  In 2017 IEEE International Conference on Computer Vision (ICCV), pages 4137–4145. IEEE, 2017.
  * [18] Philip HS Torr and Andrew Zisserman.  Mlesac: A new robust estimator with application to estimating image geometry.  Computer vision and image understanding, 78(1):138–156, 2000.
  * [19] Lijin Deng, Yan Piao, and Shuo Liu.  Research on sift image matching based on mlesac algorithm.  In Proceedings of the 2nd International Conference on Digital Signal Processing, pages 17–21. ACM, 2018.
  * [20] Lorenzo Torresani, Vladimir Kolmogorov, and Carsten Rother.  Feature correspondence via graph matching: Models and global optimization.  In European Conference on Computer Vision, pages 596–609. Springer, 2008.
  * [21] David G Lowe.  Distinctive image features from scale-invariant keypoints.  International journal of computer vision, 60(2):91–110, 2004.
  * [22] Chavdar Papazov and Darius Burschka.  An efficient ransac for 3d object recognition in noisy and occluded scenes.  In Asian Conference on Computer Vision, pages 135–148. Springer, 2010.
  * [23] A Johnson and M Hebert.  Surface matching for object recognition in complex 3-d scenes. to appear in.  Image and Vision Computing, 1998.
  * [24] Paul VC Hough.  Method and means for recognizing complex patterns, December 18 1962.  US Patent 3,069,654.
  * [25] Saket Navlakha.  3d scans of plant shoot architectures, Jul 2017.
  * [26] Federico Tombari, Samuele Salti, and Luigi Di Stefano.  Performance evaluation of 3d keypoint detectors.  International Journal of Computer Vision, 102(1-3):198–220, 2013\.
  * [27] Ajmal S Mian, Mohammed Bennamoun, and Robyn A Owens.  A novel representation and feature matching algorithm for automatic pairwise registration of range images.  International Journal of Computer Vision, 66(1):19–40, 2006.


# Perspectives on Reproducibility and Sustainability of Open-Source Scientific
Software from Seven Years of the Dedalus Project
A white paper submitted to the National Academies of Sciences, Engineering,
and Medicine’s Best Practices for a Future Open Code Policy for NASA Space
Science

missing  Bates College Jeffrey S. Oishi<EMAIL_ADDRESS>Benjamin P. Brown
University of Colorado, Boulder  Keaton J. Burns  MIT  Daniel Lecoanet
Princeton  Geoffrey M. Vasil  University of Sydney

Jeffrey S. Oishi

##  1 Introduction

As the Science Mission Directorate contemplates establishing an open code
policy, we consider it timely to share our experiences as the developers of
the open-source partial differential equation solver Dedalus. Dedalus is a
flexible framework for solving partial differential equations. Its development
team primarily uses it for studying stellar and planetary astrophysics.
Dedalus was developed originally for astrophysical fluid dynamics (AFD),
though it has found a much broader user base, including applied
mathematicians, plasma physicists, and oceanographers. Here, we will focus on
issues related to open-source software from the perspective of AFD. We use the
term AFD with the understanding that astrophysics simulations are inherently
multi-physics: fluid dynamics coupled with some combination of gravitational
dynamics, radiation transfer, relativity, and magnetic fields. In practice, a
few well-known open-source simulation packages represent a large fraction of
published work in the field. However, we will argue that an open-code policy
should encompass not just these large simulation codes, but also the input
files and analysis scripts. It is our interest that NASA adopt an open-code
policy because without it, reproducibility in computational science is
needlessly hampered.

With the caveat that we have not performed bibliometric analyses, it seems
reasonable to assert that the majority of papers published in AFD are done
using open-source codes, especially in star formation, accretion disk physics,
and large-scale structure formation. _Athena_ , _Enzo_ , _Pencil Code_ ,
_Pluto_ , and _Gadget 2_ are all widely used open-source tools. In stellar
AFD, there are fewer examples of open-source tools; MagIC and MAESTRO are
important exceptions. For many of these simulation codes, the preferred method
of in-depth data analysis is via yt, a pioneering open-source project that
introduced many of the important development and community-building concepts
we use in Dedalus. These codes tend to be very large: Athena has $\simeq
80,000$ SLOC; Enzo is over $120,000$ lines, and MAESTRO is over $161,000$
lines111All source lines of code reported here are measured using the the sloc
tool on publicly available repositories of each code as of 3 January 2018!
Such complexity is not unique to open-source tools; two of us have significant
experience with closed-source astrophysics codes (_ORION_ and _ASH_) of
comparable size. The sheer size of these codes presents challenges to
reproducibility: it takes a significant amount of effort to even be able to
build and run them.

The sub-community of astrophysical fluid dynamics modeling thus tends to
embrace open-source even without an official policy. However, this is, in our
view, not sufficient, especially as it can vary from field to field: many of
the most important spherical dynamo codes for studying stellar and planetary
magnetic fields remain closed-source (e.g. _Rayleigh_ 222At the time this
paper was submitted on 12 January 2018, _Rayleigh_ was closed source. Since
that time, it has been made available under the GPL at
https://github.com/geodynamics/Rayleigh, _ASH_ , various derivatives of the
Glatzmaier/Gilman code). The simulation code is only one part of a long
pipeline that starts with an idea and ends with a published result. An open-
code mandate should eliminate any closed-source bottlenecks that could hamper
what is, in our view, the cornerstone of the scientific computing enterprise:
repeatability. That is, simply because the largest piece of software is
available, if the input files and analysis scripts are not also public,
reproducibility becomes impossible.

Dedalus is written primarily in python; it has a very small code base of
approximately 7000 source lines of code (SLOC). It originated in an effort by
one of us to develop a easily modifiable, modern fluid dynamics solver for
studying magnetohydrodynamic turbulence problems in AFD. The first commit to
the project was in 2010. In 2012, the project expanded from two people to
five. Since then, nine people have made commits to the code base, including
four not affiliated with the core development team. We chose the term “core
team” to represent not a privileged set of developers but a reflection of the
current number of highly active participants in the development and future of
Dedalus. We have logged 46 pull requests, representing peer-reviewed code
contributions. The code is licensed under the GNU General Public License
Version 3 (GPL). We chose this license to ensure that the source code always
be available for inspection and auditing by any interested parties.

We do not ask users to register before receiving the code. Doing so would, in
our view, constitute a barrier to access without providing any substantive
estimate on the number of users. Instead, we can estimate the size of our
community by a number of metrics. We have 177 unique users participating in
the dedalus-users mailing list for discussions of using the code, and 13
unique users participating in the dedalus-dev for discussions of development
of new features. Dedalus has been used in at least thirteen peer-reviewed
publications (Bordwell et al., 2018; Anders & Brown, 2017; Couston et al.,
2017; Clark & Oishi, 2017a, b; Lecoanet et al., 2017, 2016a; Vasil et al.,
2016; Davidovits & Fisch, 2016a, b; Lecoanet et al., 2016b, 2015, 2014).

Here, we argue that the benefits of open-source projects far exceeds their
costs; we outline our view on exactly what those benefits and costs are, and
we suggest that in order to maximize scientific returns on software
development, _open-source_ is not enough. The best use comes from _community-
driven_ software tools that blur lines between users and developers. We
strongly recommend that Science Mission Directorate adapt an open code policy
that embraces _explicitly licensed_ , discoverable software.

##  2 Reproducibility

Reproducibility is the cornerstone of scientific research, but there are
serious concerns regarding the ability to reproduce computational results (see
Stodden et al., 2012, and references therein). There are many definitions of
reproducibility; here we focus on the ability of outside researchers to
replicate published results. Open-source codes support this reproducibility by
allowing any researcher the ability to, in theory, use the same tool as the
original study. However, as mentioned in the introduction, modern AFD codes
are extremely large; though Dedalus is an order of magnitude smaller than
others, it gives users much more freedom in specifying problems and thus
remains similarly complex to use. This complexity makes the cultivation of a
_community_ of contributors important in ensuring the ability of outside
researchers to reproduce results.

We actively work to build a community of practice influenced by the _Enzo_ and
_yt_ communities, as outlined in Turk (2013). This gives researchers a way not
only to access the source code but also ask questions and contribute to the
code. However, community of this sort is still not enough to ensure
reproducibility. In addition, several members of our community have made
publicly accessible the version control repositories containing _all_ material
required to reproduce the results published in Clark & Oishi (2017a,
b)333https://bitbucket.org/jsoishi/weakly_nonlinear_mri, Anders & Brown
(2017)444https://bitbucket.org/exoweather/polytrope, and Bordwell et al.
(2018)555https://bitbucket.org/exoweather/polytrope. By publishing analysis
scripts and input files, these repositories allow outside researchers to
download Dedalus, run it using the same inputs, and then analyze the output to
produce the plots found in those papers. We believe that these are the first
steps in a process that can improve the quality and reproducibility of the
science results while offering an opportunity for outside researchers to build
on those results. Researchers can test the claims made in each paper and then
extend those claims to answer new questions that the authors may not have
envisioned. For large simulations, the amount of computational resources
required to _directly_ re-run all the simulations may be prohibitive. In this
case, mechanisms to allow access to the output data from the simulations
should be made available. We discuss infrastructure to do this in section 5.

In addition to reproducing published results–an action we suggest is far too
rarely performed–open-source policies help promote verification and validation
testing. In AFD, verification tests are plentiful; there is a strong tradition
of benchmark problems and inter-code comparisons in AFD (e.g. Birn et al.,
2001; Marti et al., 2014; Kim et al., 2014; Lecoanet et al., 2016b, among many
others). However, direct validation is nearly impossible for AFD: there are
precious few astrophysical experiments that can be done repeatably in the lab.
In lieu of this, AFD codes can perform validation tests against standard
terrestrial fluid dynamics problems. Here again, Dedalus has been
significantly aided by its community: by eliminating boundaries between users
and developers, we actively encourage the development of validation problems
which are then included in the code base (e.g. Rayleigh-Benard Convection). We
have found that these validation problems are often the starting point for new
users who use them as examples and extend them to their problem. Thus,
validation problems also communicate best practices for using the tool in a
way that has been very beneficial for the entire community.

We would suggest that any open code policy adopted by the SMD explicitly
consider the importance of releasing input files, analysis scripts, and build-
time parameters. We elaborate on the need to consider the software stack upon
which the simulation and analysis tools rely in section 4.

##  3 Maintenance and Support

Providing open-source code is not without cost. Direct costs are quite
minimal; they are limited to web hosting. The source for Dedalus is hosted at
Bitbucket, a commercial provider of distributed version control services that
allows open-source projects to be hosted for free. An often-discussed cost of
open-source research products is the “need” to provide user support for them.
However, we have found that, despite a community that spans many scientific
disciplines, the amount of direct support we do is small and it occasionally
leads to significant improvements to the code base itself. Indirect support,
in the form of responding to bug reports and implementing enhancements
suggested by community members is also small, amounting to no more than a few
person hours per month. Moreover, by encouraging a community of practice in
which there is not a clear distinction between “users” and “developers”, we
have mechanisms for anyone to provide fixes and receive credit in the form of
having their identity embedded along with the changes in the version control
log.

Figure 1 shows the number of messages per month sent to the the two mailing
lists we maintain, dedalus-users and dedalus-dev. dedalus-users shows a steady
(roughly linear) average growth over the past two and a half years, with a
peak of just over 80 messages in one month. Making a rough assessment of 15
minutes per message, this represents approximately 10 person-hours per month
on the part of the core team, if we assume that half of the messages are
questions and the other half responses. Given a core team of 5 researchers,
this represents an upper limit of 2 hours per month per developer. However,
non-core members of our community often _answers_ questions on this list in
addition to asking them (and core team members also _ask_ questions too!). The
groups themselves are archived and searchable, making them a knowledge base
for users. Finally, we note that often questions on the user mailing list lead
to the development of new example problems or features: recently, a question
on dedalus-users about how to handle the arguments to a certain kind of
function
object666https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!msg/dedalus-
users/Q63KhHSlfrI/fcFULv12AAAJ led to the development of a new tutorial
problem777https://bitbucket.org/dedalus-project/dedalus/pull-requests/44. This
means that at least some of the time spent participating in answering
questions on the mailing list is not a pure cost in terms of an additional
burden on the core team–it can yield significant returns for the entire
community.

Figure 1: Message traffic on the dedalus-users and dedalus-dev google groups.

##  4 Software stack

Complex AFD software relies on a non-trivial software stack to run. The stack
can include compilers (e.g. C, C++, FORTRAN), interpreters (python), and
libraries (FFTW, HDF5, BLAS, LAPACK). Besides providing input parameter files,
analysis scripts, and the simulation code itself, one must also provide the
entire software stack to outside researchers. In the case of Dedalus, the
software stack is entirely open-source. This has both positive and negative
effects on usability and workflow. It ensures there are no closed source or
expensive, proprietary bottlenecks, and it ties to a rich ecosystem of tools,
including Jupyter notebooks, docker containers, and integration with cloud
computing services (e.g. CoCalc and Amazon Web Services (AWS)).

However, a significant drawback is the complexity of _installing_ this stack.
One of the greatest challenges we face is the distribution of software on
heterogeneous computing environments, including laptops running macOS,
university clusters from hundreds to thousands of cores, and national class
supercomputing facilities in the United States (NASA, NSF, DoE) and Europe
(PRACE). In addition to merely installing Dedalus on these systems, an
additional challenge lies in optimizing for performance and verification
testing to ensure accuracy.

To support this effort, we have developed simple scripts to install both
Dedalus and its stack. This approach is cumbersome to maintain and inefficient
in both install time and disk space, and so we are pursuing several avenues to
replace the scripts. A non-core team member has contributed a series of tools
to allow Dedalus to integrate with the popular open-source miniconda package
management system888https://github.com/pkgw/dedalus-builder.

In order to allow users to assess optimization and performance, we have
published test problems with scaling results and expected overall performance
on key machines including the NASA Pleiades
system999https://bitbucket.org/exoweather/incompressible_ns_tg.

##  5 Looking to the Future

As mentioned in section 2, there are often published results relying on
calculations at a scale far too large to allow outside researchers to
reproduce directly. In this case, we believe that the solution is to make
datasets available along with analysis scripts and parameter files. The latter
are particularly important should Moore’s law continue to provide exponential
growth in computational power; what is infeasible to reproduce today may not
be in a few years’ time. An open-code policy adapted by NASA should consider
not only the requirement to make this available but also provide guidance on
the use of _middleware_ designed to ease the burden of doing so. One
particularly promising effort is the Whole Tale Project, which promises to
link data, analysis scripts, and computing frontends (e.g. Jupyter notebook)
to allow “living publications” in which any researcher can re-run analysis on
datasets to reproduce results.

Second, an open-code policy should also be focused on a future where computing
hardware is increasingly fungible. Here, we imagine an emphasis on
containerization solutions for software stacks including the open-source AFD
simulation tools themselves. The use Docker containers for this purpose in
industry is widespread, and astrophysics should follow suit. While some
challenges remain for high-performance computing applications, Shifter and
Singularity have made significant steps toward remedying them.

##  6 Recommendations

We conclude by enumerating our recommendations for an open code policy:

  1. 1.

Require the _explicit_ use of an approved Open Source Initiative license (GPL,
MIT, Apache, etc.)

  2. 2.

Require the release of input files, analysis scripts, and build-time
parameters to ensure reproducibility

  3. 3.

Emphasize the importance of developing communities of practice around tools

  4. 4.

Consider the future of containerized distribution and encourage the
development of robust middleware to allow researchers to run analyses on
remote data sources provided as part of publication.

## References

  * Anders & Brown (2017) Anders, E. H., & Brown, B. P. 2017, Physical Review Fluids, 2, 083501, doi: 10.1103/PhysRevFluids.2.083501
  * Birn et al. (2001) Birn, J., Drake, J. F., Shay, M. A., et al. 2001, J. Geophys. Res., 106, 3715, doi: 10.1029/1999JA900449
  * Bordwell et al. (2018) Bordwell, B. R., Brown, B. P., & Oishi, J. S. 2018, ApJ, accepted
  * Clark & Oishi (2017a) Clark, S. E., & Oishi, J. S. 2017a, ApJ, 841, 1, doi: 10.3847/1538-4357/aa6ff1
  * Clark & Oishi (2017b) —. 2017b, ApJ, 841, 2, doi: 10.3847/1538-4357/aa6ff6
  * Couston et al. (2017) Couston, L.-A., Lecoanet, D., Favier, B., & Le Bars, M. 2017, Phys. Rev. Fluids, 2, 094804, doi: 10.1103/PhysRevFluids.2.094804
  * Davidovits & Fisch (2016a) Davidovits, S., & Fisch, N. J. 2016a, Phys. Rev. E, 94, 053206, doi: 10.1103/PhysRevE.94.053206
  * Davidovits & Fisch (2016b) —. 2016b, Physical Review Letters, 116, 105004, doi: 10.1103/PhysRevLett.116.105004
  * Kim et al. (2014) Kim, J.-h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14, doi: 10.1088/0067-0049/210/1/14
  * Lecoanet et al. (2014) Lecoanet, D., Brown, B. P., Zweibel, E. G., et al. 2014, ApJ, 797, 94, doi: 10.1088/0004-637X/797/2/94
  * Lecoanet et al. (2015) Lecoanet, D., Le Bars, M., Burns, K. J., et al. 2015, Phys. Rev. E, 91, 063016, doi: 10.1103/PhysRevE.91.063016
  * Lecoanet et al. (2017) Lecoanet, D., Vasil, G. M., Fuller, J., Cantiello, M., & Burns, K. J. 2017, MNRAS, 466, 2181, doi: 10.1093/mnras/stw3273
  * Lecoanet et al. (2016a) Lecoanet, D., Schwab, J., Quataert, E., et al. 2016a, ApJ, 832, 71, doi: 10.3847/0004-637X/832/1/71
  * Lecoanet et al. (2016b) Lecoanet, D., McCourt, M., Quataert, E., et al. 2016b, MNRAS, 455, 4274, doi: 10.1093/mnras/stv2564
  * Marti et al. (2014) Marti, P., Schaeffer, N., Hollerbach, R., et al. 2014, Geophysical Journal International, 197, 119, doi: 10.1093/gji/ggt518
  * Stodden et al. (2012) Stodden, V., Bailey, D. H., Borwein, J., et al. 2012, ICERM Report
  * Turk (2013) Turk, M. J. 2013, ArXiv e-prints.  https://arxiv.org/abs/1301.7064
  * Vasil et al. (2016) Vasil, G. M., Burns, K. J., Lecoanet, D., et al. 2016, Journal of Computational Physics, 325, 53, doi: 10.1016/j.jcp.2016.08.013