Morpheus: A Deep Learning Framework For Pixel-Level Analysis of Astronomical Image Data
Ryan Hausen Brant Robertson




Ryan Hausen
rhausen@ucsc.edu, brant@ucsc.edu

0000-0002-8543-761X]Ryan Hausen
Department of Computer Science and Engineering, University of
             California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064 USA

0000-0002-4271-0364]Brant E. Robertson
Department of Astronomy and Astrophysics, University of California,
             Santa Cruz, 1156 High Street, Santa Cruz, CA 95064 USA
Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ
             08540 USA




We present , a new model for generating pixel-level
morphological classifications of astronomical sources.  leverages
advances in deep learning to perform source detection, source segmentation, and
morphological classification pixel-by-pixel via a semantic segmentation
algorithm adopted from the field of computer vision. By utilizing morphological
information about the flux of real astronomical sources during object detection,
 shows resiliency to false-positive identifications of
sources. We evaluate  by performing source detection, source
segmentation, morphological classification on the Hubble Space Telescope
data in the five CANDELS fields with a focus on the GOODS South field,
and demonstrate a high completeness in recovering known GOODS South
3D-HST sources with H<26 AB. We release the code publicly, provide
online demonstrations, and present an interactive visualization of the
 results in GOODS South.













§ INTRODUCTION
 

Morphology represents the structural end state of the galaxy formation process.
Since at least <cit.>, astronomers have connected the morphological
character of galaxies to the physics governing their formation. Morphology can
reflect the initial conditions of galaxy formation, dissipation, cosmic
environment and large-scale tidal fields, merger and accretion history, internal
dynamics, star formation, the influence of supermassive black holes, and a range
of other physics
<cit.>.
The development of morphological measures for galaxies, therefore,
comprises an important task in observational astronomy. To help realize the
potential of current and future surveys for understanding galaxy formation
through morphology, this work presents , a deep learning-based model
for the simultaneous detection and morphological classification of objects
through the pixel-level semantic segmentation of large astronomical image
datasets.

The established connections between morphology and the physics of galaxy
formation run deep, and the way these connections manifest themselves
observationally depends on the measures of morphology used. Galaxy size
and surface brightness profile shape have served as common
proxies for morphology, as quantitatively measured from the light distribution
of objects <cit.>. Size, radial
profile, and isophotal shape or ellipticity vary with stellar mass and
luminosity
<cit.>.
When controlled for other variables, these measures of galaxy morphology may
show variations with cosmic environment
<cit.>,
redshift
<cit.>,
color <cit.>, star formation rate or quiescence
<cit.>,
internal dynamics <cit.>, the presence of active galactic nuclei
<cit.>, and stellar age
<cit.>. The presence and size of bulge, disk, and bar components
also vary with mass and redshift
<cit.>, and provide
information about the merger rate <cit.>.
Galaxy morphology encodes a rich spectrum of physical processes and
can augment what we learn from other galaxy properties.

While complex galaxy morphologies may be easily summarized with qualitative
descriptions (e.g., “disky”, “spheroidal”, “irregular”), providing
quantitative descriptions of this complexity represents a long-standing goal for
the field of galaxy formation and has motivated ingenuitive analysis methods
including measures of galaxy asymmetry, concentration, flux distribution
<cit.>,
shapelet decompositions <cit.>,
morphological principal component analyses <cit.>, and
unsupervised morphological hierarchical classifications <cit.>.
These measures provide well-defined characterizations of the surface brightness
distribution of galaxies and can be connected to their underlying physical state
by, e.g., calibration through numerical simulation <cit.>.
The complementarity between these quantitative measures and qualitative
morphological descriptions of galaxies means that developing both classes of
characterizations further can continue to improve our knowledge of galaxy
formation physics.

Characterizing large numbers of galaxies with descriptive classifications
simultaneously requires domain knowledge of galaxy morphology (“expertise”),
the capability to evaluate quickly each galaxy (“efficiency”), a capacity to
work on significant galaxy populations (“scalability”), some analysis of the
data to identify galaxy candidates for classification (“pre-processing”), a
presentation of galaxy images in a format that enables the characteristic
structures to be recognized (“data model”), and an output production of
reliable classifications (“accuracy”). Methods for the descriptive
classification of galaxy morphology have addressed these challenges in
complementary ways.

Perhaps the most important and influential framework for galaxy morphological
classification to date has been the Galaxy Zoo project
<cit.>, which enrolls the public in the
analysis of astronomical data including morphological classification. This
project has addressed the expertise challenge by training users in the
classification of galaxies and statistically accounting for the distribution of
users' accuracies. The efficiency of users varies, but by leveraging the power
of the public interest and enthusiasm, and now machine learning
<cit.>, the project can use scalability to offset
variability in the performance of individual users. The pre-processing and
delivery of suitable images to the users has required significant investment and
programming, but has led to a robust data model for both the astronomical data
and the data provided by user input. Science applications of Galaxy Zoo include
quantitative morphological descriptions of ∼50,000 galaxies
<cit.> in the CANDELS survey <cit.>,
probes of the connection between star formation rate and morphology in spiral
galaxies <cit.>, and measuring galaxy merger rates
<cit.>.

Other efforts have emphasized different dimensions of the morphological
classification task. <cit.> organized the visual classification
of ∼10,000 galaxies in CANDELS by a team of dozens of professional
astronomers. This important effort performed object detection and source
extraction on the CANDELS science data, assessed their completeness, and
provided detailed segmentation maps of the regions corresponding to classified
objects. The use of high expertise human classifiers leads to high accuracy, but
poses a challenge for scalability to larger samples. The work of
<cit.> also leveraged a significant investment in the
preprocessing and presentation of the data to their users with a custom
interface with a high-quality data model for the results.

Leveraging human classifiers, be they highly expert teams or well-calibrated
legions, to provide descriptive morphologies for forthcoming datasets will prove
challenging. These challenges motivate a consideration of other
approaches, and we present two salient examples in James Webb Space
Telescope <cit.> and the Large Synoptic Survey
Telescope <cit.>.

 enables both sensitive infrared imaging with NIRCam and
multiobject spectroscopy with NIRSpec free of atmospheric attenuation. The
galaxy population discovered by  will show a rich range of morphologies,
star formation histories, stellar masses, and angular sizes
<cit.>, which makes identifying NIRCam-selected samples for
spectroscopic follow-up with NIRSpec challenging. The efficiency gain of
parallel observations with NIRCam and NIRSpec will lead to programs
where the timescale for constructing NIRCam-selected samples will be very
short (∼2 months) to enable well-designed parallel survey geometries. For
this application, the ability to generate quick morphological classifications
for thousands of candidate sources will enhance the spectroscopic target
selection in valuable space-based observations.

 presents a challenge of scale, with an estimated 30 billion astronomical
sources, including billions of galaxies over ∼17,000 ^2
<cit.>. The morphological classification of these galaxies will
require the development of significant analysis methods that can both scale to
the enormity of the  dataset and perform well enough to allow imaging
data to be reprocessed in pace with the  data releases. Indeed,
morphological classification methods have been identified as keystone
preparatory science tasks by in the  Galaxies Science Roadmap
<cit.>.

Recently, advances in the field of machine learning called deep learning have
enjoyed success in morphological classification. <cit.> (D15) and
<cit.> use deep learning to classify the Galaxy Zoo Survey.
<cit.> used a deep learning model derived from D15 and the
classifications from K15 to classify the CANDELS survey. <cit.>
used deep learning to perform galaxy detection and morphological classification,
an approach that has also been used to characterize Dark Energy Survey galaxy
morphologies <cit.>. Deep learning models have been further
applied to infer the surface brightness profiles of galaxies
<cit.> and measure their fluxes  <cit.>, and now to
simulate entire surveys <cit.>.

Here, we extend previous efforts by applying a semantic segmentation algorithm
to both classify pixels and identify objects in astronomical images using our
deep learning framework called . The software architecture of the
 framework is described in Section <ref>, with the
essential convolutional neural network and deep learning components reviewed in
Appendix <ref>. The  framework has been engineered by
using TensorFlow <cit.> implementations of these components to
perform convolutions and tensorial operations, and is not a port of existing
deep learning frameworks or generated via “transfer
learning”<cit.> of existing frameworks pre-trained on
non-astronomical data such as ImageNet <cit.>.

We train  using multiband Flexible Image Transport System
<cit.> images of CANDELS galaxies visually classified by
<cit.> and their segmentation maps derived from standard sextractor analyses <cit.>. The training procedure is described in
Section <ref>, including the “loss function” used to optimize
the  framework. Since  provides local estimates of whether
image pixels contain source flux, the  output can be used to perform
source segmentation and deblending. We present fiducial segmentation and
deblending algorithms for  in Section <ref>.

We then apply  to the Hubble Legacy Fields <cit.>
reduction of the CANDELS and GOODS data in the GOODS South region, the
v1.0 data release <cit.> for the other four CANDELS
regions, and generate FITS data files of the same pixel format as the input FITS
images, each containing the pixel-by-pixel model classifications of the image
data into spheroid, disk, irregular,
point source/compact, and background classes, as described in
Section <ref>. We release publicly these 
pixel-level classification data products and detailed them in Appendix
<ref>. We evaluate the performance of  in
Section <ref>, including tests that use the catalog of
3D-HST photometric sources <cit.> to measure the
completeness of  in recovering sources as a function of source
magnitude. We find that  is highly complete (>90%) for sources
up to one magnitude fainter than objects used to train the model. Using the
 results, we provide estimates of the morphological classification of
3D-HST sources as a public value-added catalog, described in Section
<ref>. In Section <ref>, we discuss applications
of  and semantic segmentation, which extend well beyond morphological
classification, and connect the capabilities of  to other research
areas in astronomical data analysis. We publicly release the  code,
provide on-line tutorials for using the framework via Jupyter notebooks, and
present an interactive website to visualize the  classifications and
segmentation maps in the context of the HLF images and 3D-HST catalog. These
software and data releases are described in Appendices
<ref>, <ref>, and
<ref>. A summary of our work is presented with some
conclusions in Section <ref>. Throughout the paper, we have used
the AB magnitude system <cit.> and assumed a flat ΛCDM
universe (Ω_m=0.3, Ω_Λ=0.7) with a Hubble parameter H_0=70
km/s/Mpc when necessary.







§  DEEP LEARNING FRAMEWORK







 provides a deep learning framework for analyzing astronomical images
at the pixel level. Using a semantic segmentation algorithm, 
identifies which pixels in an image are likely to contain source flux and
separates them from “background” or sky pixels. ,
therefore, allows for the definition of corresponding segmentation
regions or “segmentation maps” by finding contiguous regions of source pixels
distinct from the sky. Within the same framework,  enables for
further classification of the source pixels into additional “classes”. In this
paper, we have trained  to classify the source pixels into
morphological categories (spheroid, disk, irregular, point source/compact, and background) approximating the visual
classifications performed by the CANDELS collaboration in K15. These source
pixel classes identified by  could, in principle, be
trained to reproduce other properties of the galaxies, such as, e.g.,
photometric redshift, provided a sufficient training dataset is available. In
the sections below, we describe the architecture of the  deep
learning framework. Readers unfamiliar with the primary computational elements
of deep learning architectures may refer to Appendix <ref> where
more details are provided.



 §.§ Input Data



We engineered the  deep learning framework to accept astronomical
image data as direct input for pixel-level analysis.  operates on
science-quality FITS images, with sufficient pipeline processing (e.g., flat
fielding, background subtraction, etc.) to enable photometric analysis.
 accepts multiband imaging data, with a FITS file for each of the
n_b bands used to train the model (see Section <ref>). The
pixel format of the input FITS images (or image region) matches the format of
FITS images used to perform training, reflecting the size of the convolutional
layers of the neural network determined before training.  allows for
arbitrarily large images to be analyzed by subdividing them into regions that
the model processes in parallel, as described in Section
<ref> below.

For the example application of morphological classification presented in this
paper, we use the F606W(V), F850LP(z), F125W(J), and F160W(H)
band images from
Hubble Space Telescope for training, testing, and our final analysis. Our
training and testing images were FITS thumbnails and segmentation maps provided
by <cit.>. Once trained,  can be applied to
arbitrarily large images via a parallelization scheme described below in Section
<ref>. We have used the CANDELS public release data
<cit.> in additional performance tests and the
Hubble Legacy Fields v2.0 data <cit.> for our  data
release.

We note that the approach taken by  differs from deep learning models
that use traditional image formats, e.g., three-color Portable
Network Graphics (PNG) or Joint Photographic Experts Group (JPEG) images as
input. Using PNG or JPEG files as input is convenient because deep learning
models trained on existing PNG or JPEG datasets, such as ImageNet
<cit.>, can be retrained via transfer learning  to
classify galaxies. However, the use of these inputs requires additional
pre-processing beyond the science pipeline, including arbitrary decisions
about how to weight the FITS images to represent the channels of the multicolor
PNG or JPEG. With the goal of including  framework analyses as part
of astronomical pipelines, we have instead used FITS images directly as input to
the neural network.




 §.§ Neural Network



 uses a neural network inspired by the U-Net architecture
<cit.>
and is implemented using Python 3 <cit.> and the TensorFlow library
<cit.>. We construct  from a series of “blocks”
that combine multiple operations used repeatedly by the model.
Each block performs a sequence of “block operations”. Figure <ref>
provides an illustration of a  block and its block operations. Block
operations are parameterized by the number Q of convolved output images, or
“feature maps”, they produce, one for each convolutional artificial neuron in
the layer. We describe this process in more detail below.

Consider input data 𝐗, consisting of K layers of images with
N× M pixels. We define a block operation on 𝐗 as


    OP_Q(𝐗) = ReLU(CONV_Q(BN(𝐗)),

where ReLU is the Rectified Linear Unit activation function <cit.>, CONV_Q is a
convolutional layer (see Appendix <ref>) with a number Q
convolutional artificial neurons (see Appendix <ref>), and
BN is the batch normalization procedure <cit.>. Note that the values of Q
appearing in OP_Q and CONV_Q are equal. For example,
OP_4 would indicate that the convolutional layer within the
OP_4 function has 4 convolutional artificial neurons. Unless stated
otherwise, all inputs into a convolutional layer are zero-padded to
preserve the width and height of the input, and all convolutional
artificial neurons have kernel dimensions 3 × 3. Given Equation
<ref>, for an input 𝐗 with dimensions N × M ×
K the output of the function OP_4(𝐗) would have dimensions N
× M × 4.

Equation <ref> allows for a recursive definition of a
function describing a series of block operations, where the input data to one
block operation consist of the output from a previous block operation. This
recursion can be written as

    OP^P_Q(𝐗) =
        𝐗,     if  P = 0 
    ReLU(CONV_Q(BN(OP_Q^P-1(𝐗)))     if  P > 0
        .

Equation <ref> introduces a new parameter P, shown with
a superscript in OP_Q^P. The parameter P establishes the conditions
of a base case for the recursion. Note that in Equation
<ref> the input 𝐗 is processed directly when
P=1, and when P>1 the input to the OP_Q^P function is the output
from OP_Q^P-1. It can be seen from the formulation of Equations
<ref> and <ref> that OP_Q(𝐗) =
OP_Q^1(𝐗).

Since a block performs a number P block operations, we
can define a block mathematically as

    BLOCK(Q, P, 𝐗) = OP^P_Q(𝐗).

An example block and its block operations can be seen diagrammatically in Figure
<ref>. With these definitions, we can present the neural network
architecture used in .

Like the U-Net architecture, the  architecture consists of a
contraction phase and an expansion phase. The contraction phase consists of
three blocks with parameters (P=4, Q=8), (P=4, Q=16), and (P=4, Q=32).
Each block is followed by a max-pooling operation with size=(2 ×
2) (see Section <ref>), halving the width and height of its
input. After the contraction phase there is a single intermediary block
preceding the expansion phase with the parameters (P=1, Q=16). The expansion
phase consists of three blocks with the parameters (P=2, Q=8), (P=2, Q=16),
(P=2, Q=32). Each block is preceded by a bicubic interpolation operation that
doubles the width and the height of its input. Importantly, the output from each
block in the contraction phase is concatenated (see Section
<ref>) with the output from the bicubic interpolation
operation in the expansion phase whose output matches its width and height (see
Figure <ref>). The output from the final block in the expansion
phase is passed through a single convolutional layer with 5 convolutional
artificial neurons. A softmax operation (see Equation <ref>) is
performed on the values in each pixel, ensuring the values sum to unity. The
final output is a matrix with the same width and height as the input into the
network, but where the last dimension, 5, now represents a classification
distribution describing the confidence the corresponding pixel from the input
belongs to one of the 5 specified morphological classes.




The blocks in  are organized into the U-Net structure, shown in
Figure <ref>. The model proceeds clockwise, starting from
“Input” on the upper left through to “Output” on the lower left. The very
first step involves the insertion of the input FITS images into the model. Each
FITS image is normalized to have a mean of 0 and unit variance before processing
by . We will refer to the number of input bands as n_b, and in the
application presented here, we take n_b=4 (i.e., VzJH). The input
images each have pixel dimensions N× M, and we can,
therefore, consider the astronomical input data to have dimensions
N× M × n_b. Only the first block operation takes the FITS images as
input, and every subsequent block operation in the model takes the
output from previous blocks as input.

The first convolution in the first block operation convolves the normalized
N× M × n_b astronomical data with three-dimensional kernels of size
n_k^2 × n_b, and each element of the kernel is a variable parameter of
the model to be optimized. The convolutions operate only in the two pixel
dimensions, such that n_b convolutions are performed, one for each N× M
pixel image, using a different n_k × n_k kernel for each convolution. The
n_b convolved images are then summed pixel by pixel to create an output
feature map of size N× M. The convolutional layer repeats this process
Q times with different kernels, generating Q output feature maps and an
output dataset of size N × M × Q. For the first block in 
we use Q=8 (see Figure <ref>). After the first convolution on
the astronomical data, every subsequent convolution in the first block has both
input and output data of size N × M × Q.

Each block performs a number P block operations, resulting in output data with
dimensions of N× M × Q emerging from the block. The number of
feature maps Q changes with each block. For a block producing Q filters, if
the data incoming into the block has size N× M × Q' with Q' Q,
then the first convolutional layer in the first block operation will have Q
kernels of size n_k^2 × Q'. All subsequent convolutional layers in the
block will then ingest and produce data of size N× M × Q by using
kernels of size n_k^2 × Q.

We can apply further operations on the data in between the blocks, and the
character of these operations can affect the dimensions of the data. The first
half of the model is a contraction phase, where each block is followed by a
max-pooling operation <cit.>. The max-pooling is applied to each
feature map output by the block, taking the local maximum over small areas
within each feature map (in the version of  presented here, a
2×2 pixel region) and reducing the size of the data input to the next
block by the same factor. For this paper, the contraction phase in the
 framework uses three pairs of blocks plus max-pooling layers.

After the contraction phase, the model uses a series of blocks,
bicubic interpolation layers, and data concatenations in an expansion phase to
grow the data back to the original format. Following each block in the expansion
phase, a bicubic interpolation layer expands the feature maps by the same areal
factor as the max-pooling layers applied in the contraction phase
(2×2 in the version of  presented here). The output feature
maps from the interpolation layers are concatenated with the output feature maps
from the contraction phase blocks where the data have the same format. Finally,
the output from the last block in the expansion phase is input into a
convolutional layer that produces the final output images that we call
“ classification images”, one image for each class. The pixel
values in these images contain the model estimates for their classification,
normalized such that the element-wise sum of the classification images equals
unity.  For this paper, where we are performing galaxy morphological
classification, there are five classification images (spheroid, disk, irregular, point source /compact, and background).

As the data progresses through the model, the number of feature maps and their
shapes change owing to the max-pooling and interpolation layers. For
reference, in Table <ref>, we list the dimensions of
the data at each stage in the model, assuming input images in n_b
bands, each with N× M pixels, and a total of n_c classification images
produced by the model.









 §.§ Parallelization for Large Images



While the  neural network performs semantic segmentation on pixels in
FITS images with a size determined by the training images, the model can process
and classify pixels in arbitrarily large images. To process large
images,  uses a sliding window strategy by breaking the
input FITS files into thumbnails of size N× M (the size of the training
images) and classifying them individually.  proceeds through the
large format image, first column by column, and then row by row, shifting
the active N× M window by a unit pixel stride and then recomputing the
classification for each pixel.

As the classification process continues with unit pixel shifts, each
 pixel is deliberately classified many times. We noticed
heuristically that the output  classification of pixels depended on
their location within the image, and that the pixel classifications were more
accurate relative to our training data when they resided in the inner n_p =
(N-B)×(M-B) region of the classification area, where the lesser accuracy
region consisted of a border about B∼5 pixels wide on each side. Outside of
the very outer B pixels in the large format image,  classifies each
pixel n_p times. For the large FITS data images used in this paper, this
repetition corresponds to n_p = 900 separate classifications per pixel per
output class, where each classification occurs when the pixel lies at a
different location within the active window. This substantial additional
information can be leveraged to improve the model, but storing the full
“distribution” of classifications produced by this method would increase our
data volume by roughly three orders of magnitude.

While  would enable  full use of these distributions,
for practical considerations, we instead record some statistical
information as the computation proceeds and do not store the entire set of n_p
samples. To avoid storing the full distribution, we track running estimates of
the mean and variance of the distribution[See, e.g.,
http://people.ds.cam.ac.uk/fanf2/hermes/doc/antiforgery/stats.pdf for an example
of running mean and variance estimators.]. Once the mean for each class for each
pixel is computed, we normalize the means across classes to sum to unity. We
further record a statistic we call rank voting, which is a tally of the
number of times each output class was computed by the model to be the top class
for each pixel. The sum of rank votes across classes for a single pixel equals
the number of times  processed the pixels (i.e., n_p for most
pixels). After the computation, the rank votes are normalized to sum to unity
across the classes for each pixel.

The strips of classified regions produce fifteen output images, containing the
mean and variance estimators for the classification distribution and normalized
rank votes for each class. This striped processing of the image can be performed
in parallel across multiple  instances and then stitched back
together. The weak scaling of this processing is, in principle,
trivial and is limited only by the number of available GPUs and the
total memory of the computer used to perform the calculation.







§ MODEL TRAINING



The training of deep learning frameworks involves important decisions
about the training data, the metrics used to optimize the network, numerical
parameters of the model, and the length of training. We provide some rationale
for these choices below.





 §.§ Training Data



To train a model to perform semantic segmentation, we require a dataset that
provides both information on the segmentation of regions of interest and
classifications associated with those regions. For galaxy morphological
classification, we use 7,629 galaxies sampled from the K15 dataset. Their
2-epoch CANDELS data provide an excellent combination of multiband FITS
thumbnails, segmentation maps in FITS format, and visually-classified
morphologies in tabulated form. The K15 classifications consisted of votes by
expert astronomers, between 3-60 per object, who inspected images of galaxies
and then selected from several morphological categories to
assign to the object. The number of votes for each category for each object are
provided, allowing  to use the distribution of votes across
classifications for each object when training. We downloaded and used the
publicly available K15 thumbnail FITS files for the F606W, F850LP, F125W,
and F160W bands as input into the model for training and testing. In
training Morpheus to reproduce the K15 classifications, multiband data approximates
the information provided to the astronomers who performed
the K15 classifications. Morpheus is trained using the same V, z, J, and
H-band image thumbnails used in the K15 classification process. Other bands
or different numbers of bands could be used for training as necessary, and
 allows for reconfiguration and retraining depending on the available
training images. Of the K15 data set, we used 80% of the objects to form our
training sample and 20% to form our test sample. Various statistical
properties of the test and training samples are described throughout the rest of
the paper.

The primary K15 classifications spheroid, disk, irregular, and
point source/compact were used in the example 
application presented here. We added one additional class, background,
to represent sky pixels absent significant source flux. We classify pixels as
belonging to the background category if those pixels fell outside the K15
segmentation maps. Pixels inside the segmentation maps were assigned the
distribution of classifications provided by the K15 experts.

The K15 classification scheme also included an unknown class for
objects. Since  works at the pixel level and could provide individual
pixel classifications that were locally accurate within a source but that
collectively could sum to an object whose morphology expert astronomers might
classify as unknown, we were posed with the challenge of how to treat
the K15 unknown class. Given our addition of the background
class constructed from large image regions dominated by sky, one might expect
overlap in the features of regions that are mostly noise and amorphous regions
classified as unknown. Since one might also expect overlap between
unknown and irregular classifications, we wanted to preserve
some distinction in the object classes. We, therefore, removed the
unknown class by removing any sources that had unknown as
their primary classification from the training sample (213 sources). For any
sources where the non-dominant K15 classifications included unknown, we
redistributed the unknown votes proportionally to the other classes.







 §.§ Data Augmentation



To increase the effective size of the training dataset,  uses a data
augmentation method. Augmentation supplements the input training data set by
performing transformations on the training images to alter them with the intent
of adding similar but not identical images with known classifications.
Augmentation has been used successfully in the context of galaxy morphological
classification <cit.>, and  adopts a
comparable approach to previous implementations.

During training,  produces a series of 40×40 pixel augmented
versions of the training images. The augmentation approach is illustrated in
Figure <ref>. For each band in the original training image, the image
is collectively rotated by a random angle ϕ∈[0,2π], flipped
horizontally with a random 50% probability, and then flipped vertically with
a random 50% probability. A crop of the inner 60× 60 pixels of the
resulting image is produced, and then a random 40× 40 pixel subset of the
image is selected and passed to the model for training. This method allows us to
increase the effective number of images available for training by a factor of
∼574,400 and helps ameliorate over-training on the original
training image set.






 §.§ Loss Function



A standard method for training deep learning frameworks is to define a loss
function that provides a statistic based on the output classifications to
optimize via stochastic gradient descent with gradients computed using
back-propagation <cit.>. Here, we describe how the
 loss function is constructed.

The first task is to assign a distribution of input classifications on a
per-pixel basis, choosing between the n_c classes available to the
 model. For this work, we choose n_c=5 (background, disk, spheroid, irregular, and point source/compact), but
 can adopt an arbitrary number of classes. We use the index k to
indicate a given class, with k∈[1,n_c]. Consider an N × M image of an
astronomical object that has been visually classified by a collection of
experts, and a segmentation map defining the extent of the object in the image.
Outside the segmentation map of the object, the pixels are assumed to belong to
the sky and are assigned the background class. Inside the segmentation
map, pixels are assigned the distribution of disk, spheroid, irregular, and point source/compact classifications determined by the
experts for the entire object. For each pixel ij, with i∈[1,N] rows and
j∈[1,M] columns, we then have the vector _ij whose elements q_ijk
contain the input distribution of classifications. Here, the index k runs over
the number of classes n_c and ∑_k q_ijk = 1 for each pixel with indices
ij. The goal of the model is to reproduce this normalized distribution
_ij of discrete classes for each pixel of the training images. We wish to
define a total loss function L_tot that provides a single per-image
statistic for the model to optimize when attempting to reproduce _ij.
 combines a weighted cross entropy loss function with a Dice
loss <cit.> for its optimization statistic,
which we describe below.


At the end of the  data flow, as outlined in Figure
<ref>, the raw output of the model consists of N× M vectors
_ij with n_c elements per-pixel estimates that represent unnormalized
approximations to the input per-pixel distributions _ij. The model
outputs _ij for each pixel are then normalized to form a probability
distribution _ij using the softmax function

    _ij=exp(_ij)/∑_k=1^n_cexp(x_ijk), for  k∈[1,n_c].

The distribution _ij then represents the pixel-by-pixel classifications
computed by  for each of the k∈[1,n_c] classes. For a pixel with
indices ij, we can define the per-pixel cross entropy loss function as

    L_ij(_ij,_ij)=-∑_k=1^n_c p_ijklog(q_ijk)

where _ij and _ij are again the two per-pixel probability
distributions, with _ij representing the true distribution of the input
classifications for the pixel ij and _ij representing the model output.

Equation <ref> provides the per-pixel contribution to the
entropy loss function. However, for many images, the majority of pixels
lie outside the segmentation maps of sources identified in the training data and
are therefore labeled as background. To overcome this imbalance and
disincentivize the model from erroneously learning to classify pixels containing
source flux as background, we apply a weighting to the per-pixel loss. We
define an index k_ij^max, that indicates which class is the maximum of
the input classification distribution for each pixel, written as

    k_ij^max,  = argmax _ij

with 1≤ k_ij^max, ≤ n_c. For each class k, we then define a
weight w_k that is inversely proportional to the number of pixels with
k_ij^max,  = k. We can write

    w_k =[∑_i=1^N ∑_j=1^M max(_ij) δ_k,k_ij^max, ]^-1.

Here, δ_i,j is the Kronecker delta function. The vector  has size
n_c and each of its elements w_k contain the inverse of the sum of
max(_ij) for pixels with k_ij^max, =k. In a given image, we
ignore any classes that do not appear in the input classification distribution
(i.e., any class k for which ∑_i∑_j q_ijk = 0).

Using , we define a weighted cross entropy loss for each pixel as

    L^w_ij = w_k_ij^max L_ij(_ij,_ij).

A mean weighted loss function is then computed by averaging Equation
<ref> over all pixels as

    L̅^w = 1/N× M∑_i=1^N∑_j=1^M L^w_ij.

This mean weighted loss function serves as a summary statistic of the cross
entropy between the output of  and the input classification
distribution.

When segmenting images primarily comprised of background pixels, the
classification distributions of the output pixels should be highly
unbalanced, with the majority having background≈1. In this
case, the mean loss function statistic defined by Equation
<ref> will be strongly influenced by a single class. A
common approach to handle unbalanced segmentations is to employ a Dice loss
function to supplement the entropy loss function
<cit.>. The Dice loss function used by
 is written as

    L^D(b, m)=1 - 2∑_i∑_j(S(b)
                          ∘m)_ij/∑_i∑_j(S(b) + m)_ij.

Here, S(b)=(1+exp(-b))^-1 is the sigmoid function (see Equation
<ref>) applied pixel-wise to the background classification
image output by the model. The image m is the input mask with values
m=1 denoting background pixels and m=0 indicating source pixels,
defined, e.g., by a segmentation map generated using sextractor. The
∘ symbol indicates a Hadamard product of the matrices S(b) and
m. Note that the output background matrix b has not yet been
normalized using a softmax function, and so b_ij∈[-∞,∞] and
S(b_ij)∈[0,1]. The Dice loss then ranges from L^D=0 if
S(b)≈m and L^D∼1 when S(b) and m differ
substantially. The addition of this loss function helps to maximize the spatial
coincidence of the output background pixels assigned b_ij≈1
with the non-zero elements of the input segmentation mask m.

To define the total loss function optimized during the training of ,
the cross entropy and Dice losses are combined as a sum weighted by two
parameters λ_w and λ_D. The total loss function is written
as

    L_tot = λ_wL^w + λ_DL^D

For the implementation of  used in this paper, the entropy and Dice
loss functions are weighted equally by setting λ_w =1 and
λ_D=1.




 §.§ Optimization Method


To optimize the model parameters, the Adam stochastic gradient descent
method <cit.> was used. The Adam algorithm uses the first and
second moments of first-order gradients computed via backpropagation to find
the minimum of a stochastic function (in this case, our loss function,
see Section <ref>, which depends on the many parameters of
the neural network). The Adam optimizer, in turn, depends on
hyper-parameters that determine how the algorithm iteratively finds a minimum.
Since the loss function is stochastic, the gradients change each
iteration, and Adam uses an exponential moving average of the
gradients (m̂) and squared gradients (v̂) when searching for a
minimum. Two dimensionless hyper-parameters (β_1 and β_2) set the
decay rates of these exponential averages <cit.>. As the parameters θ of the function being optimized are
iterated between steps t-1 and t, they are updated according to

    θ_t ←θ_t-1 - α·m̂_t / (√(v̂_t) + ϵ).

Here, ϵ is a small, dimensionless safety hyper-parameter that prevents
division by zero, and α is a small, dimensionless hyper-parameter that
determines the magnitude of the iteration step. Table <ref>
lists the numerical values of the Adam optimizer hyper-parameters used by
. We use the default suggested values for β_1, β_2, and
ϵ. After some experimentation, we adopted a more conservative step size
for α than used by <cit.>.






 §.§ Model Evaluation





As training proceeds, the performance of the model can be quantified using
various metrics and monitored to determine when training has effectively
completed. The actual performance of  will vary depending on the
classification scheme used, and here we report the performance of the model
relative to the CANDELS images morphologically classified in K15. Performance
metrics reported in this Section refer to pixel-level quantities, and we discuss
object-level comparisons of morphological classifications relative to K15 in
Section <ref>.

While the model training proceeds by optimizing the loss function defined in
Section <ref>, we want to quantify the accuracy of the
model in recovering the per-pixel classification and the overlap of contiguous
regions with the same classification. First, we will need to define the index
k_ij^max with maximum probability to reflect either the input
classification _ij or the output classification _ij. We define an
equivalent of Equation <ref> for _ij as

    k_ij^max,  = argmax _ij.

We can then define a percentage accuracy

    = 100/N × M∑_i=1^N∑_j=1^Mδ_k_ij^max, , k_ij^max, .

The accuracy  then provides the percentage of pixels for which the maximum
probability classes of the input and output distributions match.

In addition to accuracy, the intersection-over-union  of pixels with background probabilities above some threshold is computed between the input
_ij and output _ij distributions. If we define the index b to
represent the background class, we can express the input background
probabilities as _b=q_ijb for i∈[1,N] and j∈[1,M], and
the equivalent for the output background probabilities _b. We can
refer to _b and _b as the input and output background images,
and the regions of these images with values above some threshold  B as
_b(>B) and _b(>B), respectively. Note that the input _b only
contains values of zero or one, whereas the output _b has continuous values
between zero and one. We can then define the  metric for threshold B as

    (B) = _b(>B)∩_b(>B)/_b(>B)∪_b(>B).

Intuitively, this  metric describes how well the pixels assigned by
 as belonging to a source match up with the input source segmentation
maps. A value of =1 indicates a perfect match between source pixels
identified by  and the input segmentation maps, while a value of
=0 would indicate no pixels in common between the two sets.

As training proceeds, the accuracy  and intersection-over-union  are
monitored until they plateau with small variations. For the K15 training data,
the model plateaued after about 400 epochs. The training then continues for
another 100 epochs to find a local maximum in  and , and the model
parameters at this local maximum adopted for testing. Table
<ref> summarizes the per-pixel performance of  in
terms of  for each class separately,  for all classes, and (B)
for B=[0.5,0.6,0.7,0.8,0.9]. We also report the performance
of the training and testing samples separately. The
pixel-level classifications are 70-90% accurate depending on the class, and
the intersection-over-union is ∼0.9 for all thresholds B≥0.5.
The model shows some evidence for overfitting as accuracy declines slightly from
the training to test sets for most classes.







§ SEGMENTATION AND DEBLENDING



To evaluate the completeness of  in object detection and to compute
an object-level classification, segmentation maps must be constructed and then
deblended from the  pixel-level output.  uses the background class from the output of the neural network described in Section
<ref> to create a segmentation map. The segmentation algorithm uses a
watershed transform to separate background pixels from source
pixels and then assigns contiguous source pixels a unique label.
The deblending algorithm uses the flux from the input science images and the
output of the segmentation algorithm to deblend under-segmented regions
containing multiple sources. We summarize these procedures as Algorithms
<ref> and <ref>. Figure <ref>
illustrates the process for generating and deblending segmentation maps.



 §.§ Segmentation



The segmentation algorithm operates on the output background
classification image and identifies contiguous regions of low
background as sources. The algorithm begins with the background image
≡_b defined in Section <ref> and an initially
empty mask =0 of the same size. For every pixel in the image,
if b_ij=1 we set m_ij=1 and if b_ij=0 we set m_ij=2. The background mask  then indicates extreme regions of . The
<cit.> algorithm is applied to the background image  to
produce a Sobel edge image .  then applies the watershed
algorithm of <cit.>, using the Sobel image  as the “input
image” and the background mask  as the “marker set”. We refer
the reader to <cit.> for more details on the watershed algorithm,
but in summary, the watershed algorithm collects together regions with
the same marker set value within basins in the input image. The Sobel image
 provide these basins by identifying edges in the
background, and the background mask  provides the marker
locations for generating the individual sheds. The output of the watershed
algorithm is then an image  containing distinct regions generated from
areas of low background that are bounded by edges where the background is changing quickly. The algorithm then visits each of the distinct
regions in  and assigns them a unique id, creating the segmentation map
 before deblending.







 §.§ Deblending



The algorithm described in Section <ref> provides a
collection of segmented regions of contiguous areas, each with a unique index.
Since this algorithm identifies contiguous regions of low background,
neighboring sources with overlapping flux in the science images will be blended
by the segmentation algorithm. The deblending algorithm used in  is
ad hoc and is primarily designed to separate the segmented
regions into distinct subregions containing a single pre-defined object. The
locations of these objects may be externally specified, such as catalog entries
from a source catalog (e.g., 3D-HST sources), or they may be internally derived
from the science images themselves (e.g., local flux maxima).

The deblending algorithm we use applies another round of the watershed operation
on each of the distinct regions identified by the segmentation algorithm, using
the local flux distributions from the negative of a science image (e.g.,
F160W) as the basins to fill and object locations as the marker set. We assign
the resulting subdivided segmentations a distinct subid in addition to their
shared id, allowing us to keep track of adjacent deblended regions that share
the same parent segmentation region. The subid of deblended sources is
indicated by decimal values and the parent id is indicated by the whole number
of the id. For example, if a source with id=8 was actually two sources,
after deblending the two deblended sources would have id values 8.1 and 8.2.

In testing , we find that the deblending algorithm may shred
extended sources like large disks or point source diffraction spikes. However,
the Morpheus algorithm successfully deblends some small or faint sources
proximate to bright sources that are missing from the 3D-HST catalog.










§ OBJECT-LEVEL CLASSIFICATION



While  uses a semantic segmentation model to enable pixel-level
classification of astronomical images using a deep learning framework, some
applications, like the morphological classification of galaxies, additionally
require object-level classification.  aggregates pixel-level
classifications into an object-level classification by using a flux-weighted
average.




Figure <ref> shows the results of the  pixel-level
classification for an example area of the CANDELS region of GOODS South. The
leftmost panel shows a three-color VzH
composite of the example area for reference, though  operates
directly on the science-quality VzJH FITS images. The central panels show the
output pixel classifications (i.e.,  from Section
<ref>) for the background, spheroid, disk, irregular, and point source/compact classes, with the
intensity of each pixel indicating the normalized probability
q_ijk∈[0,1]. The segmentation map resulting from the algorithms
described in Section <ref> is also shown in as a central
panel. The rightmost panel shows a color composite
of the  pixel-level classification, with the color of each pixel
indicating its the dominant class and the saturation of the pixel being
proportional to the difference Δ between the dominant and second most
dominant class. White pixels then indicate regions where the model did not
strongly distinguish between two classes, such as in transition regions in the
image between two objects with different morphological classes. The pixel
intensities in the pixel-level classification image are set to 1-background and are not flux-weighted. The dominant classification
for each object, as determined by , is often clear
visually. The brightest objects are well-classified and agree with the intuitive
morphological classifications an astronomer might assign based on the VzH
color composite image. Faint objects in the image have less morphological
information available and are typically classified as
point source/compact, in rough agreement with their classifications in the
K15 training set. However, these visual comparisons are qualitative, and we now
turn to quantifying the object-level classification from the pixel values.

Consider a deblended object  containing a total of n_o contiguous pixels
of arbitrary shape within a flux image, and a single index i=[1,n_o]
scanning through the pixels in . Each class k∈[1,n_c] in the
distribution of classification probabilities Q for the object is computed
as

    Q_k = ∑_i=1^n_o y_i q_ik/∑_i=1^n_o y_i.

Here,  represents the pixel region in a science image assigned to the
object, and y_i is the flux in the ith pixel of the object. The quantity
q_ik is the kth classification probability of the ith pixel in .
Equation <ref> represents object-level classification
computed as the flux-weighted average of the pixel-level classifications in the
object.








§  DATA PRODUCTS



Before turning the quantifications of the object-level performance, we provide a
brief overview of the derived data products produced by . A more
detailed description of the data products is presented in Appendix
<ref>, where we describe a release of pixel-level
morphologies for the 5 CANDELS fields and
3D-HST value-added catalog, including object-level
morphologies. The Hubble Legacy Fields <cit.> GOODS South
v2.0 release and 3D-HST survey <cit.> are the primary focus of the
analysis of the ' performance owing to their depth and completeness.

As described in Section <ref>,  produces a
set of n_c “classification images” that correspond to the pixel-by-pixel
model estimates _ij for each class, normalized across classes such that
∑_k q_ijk = 1. The value of each pixel is, therefore,
bounded (q_ijk∈[0,1]). The classification images are stored in FITS
format, and inherit the same (N× M) pixel dimensions as the input
FITS science images provided to . When presenting classification
images used in this paper, we represent background images in negative
grayscale, spheroid images in black-red, disk images in black-blue, irregular images in black-green, and point source/compact images in black-yellow color scales. Figure
<ref> shows spheroid, disk, irregular, point
source/compact, and background images (central panels) for a region of
CANDELS GOODS South.

Given the separate classification images, we can construct what we deem a
“ morphological color image" that indicates the local dominant class
for each pixel. To produce a Red-Blue-Green false color image to represent the
morphological classes visually, we use the Hue-Saturation-Value (HSV) color
space and convert from HSV to RGB via standard conversions. In the HSV color
space, the Hue image indicates a hue on the color wheel, Saturation provides the
richness of the color (from white or black to a deep color), and Value sets the
brightness of a pixel (from dark to bright). On a color wheel of hues,
ℋ∈[0,360] ranges from red (ℋ=0) to red
(ℋ=360) through yellow (ℋ=120), green
(ℋ=180), and blue (ℋ=240), we can assign Hue pixel
values corresponding to the dominant morphological class (spheroid as red,
disk as blue, irregular as green, and point source/compact
as yellow). We set the Saturation of the image to be the Δ q_ijk
between the dominant class and the second most dominant class, such that cleanly
classified pixels (q_ijk_ij^max≈1, Δ q_ijk≈1)
appear as deep red, blue, green, or yellow, and pixels where 
produces an indeterminate classification (Δ q_ijk≈0) appear as
white or desaturated. The Value channel is set equal to 1 -background,
such that regions of low background containing sources are bright, and
regions with high background are dark. Figure <ref> also shows the
 morphological color image (far right panel) for a region of CANDELS
GOODS South.




 §.§ Morphological Images for GOODS South



As part of our data products, we have produced  morphological images
of the Hubble Legacy Fields <cit.> reduction of
GOODS South. These data products are used in Section <ref>
to quantify the performance of  relative to standard astronomical
analyses, and we, therefore, introduce them here. The 
morphological classification images for the HLF were computed as described in
Section <ref>, feeding  subregions of the HLF
VzJH images for processing and then tracking the distribution of output pixel
classifications to select the best classification for each. The ∼10^8
pixels in each classification image are then stitched back together to produce
contiguous background, spheroid, disk, irregular, and
point source/compact images for the entire HLF GOODS South.





  §.§.§ Background Image



Figure <ref> shows the background image for the 
analysis of the HLF reduction of GOODS South. The background
classification for each pixel is shown in negative gray scale, with black
corresponding to background=1 and white regions corresponding to background=0. The background image is used throughout Section
<ref> to quantify the performance of  in object
detection.








  §.§.§ Spheroid Image



Figure <ref> shows the spheroid image for the 
analysis of the HLF reduction of GOODS South. The spheroid classification
for each pixel is shown on a black-to-red colormap, with black corresponding
to spheroid=0 and red regions corresponding to spheroid=1.








  §.§.§ Disk Image



Figure <ref> shows the disk image for the  analysis of
the HLF reduction of GOODS South. The disk classification for each pixel
is shown on a black-to-blue colormap, with black corresponding to disk=0 and blue regions corresponding to disk=1.











  §.§.§ Irregular Image




Figure <ref> shows the disk image for the 
analysis of the HLF reduction of GOODS South. The irregular classification
for each pixel is shown on a black-to-green colormap, with black corresponding
to irregular=0 and green regions corresponding to irregular=1.






  §.§.§ Point source/compact Image



Figure <ref> shows the point source/compact image for
the  analysis of the HLF reduction of GOODS South. The point
source/compact classification for each pixel is shown on a black-to-yellow
colormap, with black corresponding to point source/compact=0 and
yellow regions corresponding to point source/compact=1.







  §.§.§ Morphological Color Image



Figure <ref> shows the morphological color image for the
 analysis of the HLF reduction of GOODS South. The false color image
is constructed following Section <ref>, with the pixel
intensities scaling with 1-background, the pixel hues set according to
the dominant class, and the saturation indicating the indeterminacy of the pixel
classification. Pixels with a single dominant class appear as bright red, blue,
green, or yellow for spheroid, disk, irregular, or point
source/compact classifications, respectively. Bright white pixels indicate
regions of the image where the model results were indeterminate in selecting a
dominant class. Dark regions represent pixels the model classified as background. We note that the pixel intensities are not scaled with the flux in
the image, and the per-object classifications require a local flux weighting
following Equation <ref> and the process described in
Section <ref>. This flux weighting usually results in a
distinctive class for each object, since the bright regions of objects often
have a dominant shared pixel classification. The outer regions of objects with
low flux show more substantial variation in the per-pixel classifications, but
these regions often do not contribute strongly to the flux-weighted per-object
classifications computed from this morphological color image.










§  PERFORMANCE



Given the data products generated by , we can perform a variety of
tests to quantify the performance of the model. There are basic performance
metrics relevant to how the model is optimized, reflecting the relative
agreement between the output of the model and the training data classifications.
However, given the semantic segmentation approach of  and the
pixel-level classification it provides, there are additional performance metrics
that can be constructed to mirror widely-used performance metrics in more
standard astronomical analyses including the completeness of sources detected by
 as regions of low background. In what follows, we attempt to address
both kinds of metrics and provide some ancillary quantifications to enable
translations between the performance of  as a deep learning framework
and as an astronomical analysis tool. In particular, we focus our analysis
on the 3D-HST catalog and HLF
reduction of the GOODS South region in the CANDELS Survey.




 §.§ Object-Level Morphological Classifications



The semantic segmentation approach of  provides classifications for
each pixel in an astronomical image. These pixel-level classifications can then
be combined into object-level classifications p using the flux-weighted
average described by Equation <ref>. The 
object-level classifications can then be compared directly with a test set of
visually-classified object morphologies provided by <cit.>.

To understand the performance of  relative to the K15 visual
classifications, we present some summary statistics of the training and test
sets pulled from the K15 samples. During training, the loss function used by
 is computed relative to the distribution of input K15
classifications for each object and not only their dominant classification. The
goal is to retain a measure of the uncertainty in visual classifications for
cases where the morphology of an object is not distinct.



  §.§.§ Distribution of Training Sample Classifications


Galaxies in the K15 training set have been visually classified by multiple
experts, providing a distribution of possible classifications for each object in
the sample. Figure <ref> presents histograms of the fraction of
K15 classifiers recording votes for spheroid, disk, irregular,
and point source/compact classes for each object. Only classes with more
than one vote are plotted.





  §.§.§ Classification Agreement in Training Sample


To aid these comparisons, we introduce the agreement statistic

    a(p) = 1 - H(p)/log(n_c)

where p is the distribution of classifications and n_c is the number of
classes. The quantity

    H(p)≡ -∑_k=1^n_c p_klog p_k

is the self entropy. According to these definitions, H(p)∈ [0,log
n_c] and a(p)∈[0,1]. The agreement a(p)→1 when the
distribution of classifications p is concentrated in a single class, and
a(p)→0 when the classifications are equally distributed. For reference,
a(p)≈0.57 for two equal classes and a(p)≈0.8 for a
90% / 10% split between two classes for n_c=5 possible classes.



  §.§.§ Training and Test Set Statistics


The K15 classifications have substantial variation in their agreement
a(p). Figure <ref> shows histograms and the
cumulative distribution of a(p) for objects with spheroid, disk, irregular, and point source/compact dominant classes. These
distributions of a(p) are roughly bimodal, consisting of a single peak
near a(p)=1 and a broader peak near a(p)≈0.5 with a tail
to larger a(p). As the cumulative distributions indicate, roughly
20%-60% of objects in the K15 sample had perfect agreement in their
morphological classification, with disk and point source/compact
being the most distinctive classes.



















The breadth in the agreement statistic for the input K15 data indicates
substantial variation in how expert astronomers would visually classify
individual objects. As these data are used to train , understanding
exactly what  should reproduce requires further analysis of the K15
data. An important characterization of the input K15 data is the confusion
matrix of object classifications. This matrix describes the typical
classification distribution for objects of a given dominant class. Figure
<ref> presents the confusion matrix for the K15
classifications, showing the typical spread in classifications for objects
assigned spheroid, disk, irregular, or point
source/compact dominant morphologies. For reference, a confusion matrix for a
distribution with perfect agreement is the identity matrix. Figure
<ref> provides some insight into the natural degeneracies
present in visually-classified morphologies. Objects with a dominant disk
classification are partially classified as spheroid (10%) and irregular (11%). The irregular objects frequently receive an alternative
disk classification (19%). The point source/compact objects also
are assigned spheroid classifications (14%). Objects with a dominant spheroid class have the highest variation and receive substantial
disk (18%) and point source/compact (11%) classifications. This
result is consistent with Figure <ref>, which shows a
relatively large disagreement for objects with a dominant spheroid
classification.


Since  is trained to reproduce the distribution of K15
classifications, the confusion matrix between the dominant 
classifications and the K15 classification distributions should be similar to
Figure <ref>. Indeed, Figure <ref> shows
the distribution of K15 classifications for objects with a given dominant
 classification agrees well with the input K15 distributions shown in
Figure <ref>. This result demonstrates 
reproduces well the intrinsic uncertainty in the K15 classifications, as
measured by the distribution of morphologies, recovered for a given K15
dominant classification.

The ability of  to reproduce the distribution of K15 classifications
is not the only metric of interest, as it does not indicate whether the
object-by-object  classifications agree with the K15 classifications
for objects with distinctive morphologies. Figure <ref>
shows that 20-60% of objects in the K15 classifications have an agreement
a(p)=1, meaning that all K15 visual classifiers agreed on the object
morphology. The confusion matrix for these distinctive objects constructed from
the K15 data is diagonal, and the confusion matrix for these objects constructed
from the  output should also be diagonal if  perfectly
reproduced the object-by-object K15 classifications. Further, to ensure
that  captures the distribution of the K15 morphologies, the
cumulative distribution of dominant K15 morphologies and dominant 
morphologies as a function of color were compared using a two-sample
Kolmogorov–Smirnov test. For each morphology, the p-values (p=0.3-0.99)
indicate consistency between the Morpheus and K15 distributions as a function of
color. These results suggest that  accurately captures the K15
representation of morphology without significant color bias.

To characterize the performance of  for the a(p)=1 K15
subsample, we used the  output classification images computed from
the HLF GOODS South images. The flux-weighted  morphological
classifications were computed following Equation
<ref> and using the K15 segmentation maps to
ensure the same pixels were being evaluated. Figure
<ref> presents the resulting confusion matrix showing
the  dominant classification for each object's dominant
classification determined by K15. As Figure <ref>
demonstrates,  achieves extremely high agreement (≥90%) with K15
for spheroid and point source/compact objects, and good agreement
(≥80%) for disk and irregular objects with some mixing
∼15% between them. This performance is comparable to other
object-by-object morphological classifications in the literature
<cit.>, but is constructed directly from a
flux-weighted average of pixel-by-pixel classifications by  using
real FITS image data of differing formats and depth.



  §.§.§ Redshift Evolution of Morphology in CANDELS Galaxies



To illustrate the scientific applications of ,
we examine the morphological distribution of ∼54,000
3D-HST sources in the five CANDELS fields as a function of redshift and stellar
mass (Figure <ref>). We combine together the
flux-weighted  classifications of galaxies identified in CANDELS with
the 3D-HST stellar masses and
redshift, dividing the sample into coarse redshift bins. The fraction of objects
N/N_tot with a flux-weighted classification of spheroid (red),
disk (blue), or irregular (green) are shown as a function of
stellar mass for each redshift bin, along with Poisson uncertainties on the
binned values. The well-established trends of increasing fractions of irregular
objects at small masses and high redshifts are correctly reproduced by
, as well
as the growth of the disk population at low redshifts.
These results can be compared
with the results
reported in Figure 3 of <cit.>.
To ensure comparable samples
between HC16 and this work, the -classified samples in
Figure <ref> are limited to
objects with H<24.5AB. Since HC16 and Morpheus use similar
but not
identical morphological classifications, we adapt the sample definitions
used by HC16 to the  classification scheme. To be counted as a
part of a morphological class, each galaxy's flux-weighted confidence value assigned by
 must be
greater than 0.7. This threshold ensures each classification is
mutually exclusive, but
low enough to ensure a comparable sample size to HC16.

The trends in Figure <ref> agree with those found by
HC16 in two important aspects. First, at lower redshifts, disks tend to dominate
spheroids, and as redshift increases, spheroids tend to dominate disks. Second,
irregular sources are a larger portion of the population than spheroids and
disks at lower stellar masses and more become less abundant than spheroids and
disks as stellar mass increases. The agreement between  and the
results of HC16, which were based on object-level classifications, confirms
the ability of  to capture source-level
morphologies by aggregating pixel-level classifications.






 §.§ Simulated Detection Tests



The  framework enables the detection of astronomical objects by
producing a background classification image, with source locations
corresponding to regions where background<1. If generating points in
the form of a source catalog are not supplied, the
segmentation algorithm of  uses an even more restrictive condition
that regions near sources must contain pixels with background=0. Given
that the semantic segmentation algorithm of  was trained on the K15
sample that has a completeness limit, whether the regions identified by
 to have background=0 correspond to an approximate flux
limit should be tested. Similarly, whether noise fluctuations lead to regions
assigned background≈0 in error should also be evaluated.

Below, we summarize detection tests for  using simulated images. For
these tests, a simulated sky background was generated using
Gaussian random noise with RMS
scatter measured in 0.5” apertures after convolving with a model HST PSF and
scaled to that measured from the K15 training images. The Tiny Tim software
<cit.> software was used to produce the PSF models appropriate for
each band.



  §.§.§ Simulated False Positive Test



Provided a large enough image of the sky, random sampling of the noise could
produce regions with local fluctuation some factor f above the
RMS background σ and lead to a false positive
detection. A classical extraction technique using aperture flux thresholds would
typically identify such regions as a SNR=f source. Here, we evaluate whether
 behaves similarly.

Using the Gaussian random noise field,
single-pixel fluctuations were added to the H-band only such that the
local flux measured in a 0.5” aperture after convolving with Tiny Tim
corresponded to SNR=[0.5, 1, 2, 3, 4, 5, 6, 7, 10]. The false signals were
inserted at well-separated locations such that  evaluated them
independently. The V, z, and J images were left as blank noise, and then
all four images were supplied to . We find that  assigns
none of these fake signals pixels with background=0. However, the
SNR=7 and SNR=10 regions have some background<1 pixels, and
while in the default algorithm,  would not assign these
regions segmentation maps, a more permissive version of the
algorithm could. An alternative test was performed by replacing the SNR=10
noise fluctuation in the H-band image with a Tiny Tim H-band PSF, added
after the convolution step with an amplitude corresponding to SNR=10 measured
in a 0.5” aperture. This test evaluates whether the shape of flux
distribution influences the detection of single-band noise fluctuations.
In this case, the minimum pixel values decreased to background≈0.05 for a single band SNR=10 fluctuation shaped like
an H-band PSF, but did not lead to a detection. We conclude that 
is robust to false positives arising from relatively large (SNR≲7)
noise fluctuations.



  §.§.§ False Negative Test



Given that  seems insensitive to false positives from noise
fluctuations, it may also miss real but low SNR sources. By performing a
similar test to that presented in Section <ref> but with
sources inserted in all bands rather than noise fluctuations inserted in a
single band, the typical SNR where  becomes incomplete for real
objects can be estimated.

Noise images were generated to have the same RMS noise as
the K15 images by convolving Gaussian random variates
with the Tiny Tim <cit.> model for the HST PSF. An array of
well-separated point sources modeled by the PSF were then inserted with a range
of SNR∈[1,25] into all four input band images. The  model was
then applied to the images, and the output background image analyzed to
find regions with background below some threshold value. Figure
<ref> shows the number of pixels below various background threshold values assigned to objects with different SNR. Below
about SNR∼15, the number of pixels identified as low background
begins to decline rapidly. We therefore expect  to show
incompleteness in real data for SNR≲15 sources. However, we emphasize
that this limitation likely depends on the training sample used. Indeed, the K15
training data set is complete to H=24.5AB in images with 5σ source
sensitivities of H≈27AB. If trained on deeper samples,  may
prove more complete to fainter magnitudes. We revisit this issue in Section
<ref> below, but will explore training  on deeper
training sets in future work.







 §.§ Morphological Classification vs. Surface Brightness Profile



In this paper, the  framework is trained on the K15 visual
classifications to provide pixel-level morphologies for galaxies. The K15
galaxies are real astronomical objects with a range of surface brightness
profiles for a given dominant morphology. Correspondingly, the typical
classification that  would assign to idealized objects
with a specified surface brightness profile is difficult to anticipate without
computing it directly. Understanding how  would classify idealized
galaxy models can provide some intuition about how the deep learning framework
operates and what image features are related to output 
classifications.

Figure <ref> shows the output  classification
distribution for simulated objects with circular <cit.> surface
brightness profiles, for objects with SNR=20, Sersic indices η∈[1,9],
and effective radii ranging from three to nine pixels. Synthetic FITS images for
each object in each band were constructed by assuming zero color gradients and a
flat f_ν spectrum, populating the image with a Sersic profile object and
noise consistent with the K15 images, and then convolving the images with a Tiny
Tim point spread function model appropriate for each input HST filter.

The results from  reflect common expectations for the typical Sersic
profile of morphological classes. Objects with η=1 were typically
classified as disk or spheroid, while intermediate Sersic index
objects (e.g., η≈2-3) were classified as spheroid. More
compact objects, with Sersic indices η≥4, were dominantly classified as
point source/compact. Also, as expected for azimuthally-symmetric
surface brightness profiles,  did not significantly classify any
objects as irregular. Figure <ref> provides a
complementary summary of the  classification of Sersic profile
objects, showing a matrix indicating the dominant classification assigned for
each pair of [η, R_e] values. The  model classifies large
objects with low η as disk, large objects with high η as spheroid, and small objects with high η as point source/compact.

Overall, this test indicates that for objects with circular Sersic profiles,
 reproduces the expected morphological classifications and that
asymmetries in the surface brightness are needed for  to return an
irregular morphological classification.












 §.§ Source Detection and Completeness



The semantic segmentation capability of  allows for the detection of
astronomical objects directly from the pixel classifications. In its simplest
form, this object detection corresponds to regions of the output 
classification images with low background class values. However,
the  object detection capability raises several
questions. The model was trained on the K15 sample, which has a reported
completeness of H=24.5AB, and given the pixel-by-pixel background
classifications computed by , it is unclear whether the
object-level detection of sources in images would match the K15 completeness. In
regions of low background, the transition to regions of high background likely depends on the individual pixel fluxes, but this transition
should be characterized.

In what follows below, we provide some quantification of the 
performance for identifying objects with different fluxes. To do this, we use
results from the 3D-HST catalog of sources for the GOODS South
<cit.>. Given the output 
background classification images computed from the HLF GOODS South FITS
images in F606W, F850LP, F125W, and F160W, we can report the
pixel-by-pixel background values and typical background values
aggregated for objects. These measurements can be compared directly with sources
in the <cit.> catalog to characterize how  detects
objects and the corresponding completeness relative to 3D-HST.






In a first test, we can locate the <cit.> catalog objects based
on their reported coordinates in the  background image, and
then record the background pixel values at those locations. Figure
<ref> shows the two-dimensional histogram of
 background value and 3D-HST source H-band AB magnitude,
along with the marginal distributions of both quantities. The figure also
indicates the reported K15 sample and 3D-HST 90% completeness flux levels. The
results demonstrate that for the majority of 3D-HST sources and for the vast
majority of bright 3D-HST sources with H<25, the local  background=0. The low background values computed by 
extend to extremely faint magnitudes (e.g., H≈29), indicating that for
some faint sources,  reports background=0 and that
background is not a simple function of the local SNR of an object. For
many objects with fluxes below the 3D-HST completeness, the  background value does increase with decreasing flux, and there is a
rapid transition between detected sources at H≈26.5 to undetected
sources at H≤27.5.

Owing to this transition in background with decreasing flux, the
completeness of  relative to 3D-HST will depend on a threshold in
background used to define a detection. Figure <ref>
shows the completeness of  in recovering 3D-HST objects as a function
of H-band source flux for different background levels
defining a  detection. The completeness flux limits for K15 and
3D-HST are indicated for reference. For magnitudes H<25AB, where 3D-HST and
K15 are complete,  is highly complete and recovers more than 99%
of all 3D-HST sources. The  completeness declines rapidly at fluxes
H>26.5AB, where  is 90% relative to 3D-HST for background
thresholds of P≤0.5. Perhaps remarkably, for all background
thresholds P≤0.01-0.5  detects some objects as faint as
H≈29, about 100× fainter in flux than the training set objects.

We further examined the detection of 3D-HST sources as a function of color
(V-H) to evaluate bias that may have been inherited as result of the training
dataset. In our tests, we found that Morpheus is not biased with respect to
color for those sources which are brighter than the K15 magnitude limit
(Figure <ref>). When considering all sources
within the 3D-HST catalog, Morpheus detects sources well, with a slight bias for
bluer sources, but performs less well for very red ((V-H) ≥ 9) and
((V-H) < 0) sources. However, it should be noted that there are very few
such sources in the training set, and with a more extensive training sample,
Morpheus could be more complete.







 §.§ Morphological Classification vs. Source Magnitude



The tests of  on simulated Sersic objects of different effective
radii and the completeness study suggest that the ability of  to
provide informative morphological information about astronomical sources will
depend on the size and signal-to-noise of the object. While these are intuitive
limitations on any morphological classification method, the distribution of
morphological classifications with source flux determined by  should
be quantified.

Figure <ref> shows the fraction of 3D-HST objects detected and
classified by  as spheroid, disk, irregular, and
point source/compact as a function of their H-band magnitude. Most of
the brightest objects in the image are nearby stars, classified as point
source/compact. At intermediate magnitudes,  classifies the objects
as primarily a mix of disk (∼50%) and spheroid
(∼30%), with contributions from irregular (∼10-30%) and
point source/compact (∼5-15%). For fainter objects, below the
completeness limit of the K15 training sample,  increasingly
classifies objects as irregular and point source/compact. This
behavior is in part physical, in that many low mass galaxies are irregular and
distant galaxies are physically compact. In part, it reflects a natural
bias in how the morphologies are defined during training. In K15, the class
point/source compact can describe bright stars and compact unresolved
sources (see Section 3.1 of K15). However, the trend also reflects how 
becomes less effective at distinguishing morphologies in small, faint objects
and returns either point source/compact and irregular for low SNR
and compact sources. While training  on fainter objects with
well-defined morphologies could enhance the ability of  to
distinguish the features of faint sources, the results of this test make sense
in the context of the completeness limit of the K15 training sample used.





 §.§ False Positives in GOODS South



The segmentation and deblending of real astronomical datasets
are challenging tasks.
An important test of the efficacy of the  segmentation and deblending algorithms
is to examine false positives generated when  is applied
to a real image. To quantify the
propensity for  to
generate false positives, the segmentation and deblending algorithms were run on
the HLF GOODS South image without the specified marker set parameter
p (See Algorithms <ref> and <ref>). For the purposes of this test, a
false positive is then defined as a set of pixels classified by the 
segmentation and deblending algorithms as a source but that
does not contain a source
from the 3DHST and CANDELS <cit.> catalogs.
Additionally, since the edges of the GOODS South
classified image are a frayed mix of pixels, to
minimize the effects of data artifacts sources
less than 20 pixels from the edge of the classified area were excluded from the
analysis. Further, we conservatively use the “default”
 algorithms that identify sources with
background=0, i.e., when  indicates a source
detection with high confidence.
With these choices, the sample used for the
false positive analysis was a total of 19,481 sources.

Among the objects classified by the segmentation and deblending
algorithms, 123 sources were not present in the CANDELS or 3D-HST catalogs. Upon
visual inspection of these sources, each can be categorized as an image
artifact, a poor deblend, a missed source, or an actual false positive. We list
the number of sources in each category in Table <ref>.

Sources in the image artifiact category are false positives
caused by image artifacts. The poor deblend category represents false
positives caused by the  deblending algorithm, where single sources
in the CANDELS or 3D-HST catalogs were shredded into multiple 
sources. The missed sources are  sources that upon visual
inspection correspond to real objects missed by the 3D-HST or CANDELS catalogs.
Sources in the actual false positive category are false positives not
associated with any image artifact or real source after visual inspection.

As Table <ref> shows, 
can identify real sources that other methods used
to generate catalogs can miss, although
the algorithms used by  can very rarely
cause actual false positives (at roughly the 0.1% rate).
Given the delicate nature of deblending, this analysis
suggests that the  deblending algorithm could be integrated with
other methods to generate more robust segmentation maps.









§ VALUE ADDED CATALOG FOR 3D-HST SOURCES WITH  MORPHOLOGIES



The  framework provides a system for performing the pixel-level
analysis of astronomical images and has been engineered to allow
for the processing of large-format scientific FITS data. As described in Section
<ref>,  was applied to the Hubble Legacy
Fields <cit.> reduction of HST imaging in GOODS
South[Some bright pixels in the released HLF images are censored with
zeros. For the purpose of computing the segmentation maps only, we replaced
these censored pixels with nearby flux values.] and a suite of morphological
classification images produced. Using the  background in GOODS
South, the detection efficiency of  relative to the
<cit.> 3D-HST catalog was computed (see Section
<ref>) and a high level of completeness demonstrated for
objects comparably bright to the <cit.> galaxy sample used to
train the model. By segmenting and deblending the HLF images,  can
then compute flux-weighted morphologies for all the 3D-HST sources.

Table <ref> provides the  morphological classifications for
50,506 sources from the 3D-HST catalog of <cit.>. This
value-added catalog lists the 3D-HST ID, the source right ascension and
declination, the F160W-band AB magnitude (or -1 for negative flux objects),
and properties for the sources computed by . The value-added
properties include a flag denoting whether and how  detected the
object, the area in pixels assigned to each source, and the spheroid, disk, irregular, point source/compact, and background
flux-weighted classifications determined by . The size of the
segmentation regions assigned to each 3D-HST object following Algorithms
<ref> and <ref> is reported for all objects. If
the segmentation region assigned to an object was smaller than a circle with a
0.36” radius, or the object was undetected, instead, use a 0.36” radius
aperture (about 109 pixels) to measure flux-weighted quantities. Only objects
with joint coverage in the HLF V, z, J, and H FITS images are classified
and receive an assigned pixel area. The full results for the 
morphological classifications of 3D-HST objects are released as a
machine-readable table accompanying this paper. Appendix
<ref> describes the  Data Release associated
with this paper, including FITS images of the classification images, the
value-added catalog, and segmentation maps generated by  for
the 3D-HST sources used to compute flux-weighted morphologies. Additionally, we
release an interactive online map at
https://morpheus-project.github.io/morpheus/https://morpheus-project.github.io/morpheus/, which provides an interface to examine the data and overlay the 3D-HST
catalog on the  classification images, morphological color images,
and segmentation maps.











§ DISCUSSION



The analysis of astronomical imagery necessarily involves pixel-level
information to be used to characterize sources. The semantic segmentation
approach of  delivers pixel-level separation between sources and the
background sky, and provides an automated classification of the source pixels.
In this paper, we trained  with the visual morphological
classifications from <cit.>. We then characterized the
performance of  in reproducing the object-level classifications of
K15 after aggregating the pixel information through flux-weighted averages of
pixels in -derived segmentation maps, and in detecting objects via
completeness measured relative to the 3D-HST catalog <cit.>. The
potential applications of  extend well beyond object-level
morphological classification. Below, we discuss some applications of the
pixel-level information to understanding the complexities of galaxy morphology
and future applications of the semantic segmentation approach of  in
areas besides morphological classification. We also comment on some features of
 specific to its application on astronomical images.



 §.§ Pixel-Level Morphology



The complex morphologies of astronomical objects have been described by both
visual classification schemes and quantitative morphological measures for many
years. Both <cit.> and <cit.> sought to
subdivide broad morphological classifications into more descriptive categories.
Quantitative morphological decompositions of galaxies <cit.>
also characterize the relative strength of bulge and disk components in
galaxies, and quantitative morphological classifications often measure the
degree of object asymmetry
<cit.>.

The object-level classifications computed by  provide a mixture of
the pixel-level morphologies from the  classification images. The
classification distributions reported in the  value-added catalog in
GOODS South provide many examples of flux-weighted measures of morphological
type. However, more information is available in the pixel-level classifications
than flux-weighted summaries provide.

Figure <ref> shows an example object for which the
 pixel-level classifications provide direct information about its
complex morphology. The figure shows a disk galaxy with a prominent central
bulge. The pixel-level classifications capture both the central bulge and the
extended disk, with the pixels in each structural component receiving dominant
bulge or disk classifications from . Note that  was not
trained to perform this automated bulge–disk decomposition, as in the training
process, all pixels in a given object are assigned the same distribution
of classifications as determined by the K15 visual classifiers. As the use
of pixel-level morphological classifications becomes wide-spread, the
development of standard datasets that include labels at the pixel-level will be
needed to evaluate the efficacy of classifiers. Simulations of galaxy formation
may be useful for generating such training datasets
<cit.>. We leave a more thorough analysis of
automated morphological decompositions with  to future work.










 §.§ Morphological Deblending



The ability of  to provide pixel-level morphological classifications
has applications beyond the bulk categorization of objects. One potential
additional application is the morphological deblending of overlapping objects,
where the pixel-level classifications are used to augment the deblending
process. Figure <ref> shows an example of two blended
objects, 3D-HST IDs 543 and 601, where the  pixel-level
classifications could be used to perform or augment star-galaxy separation. As
the figure makes clear, when  correctly assigns dominant
classifications to pixels, there exists an interface region between regions with
distinctive morphologies (in this case, spheroid and point
source/compact) that could serve as an interface between segmented regions in
the image. The deblending algorithm used in this work could include
other forms of machine learning
<cit.>
information in the deblending process. If  was trained on information
other than morphology, such as photometric redshift, those pixel-level
classifications could be used in the deblending process as well. We plan to
explore this idea in future applications of .




 §.§ Classifications Beyond Morphology



The semantic segmentation approach of  allows for complex features of
astronomical objects to be learned from the data, as long as those features can
be spatially localized by other means. In this paper, we used the segmentation
maps of K15 to separate source pixels from the sky, and then assigned
pixels within the segmentation maps the morphological classification determined
by K15 on an object-by-object basis. In principle, this approach can be extended
to identify regions of pixels that contain a wide variety of features. For
instance,  could be trained to identify image artifacts, spurious
cosmic rays, or other instrumental or data effects that lead to distinctive
pixel-level features in images. Of course, real features in images could also be
identified, such as the pixels containing arcs in gravitational lenses, or
perhaps low-surface brightness features in interacting systems and stellar
halos. These pixel-level applications of  complement machine
learning-based methods already deployed, such as those that discover and model
gravitational lenses
<cit.>. Pixel-level
photometric redshift estimates could also be adopted by  and compared
with existing methods based on SED fitting or other forms of machine learning
<cit.>.



 §.§ Deep Learning and Astronomical Imagery



An important difference in the approach of , where a purpose-built
framework was constructed from TensorFlow primitives, compared with the
adaptation and retraining of existing frameworks like Inception
<cit.> is the use of astronomical FITS images as
training, test, and input data rather than preprocessed PNG or JPG files. The
incorporation of deep learning into astronomical pipelines will benefit from the
consistency of the data format. The output data of  are also FITS
classification images, allowing pixel-by-pixel information to be easily
referenced between the astronomical science images and the  model
images. As indicated in Section <ref>, the  framework is
extensible and allows for any number of astronomical filter images to be used,
as opposed to a fixed red-blue-green set of layers in PNG or JPG files. The
 framework has been engineered to allow for the classification of
arbitrarily-sized astronomical images. The same approach also provides
 a measure of the dispersion of the classifications of individual
pixels, allowing the user to choose a metric for the “best” pixel-by-pixel
classification. The combination of these features allows for immense flexibility
in adapting the  framework to problems in astronomical image
classification.







§ SUMMARY AND CONCLUSIONS



In this paper, we presented , a deep learning framework for
the pixel-level analysis of astronomical images. The architecture of
 consists of our original implementation of a U-Net
<cit.> convolutional neural network.  applies the
semantic segmentation technique adopted from computer vision to enable
pixel-by-pixel classifications, and by separately identifying background and
source pixels  combines object detection and classification into a
single analysis.  represents a new approach to astronomical data
analysis, with wide applicability in enabling per-pixel classification of images
where suitable training datasets exist. Important results from this
paper include:


	
  *  provides pixel-level classifications of astronomical FITS
	images. By using user-supplied segmentation maps during training, the model
	learns to distinguish background pixels from pixels containing source
	flux. The pixels associated with astronomical objects are then classified
	according to the classification scheme of the training data set. The entire
	 source code has been publicly released, and a Python package
    installer for  provided. Further, we have a
    citable “frozen”
	version of code available through Zenodo <cit.>.
    
  * As a salient application, we trained  to provide
    pixel-level classifications of galaxy morphology by using the
    <cit.> visual morphological classifications of galaxies in
    the CANDELS dataset <cit.> as our training
    sample.
	
  * Applying  to the Hubble Legacy Fields
	<cit.> v2.0 reduction of the CANDELS data in GOODS South
	and the v1.0 data <cit.> for COSMOS, EGS,
	GOODS North and UDS, we generated morphological classifications for every
	pixel in the HLF mosaics.
	The resulting  morphological
	classification images have been publicly released.
	
  * The pixel-level morphological classifications in GOODS South were then
	used to compute and publicly release a “value-added” catalog of
	morphologies for all objects in the public 3D-HST source catalog
	<cit.>.
	
  * The CANDELS HLF and 3D-HST data were used to quantify the
	performance of , both for morphological classification and its
	completeness in object detection. As trained, the  code shows
	high completeness at magnitudes H≲26.5AB. We demonstrate that
	 can detect objects in astronomical images at flux levels up to
	100× fainter than the completeness limit of its training sample
	(H∼29AB).
	
  * Tutorials for using the  deep learning framework have been
	created and publicly released as Jupyter notebooks.
	
  * An interactive visualization of the  model results for
	GOODS South, including the  segmentation maps and pixel-level
	morphological classifications of 3D-HST sources, has been publicly released.


We expect that semantic segmentation will be increasingly used in astronomical
applications of deep learning, and  serves as an example framework
that leverages this technique to identify and classify objects in astronomical
images. We caution that  may be most effective at wavelengths
similar to the data
on which the model was trained (i.e., the F606W, F850LP, F125W, and F160W
bands).
However, <cit.> have shown recent success in applying
transfer learning on astronomical datasets with morphological labels. With the
advent of large imaging data sets such those provided by Dark Energy Survey
<cit.> and Hyper Suprime-Cam <cit.>, and
next-generation surveys to be conducted by Large Synoptic Survey Telescope
<cit.>, Euclid <cit.>,
and the Wide Field Infrared Survey Telescope <cit.>, pixel-level
analysis of massive imaging data sets with deep learning will find many
applications. While the details of the  neural network architecture
will likely change and possibly improve, we expect the approach of using
semantic segmentation to provide pixel-level analyses of astronomical images
with deep learning models will be broadly useful. The public release of the
 code, tutorials, and example data products should provide a basis
for future applications of deep learning for astronomical datasets.



B.E.R. acknowledges a Maureen and John Hendricks Visiting Professorship at the
Institute for Advanced Study, and NASA contract NNG16PJ25C.
R.D.H acknowledges the Eugene V. Cota-Robles Fellowship.
The authors acknowledge use of the lux supercomputer at UC Santa Cruz,
funded by NSF MRI grant AST 1828315. This work is based on observations taken by
the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS5-26555. This work is based on observations taken by the
3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS5-26555.


    Python <cit.>,
    NumPy <cit.>,
    astropy <cit.>,
    scikit-learn <cit.>,
    matplotlib <cit.>,
    Docker <cit.>,
    TensorFlow <cit.>










§ DEEP LEARNING



The  deep learning framework incorporates a variety of technologies
developed for machine learning applications. The following descriptions of deep
learning techniques complement the overview of  provided in Section
<ref>, and are useful for understanding optional
configurations of the model.





 §.§ Artificial Neuron



The basic unit of the  neural network is the artificial neuron (AN),
which transforms an input vector 𝐱 to a single output
AN(𝐱). The AN is designed to mimic the activation of a neuron,
producing a nonlinear response to an input stimulus value when it exceeds a
rough threshold.

The first stage of an AN consists of a function

    z(𝐱)  = ∑_i=1^n w_i x_i + b

that adds the dot product of the n-element vector 𝐱 with a vector
of weights 𝐰 to a bias b. The values of the 𝐰 elements
and b are parameters of the model that are set during optimization. The
function z(𝐱) is equivalent to a  linear
transformation on input data 𝐱.

In the second stage, a nonlinear function a is applied to the output of
z(𝐱). We write

    AN(𝐱) ≡ a(z(𝐱)),

where a(z) is called the activation function. The  framework allows
the user to specify the activation function, including the sigmoid

    sigmoid(z) = 1/1 + e^-z,

the hyperbolic tangent

    tanh(z) = e^z - e^-z/e^z + e^-z,

and the rectified linear unit

    relu(z) = max(0, z).

These functions share a thresholding behavior, such that the function activates
a nonlinear behavior at a characteristic value of z, but the domain of these
functions differ. For the morphological classification problem presented in this
paper, the rectified linear unit (Equation <ref>) was used as the
activation function.





 §.§ Neural Networks


Increasingly complex computational structures can be constructed from ANs.
Single ANs are combined into layers, which are collections of distinct ANs that
process the same input vector x. A collection of layers forms a
neural network (NN), with the layers ordered such that the outputs from one
layer provide the inputs to the neurons in the subsequent layer. Figure
<ref> shows a schematic of a NN and how the initial input vector
x is processed by multiple layers. As shown, these layers are
commonly called fully-connected since each neuron in a given layer receives the
outputs z from all neurons in the previous layer.








 §.§ Convolutional Neural Networks


The  framework operates on image data with a convolutional neural
network (CNN). A CNN includes at least one layer of ANs whose z
function uses a discrete cross-correlation (convolution) in place of the dot
product in Equation <ref>. For a convolutional artificial neuron (CAN), we
write

    z(𝐗) = (𝐗 * 𝐖) + b𝐉,

where 𝐗 * 𝐖 represents the convolution of an input image
𝐗 and a kernel 𝐖. The elements of the kernel 𝐖
are parameters of the model, and 𝐖 may differ in dimensions
from 𝐗. In , the dimensions of 𝐖 are set to be
3 × 3 throughout. The bias b is a scalar as before, and 𝐉
represents a matrix of 1s with the same dimensions as the result of the
convolution. In , the convolution is zero-padded to maintain the
dimensions of the input data.

The activation function of the neuron is computed element-wise after the
convolution and bias have been applied to the input. We write

    CAN(𝐗) ≡ a(z(𝐗)).

We refer to the output from a CAN as a feature map.

As with fully-connected layers, convolutional layers consist of a group of CANs
that process the same input data 𝐗. Convolutional layers can also be
arranged sequentially such that the output from one convolutional layer serves
as input to the next. ' neural network architecture,
being U-Net based, is comprised of CANs (see Figure <ref> for a
schematic). In typical convolutional neural network topologies, CANs are
used to extract features from input images. The resulting feature maps are
eventually flattened into a single vector and processed by a fully connected
layer to produce the output classification values.







 §.§ Other Functions in Neural Networks


The primary computational elements of  are a convolutional neural
network (Section <ref>) and a fully connected layer (Section
<ref>). In detail, other layers are used to reformat or summarize
the data, renormalize it, or combine data from different stages in the network.



  §.§.§ Pooling



Pooling layers (Figure <ref>) are composed of functions that
summarize their input data to reduce its size while preserving some information.
These layers perform a moving average (average pooling) or maximum (max pooling)
over a window of data elements, repeating these reductions as the window scans
through the input image with a stride equal to the window size. In the
morphological classification tasks described in this paper,  uses
2×2 windows and max pooling.





  §.§.§ Up-sampling


Up-sampling layers expand the size of feature maps by a specified factor through
an interpolation between input data elements. The up-sampling layers operate in
the image dimensions of the feature map and typically employ
bicubic and bilinear interpolation. In the morphological classification
application explored in this paper,  used 2×2 up-sampling and
bicubic interpolation.



  §.§.§ Concatenation



Concatenation layers combine multiple feature maps by appending them without
changing their contents. For instance, the concatenation of red, green, and blue
(RGB) channels into a three-color image would append three N× M images
into an RGB image with dimensions N× M × 3. This operation is used
in  to combine together data from the contraction phase with the
output from bicubic interpolations in the expansion phase (see Figure
<ref>).




  §.§.§ Batch Normalization


A common preprocessing step for neural network architectures is to normalize the
input data x using, e.g., the operation

    x̂ = (x-μ) / √(σ^2)

where x̂ is the normalized data, and μ and σ are parameters of
the model. <cit.> extended this normalization step to apply to the
inputs of layers within the network, such that activations (AN) and feature maps
(CAN) are normalized over each batch. A batch consists of a subset of the
training examples used during the training process. Simple normalization
operations like Equation <ref> can reduce the range of
values represented in the data provided to a layer, which can inhibit learning.
<cit.> addressed this issue by providing an alternative
normalization operation that introduces additional parameters to be learned
during training. The input data elements x_i are first rescaled as

    x̂_̂î = x_i - μ_𝐱/√(σ^2_𝐱+ϵ).

Here, x_i is a single element from the data output by a single AN or CAN over
a batch, μ_𝐱 is their mean, and σ^2_𝐱 is their
variance. The parameter ϵ is learned during optimization. The new
normalization BN_x̂_i is then taken to be a linear transformation

    BN_x̂_i = γ_xx̂_̂î + β_x.

The parameters γ_x and β_x are also learned during optimization.
<cit.> demonstrated that batch normalization, in the form of
Equation <ref>, can increase overall accuracy and decrease training
time, and we adopt this approach in the  framework.



 §.§ U-Net Architecture
 

The  framework uses a U-Net architecture, first introduced by
<cit.>. The U-Net architecture was originally designed for
segmentation of medical imagery, but has enjoyed success in other fields. The
U-Net takes as input a set of images and outputs a classification image
of pixel-level probability distributions. The architecture begins with a
contraction phase composed of a series of convolutional and pooling layers,
followed by an expansion phase composed of a series of
convolutional and up-sampling layers. Each of the outputs from the down-sampling
layers is concatenated with the output of an up-sampling layer when the height
and width dimensions of the feature maps match. These concatenations help
preserve the locality of learned features in the output of the NN.







§ CODE RELEASE



The code for  has been release via GitHub
(https://github.com/morpheus-project/morpheushttps://github.com/morpheus-project/morpheus).
 is also available as a python package installable via pip
(https://pypi.org/project/morpheus-astro/https://pypi.org/project/morpheus-astro/)
and as Docker images available via Docker Hub
(https://hub.docker.com/r/morpheusastro/morpheushttps://hub.docker.com/r/morpheusastro/morpheus).
 includes both a Python API and a command-line
interface, the documentation of which can be found online at
https://morpheus-astro.readthedocs.io/en/latest/https://morpheus-astro.readthedocs.io/en/latest/.



§ CODE TUTORIAL



An online tutorial demonstrating the  Python API in the form of a Jupyter
notebook can be found at

https://github.com/morpheus-project/morpheus/blob/master/examples/example_array.ipynbhttps://github.com/morpheus-project/morpheus/blob/master/examples/example_array.ipynb.
The tutorial walks through the classification of an example image. Additionally,
the tutorial explores other features of , including generating
segmentation maps and morphological catalogs.



§ DATA RELEASE



The data release associated with this work consists of multiple data
products. For each field in the CANDELS survey, we provide the following data
products: pixel-level morphological classifications, segmentation maps, and
value-added catalogs(see also Section <ref>) for the 3D-HST
catalogs. Tables <ref>-<ref> provide the URLs for each of
the data products. Each of the fields has two types of segmentation maps, a
segmentation map informed by the 3D-HST survey and a segmentation map informed
only by the background values provided by  (see Algorithm
<ref>). The classifications for the EGS and UDS fields may vary
as a result of using the F814W band in place of the F850LP due to availability.

An interactive online visualization of the HST images,  classification
images, and 3D-HST sources is available at
https://morpheus-project.github.io/morpheus/https://morpheus-project.github.io/morpheus/.






























