# An Automated Catalog of Long Period Variables using Infrared Lightcurves
from Palomar Gattini-IR

Aswin Suresh Department of Physics, Indian Institute of Technology Bombay,
Mumbai 400076, India  Viraj Karambelkar Cahill Center for Astrophysics,
California Institute of Technology, 1200 E. California Blvd. Pasadena, CA
91125, USA.  Mansi M. Kasliwal Cahill Center for Astrophysics, California
Institute of Technology, 1200 E. California Blvd. Pasadena, CA 91125, USA.
Michael C. B. Ashley School of Physics, University of New South Wales, Sydney
NSW 2052, Australia  Kishalay De NASA Einstein Fellow  MIT-Kavli Institute for
Astrophysics and Space Research, 77 Massachusetts Ave., Cambridge, MA 02139,
USA  Matthew J. Hankins Arkansas Tech University, Russellville, AR 72801, USA
Anna M. Moore  Research School of Astronomy and Astrophysics, Australian
National University, Canberra, ACT 2611, Australia  Jamie Soon Research School
of Astronomy and Astrophysics, Australian National University, Canberra, ACT
2611, Australia  Roberto Soria College of Astronomy and Space Sciences,
University of the Chinese Academy of Sciences, Beijing 100049, China  INAF-
Osservatorio Astrofisico di Torino, Strada Osservatorio 20, I-10025 Pino
Torinese, Italy  Sydney Institute for Astronomy, School of Physics A28, The
University of Sydney, Sydney, NSW 2006, Australia  Tony Travouillon Research
School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  Kayton K. Truong  Cahill Center for
Astrophysics, California Institute of Technology, 1200 E. California Blvd.
Pasadena, CA 91125, USA.

(Last up-dated ; in original form )

###### Abstract

Stars in the Asymptotic Giant Branch (AGB) phase, dominated by low to
intermediate-mass stars in the late stage of evolution, undergo periodic
pulsations, with periods of several hundred days, earning them the name Long
Period Variables (LPVs). These stars gradually shed their mass through stellar
winds and mass ejections, enveloping themselves in dust. Infrared (IR) surveys
can probe these dust-enshrouded phases and uncover populations of LPV stars in
the Milky Way. In this paper, we present a catalog of 159,696 Long Period
Variables using near-IR lightcurves from the Palomar Gattini - IR (PGIR)
survey. PGIR has been surveying the entire accessible northern sky
($\delta>-28^{\circ}$) in the _J_ -band at a cadence of 2–3 days since
September 2018, and has produced J-band lightcurves for more than 60 million
sources. We used a gradient-boosted decision tree classifier trained on a
comprehensive feature set extracted from PGIR lightcurves to search for LPVs
in this dataset. We developed a parallelized and optimized code to extract
features at a rate of $\scriptstyle\mathtt{\sim}$0.1 seconds per lightcurve.
Our model can successfully distinguish LPVs from other stars with a true
positive rate and weighted g-mean of 0.95. 73,346
($\scriptstyle\mathtt{\sim}$46%) of the sources in our catalog are new,
previously unknown LPVs.

##  1 Introduction

Stars in the final stages of stellar evolution exhibit periodic variations in
their brightness on timescales longer than 100 days, and are hence termed as
Long Period Variables (LPVs, Samus’ et al. 2017). LPVs are generally
categorized into three types: Mira variables, semiregular variables, and
OGLE111Optical Gravitational Lensing Experiment (Udalski et al., 1992) small
amplitude red giants (Soszyński et al., 2009; Wray et al., 2004). Miras can
show oscillation amplitudes larger than ten magnitudes in the optical bands,
partly due to absorption of the optical light by nearby molecules. In the
infrared, however, the brightness variation is less pronounced. Semiregular
variables, which typically display variable amplitudes and periods, and small
amplitude red giants are composed primarily of Asymptotic Giant Branch (AGB)
stars or Red Giant Branch (RGB) stars and typically show oscillations with
amplitudes less than 0.2 mag (Smak, 1966).

The advent of time domain surveys has resulted in a number of variable star
catalogs covering a wide range of optical variability in stellar populations.
The Massive Compact Halo Object (MACHO) survey (Alcock et al., 2001) found
20000 LPVs in the Large Magellanic Cloud (LMC). Wood et al. (1999) used data
from 1500 LPVs in the LMC to establish several linear sequences in the
log(period) - log(luminosity) space of LPVs. OGLE (Udalski et al., 1992)
identified 91995 LPVs in the Large Magellanic Cloud (LMC) from the OGLE-III
catalog of variable stars by examining the log(period) - log(Wesenheit index)
relation, finding 1667 Miras in the LMC (Soszyński et al., 2009). The Vista
Variables in the Vía Láctea (VVV) survey (Dekany et al., 2011) discovered 129
Miras in the Galactic bulge, along with producing a catalog of 1013 LPVs in
the Milky Way using point-spread function photometry data in the
$\text{K}_{\text{s}}$ band from 10 tiles around the Galactic bulge and studied
the period-amplitude relationship of Miras in the near-infrared (Nikzat et
al., 2022). The largest catalog of LPVs to date, consisting of over 1.7
million LPV candidates, of which 392,240 have published periods, was presented
in Lebzelter et al. (2023) using 34 months of data from the third Gaia data
release. Karambelkar et al. (2019) identified 417 extremely luminous long
period variables in nearby galaxies from the SPitzer InfraRed Intensive
Transients Survey (Kasliwal et al., 2018) with the _Spitzer Space Telescope_
and extended the period-luminosity relationship of AGB stars to longer periods
and luminosities. Several catalogs of variable stars consisting of a mixture
of short-period (<100 days) stars, such as Cepheids and RR Lyrae, as well as
longer-period variables, have been identified (see for e.g. Chen et al., 2018,
Chen et al., 2020, Pojmanski et al., 2005). The wide parameter space of
stellar variability and increasingly large datasets from time-domain surveys
present the biggest challenge in these studies to identify and characterize
variable stars.

The efficacy of Machine Learning (ML) techniques in identifying and
classifying variable stars has been demonstrated in multiple works. For
instance, Masci et al. (2014) used a random forest classifier trained on
lightcurve morphology-based features and Fourier decomposition to distinguish
between Algols, RR Lyrae, and W Ursae Majoris among 8273 variables from WISE.
Debosscher et al. (2007) employ a multivariate Gaussian mixture classifier to
classify between 35 classes of variable stars solely using Fourier amplitudes
and phases as the feature set. Sánchez-Sáez et al. (2021) classified 868,371
sources using lightcurves from the Zwicky Transient Facility (Bellm et al.,
2019) using a multi-level classifier architecture trained on 152 features.
Boone (2019) developed a gradient-boosted decision tree classifier trained on
features extracted from Gaussian process augmented multi-band lightcurves to
classify between several variable and transient classes ranging from Miras and
RR Lyrae to TDEs, SNe Ia, and kilonovae.

In this paper, we employ an ML approach to search for Long Period Variables
using data from the Palomar Gattini-IR (PGIR) survey taken between September
2018 and July 2021. PGIR’s survey of the northern sky in the Near-Infrared J
band, at a cadence of 2 days to a depth of 16 AB mag, has produced lightcurves
of more than 60 million stars. We utilize a machine learning framework to
identify LPVs in this lightcurve archive. We build a comprehensive training
set consisting of LPVs, other erratic and non-periodic variables such as RCrB
stars, Young Stellar Objects and non-variables. We train a gradient-boosted
decision tree to distinguish between these classes, and use its predictions to
generate a catalog of 159,696 LPVs. Of these, 73,346 (45.9% of the full
catalog) are newly identified LPVs. We describe the training sample and
feature set for the machine learning classifier in Section 2. We explain the
training procedure, iterative build-up of the training set, and the relative
importance of various features in Section 3. We analyze the features of LPVs
in the PGIR catalog and assess the validity, completeness, and purity of the
catalog by comparison with other catalogs in Section 4. Section 5 summarizes
our conclusions.

##  2 Training Sample

The default data processing pipeline of PGIR (De et al., 2020) produces light
curves for sources detected in multi-epoch images as part of ‘match files’.
Briefly, the reference images of each field is used to construct a master
source catalog, and spatially cross-matched to the source catalog for every
observation of the field as part of survey operations. The resulting cross-
match produces light curves for every source detected in the reference image
covering the epochs where the source is detected in the individual epochs. In
this paper, we use data from the beginning of PGIR survey operations until 14
July 2021, giving a baseline of $\scriptstyle\mathtt{\sim}$1400 days for each
lightcurve, which is well-suited for characterizing the few hundred day
periods of LPVs.

Constructing a catalog of Long Period Variables requires distinguishing LPVs
from non-periodic variable stars (referred to as “non-LPVs" hereafter) and
non-variable stars (“non-vars" hereafter) with high accuracy. The vast
parameter space spanned by these stellar lightcurves makes machine learning a
viable option to learn intricate features distinguishing the three classes.
The starting point for building the machine learning classifier is assembling
a training dataset representative of the various characteristics of LPVs, non-
LPVs and non-vars. We build this training set iteratively. We first train a
classifier on a small set of bonafide members of each category. We use this
classifier to predict the classes of a larger sample of sources, visually
examine the results and use the correctly classified sources to expand the
initial training set. The initial training set is described in Section 2.1 and
the final, expanded dataset is described in 3.2.

###  2.1 Initial bonafide training set

We begin with a relatively small training set of 1344 sources, consisting only
of LPVs and non-LPVs. These were assembled as part of ongoing searches for the
erraticaly varying R Coronae Borealis (R Cor Bor, Clayton 2012) stars and
other erratic, large-amplitude Galactic variables using PGIR (Karambelkar et
al., 2021, Earley et al. 2024, in prep.). Our “LPV" set includes stars
identified as LPVs in these studies, primarily based on visual inspection of
their lightcurves and Lomb Scargle (Lomb, 1976; Scargle, 1982) periodograms.
We also include any additional late M-type stars that were spectroscopically
classified as part of these searches. The “non-LPV" set consists of all known,
spectroscopically confirmed R Cor Bor stars that show erratic variations in
PGIR data and a smattering of other erratically varying stars such as young
stellar objects. LPVs constitute the majority sample of this initial training
set, with 1265 samples, and non-LPVs form the minority class with 79 samples
(5.87% of the full dataset). Representative examples of the two classes are
shown in Figure 1.

Figure 1: Examples of LPVs and non-LPVs from the bonafide set. LPVs display an
approximately sinusoidal morphology in their lightcurves, as can be seen in
the figure on the left, while non-LPVs show erratic behavior, as the
spectroscopically confirmed RCorBor variable shown on the right. The blue line
and shaded blue region represent the mean of the Gaussian process regression
fit and the 95% confidence interval respectively.

This training sample may not be sufficient to cover the breadth of parameters
spanned by the two classes. Hence it is used as a baseline dataset to obtain
broad classifications. The classifier will inevitably have a high
misclassification rate in the first few training iterations before a more
balanced dataset enables it to learn the characteristics of both classes. To
mitigate this, we built a bigger dataset iteratively by visually examining the
predictions of previous classifiers and re-labeling misclassified sources
manually (see Section 3.2).

###  2.2 Feature Extraction

We use a set of 16 features for initial training iterations. We begin by
fitting a Gaussian Process (GP) (Williams & Rasmussen, 1995) to model the
lightcurves. Gaussian processes have been shown to be an effective method for
astronomical lightcurve modeling (see for e.g. Boone, 2019, Qu et al., 2021,
Foreman-Mackey et al., 2017). We use the python package george (Ambikasaran et
al., 2015) to implement the GP fits. The kernel of choice is a Matern 3/2
kernel with a length scale set to the period extracted from the raw
lightcurve. A uniform time sampling of 500 points from the first observation
to the most recent observation is used to calculate the predicted J-band
magnitudes. The minimum and maximum slopes and peak-to-peak amplitudes are
calculated from the augmented lightcurve. The goodness of fit of the Gaussian
process is quantified by the $\textit{R}^{2}$ score and is saved as the “GP
score". We calculate the difference between the $90^{th}$ percentile and
$50^{th}$ percentile of the J-band magnitudes of each lightcurve as the
percentile difference.

Periodicity-based features are extracted using the Lomb-Scargle periodogram,
implemented using the gatspy library (VanderPlas & Ivezić, 2015). The
augmented lightcurve, with J band magnitude, converted to flux, is fit to
obtain the three most prominent periods and their associated Lomb-Scargle
score. A sinusoid with the best fit period is fit to the lightcurve to
calculate the reduced $\chi^{2}$ of the fit. Additionally, the reduced null
$\chi^{2}$ is calculated as

$\chi^{2}_{red,null}=\sum_{i}\left(\frac{y_{i}-\hat{y}}{\sigma_{i}}\right)^{2}$
(1)

where ${y_{i}}$ and $\sigma_{i}$ are the fluxes and errors respectively, and
$\hat{y}$ is the error weighted mean of J band flux.

The final periodicity-based features used from these calculations are the
Lomb-Scargle score of the best period - LS score, reduced $\chi^{2}$ and
reduced null $\chi^{2}$, ratio of the best fit period to the maximum period
(defined as the observation baseline) - period ratio, ratio of Lomb-Scargle
score of the two most prominent periods - Lomb-Scargle score ratio and ratio
of reduced $\chi^{2}$ and reduced null $\chi^{2}$ \- ${\chi^{2}}$ ratio.
Phase-folded lightcurves of LPVs show a clean sinusoidal evolution. The
$\chi^{2}$ value of a sinusoidal fit to the phase folded lightcurve is saved
as “phase $\chi^{2}$" and is useful to distinguish non-LPVs, non-variables,
and erratic variables, which have high phase $\chi^{2}$ values compared to
LPVs.

Another useful measure of variability is given by the Stetson L index
(Stetson, 1996):

$L=\frac{JK}{0.789}$ (2)

where,

$J=\sum_{n=1}^{N-1}{\text{sign}(\delta_{n}\delta_{n+1})\sqrt{|\delta_{n}\delta_{n+1}|}}$
(3)

and

$K=\frac{1/N\sum_{i=1}^{N}{|\delta_{i}|}}{\sqrt{{1/N}\sum_{i=1}^{N}{{\delta_{i}}^{2}}}}$
(4)

are the Stetson J and K indices (Stetson, 1996), with $\delta_{i}$ being the
difference between the magnitude of the $i^{th}$ point and the mean magnitude
of the lightcurve, scaled by the error of the $i^{th}$ point.

Figure 2: Distributions of select features for variables in the bonafide set.
Features such as the Reduced $\chi^{2}$ and the peak-to-peak amplitude can be
seen to separate the LPVs and non-LPVs while the period ratio and the von
Neumann score occupy similar values for both classes.

The final feature extracted from lightcurves is the von Neumann score (von
Neumann, 1941) calculated as

$\eta=\frac{\sum_{n=1}^{N-1}(x_{n+1}-x_{n})^{2}/(N-1)}{\sigma^{2}}$ (5)

. $\eta$ is a measure of the correlation between successive observations, and
positive correlations lead to small values of $\eta$.

Using a fully parallelized and optimized feature extraction process, our code
can calculate all features for each lightcurve in $\sim 0.1$ seconds. The
distribution of select features for LPVs and non-LPVs in the bonafide set is
shown in Figure 2

We supplement the features obtained directly from PGIR lightcurves using color
information of variable stars observed by other ground and space-based
missions. In particular, we query the J, H, and K band colors of variables
from the 2-Micron All Sky Survey (2MASS) (Skrutskie et al., 2006). The 2MASS
All-Sky Catalog of Point Sources (Cutri et al., 2003) contains over 300
million objects, including variable stars. We query the infrared photometry of
crossmatched sources and calculate the J-H and J-K colors. We note that the
J-H and J-K colors are highly correlated in the bonafide set. The Wide-field
Infrared Survey Explorer (WISE) (Wright et al., 2010) is a space-based
infrared observatory surveying the infrared sky at 3.4 (W1), 4.6 (W2), 12 (W3)
and 22 (W4) ${\mu}m$ wavelengths. We query the WISE All-sky data release
(Cutri & et al., 2012) to obtain W1-W2, W1-W4, and W3-W4 colors. For cases
where no crossmatch is present in WISE or 2MASS, the median color value of the
entire dataset is used.

Table 1: Summary of features used to train the classifier.  Features | Description
---|---
von_neumann_score | von Neumann score $\eta$
ptp | Peak-to-peak amplitude
slope_min | Minimum slope of lightcurve
slope_max | Maximum slope of lightcurve
gp_score |  $\textit{R}^{2}$ of GP fit
per_diff_90_50 |  Difference between $90^{th}$ and $50^{th}$ percentiles of J-band magnitude
LSscore | Lomb-Scargle score of best period
redchi2 | Reduced $\chi^{2}$ of sinusoid fit
rednullchi2 |  Reduced $\chi^{2}$ assuming constant expected value
chi2ratio | Ratio of redchi2 and rednullchi2
period_ratio |  Ratio of best period and maximum period
LSratio |  Ratio of Lomb-Scargle score of best two periods
phase_chisq |  $\chi^{2}$ of sinusoid fit to phase folded lightcurve
L_m | Stetson L index
J-H | 2MASS J-H color
J-K | 2MASS J-K color
W1-W2 | WISE W1-W2 color
W1-W4 | WISE W1-W4 color
W3-W4 | WISE W3-W4 color

###  2.3 Data Resampling

To overcome the class imbalance in the initial training set, we use synthetic
sampling methods implemented using the imbalanced-learn package (Lemaître et
al., 2017). Multiple strategies exist to selectively upsample the minority
class and undersample the majority class. Based on trial and error, we choose
the Adaptive Synthetic (ADASYN) sampling method (He et al., 2008) for
upsampling and all K-Nearest Neighbours (allKNN) method (Tomek, 1976) for
undersampling.

ADASYN sampling creates new synthetic samples for the minority class by
interpolating between existing minority samples. This interpolation is
performed by considering the neighboring minority samples and creating
synthetic samples that lie along the line connecting them. ADASYN is unique in
the sense that it adjusts the number and distribution of synthetic samples
based on the difficulty of classification. It focuses on generating synthetic
instances in the regions of the feature space that are challenging for the
classifier. By doing so, it helps to balance the class distribution while
maintaining the overall integrity of the data, resulting in improved model
performance on imbalanced datasets.

The allKNN downsampling algorithm works by assigning a minority density to
each sample in the minority class. The minority density is calculated as the
number of minority samples in k-nearest neighbors, where “k" is a user-
controlled parameter. Samples from the majority class are selectively removed
if they fall into a minority neighborhood with a density below a set
threshold. Samples with density below the threshold could be sparsely located
or overlap with the majority class, making it hard for the model to learn.
This process is iterated till a desirable balance is achieved. allKNN
downsampling, therefore helps create a balanced dataset while also removing
samples that are harder to learn.

This combination is used to obtain a more balanced training set used to train
the classifier.

##  3 Machine Learning Framework

Decision trees and random forest-based classifiers have been shown to be
effective for the task of classifying variable star lightcurves (see for e.g.
Sánchez-Sáez et al., 2021, Boone, 2019, Masci et al., 2014, Debosscher et al.,
2007). Decision trees (Breiman et al., 1984) are intuitive models having a
tree-like branched node structure. They recursively partition the input data
at each node to create homogeneous subgroups based on criteria set from the
features. Data can be split either to maximize information gain or minimize
impurity in partitioned sets. Complex decision trees may be prone to
overfitting. Parameters such as the maximum depth of the tree and the minimum
number of samples required to create a split can be controlled to avoid
overfitting data. Random forests (Breiman, 2001) is an ensemble method that
builds multiple decision trees and combines their predictions to assign a
probability to the predicted output. Different trees in the forest are built
from random subsets of the training set and features. This helps avoid
overfitting and gives improved performance. The importance of features in
separating various classes in the input data can also be quantified in random
forests.

We employ a gradient-boosted decision tree classifier implemented using the
LightGBM framework (Ke et al., 2017). Gradient boosting is a popular technique
for creating ensemble decision trees. It sequentially builds weak-learner
trees where each tree is trained to correct the residuals of the preceding
tree, “boosting" its performance. Gradient descent is used to minimize a loss
function, which is how the model learns. LightGBM offers an efficient
implementation of gradient-boosted decision trees and reduces memory
consumption by binning features.

For the training process, we split the bonafide dataset in a 70:30 train-test
split. It is important that the model is not tested on resampled data, hence
resampling is done only on 70% of the bonafide set. LPVs are labelled as 1 and
non-LPVs as 0. Due to class imbalance, the overall accuracy of the classifier
may be misleading. We use the following metrics to quantify model performance:

$\textrm{Recall / True positive rate
(TPR)}=\frac{\textrm{TP}}{\textrm{TP}+\textrm{FN}}$ (6)

$\textrm{True negative rate
(TNR)}=\frac{\textrm{TN}}{\textrm{TN}+\textrm{FP}}$ (7)

$\textrm{Precision (prec)}=\frac{\textrm{TP}}{\textrm{TP}+\textrm{FP}}$ (8)

$\textrm{F-measure}=\frac{2*\textrm{prec}*\textrm{recall}}{\textrm{prec}+\textrm{recall}}$
(9)

$\textrm{Weighted
g-mean}=(\textrm{TPR}^{\textrm{P}}*\textrm{TNR}^{\textrm{N}})^{\frac{1}{\textrm{P}+\textrm{N}}}$
(10)

where TP is true positive, FP is false positive, TN is true negative, FN is
false negative, P is all positives, and N is all negatives.
The true positive rate is the probability of an LPV being labeled correctly,
and the true negative rate is the probability of a non-LPV being labeled
correctly. Precision gives the fraction of all LPVs classified correctly. The
F-measure is an accuracy metric given by the harmonic mean of precision and
recall. The g-mean is the geometric mean of TPR and TNR. We observe that due
to class imbalance in the test set, the simple geometric-mean tends to favor
models that have a lower misclassification rate of non-LPVs. To mitigate this
bias, we calculate the weighted geometric mean of TPR and TNR, using the total
samples in each class as weights to get a robust measure of classifier
performance.

###  3.1 Training on the bonafide set

The LightGBM gradient-boosted decision tree has a large number of
hyperparameters that need to be tuned to get the most out of the classifier. A
grid search that chooses the optimal hyperparameters by trying each possible
combination would be time-consuming and is not feasible. In each training
iteration, we randomly sample 200 points from the hyperparameter space and
choose the set with the highest accuracy. This method has been shown to be
effective at achieving near-optimal performance with high probability
(Bergstra & Bengio, 2012). For each combination of hyperparameters, accuracy
is calculated using five-fold cross-validation. The training dataset is
partitioned into five folds, with the model being fit to four such folds and
evaluated on the remaining fold. This process is repeated to effectively
sample the training data. We evaluate the metrics on the model and create a
confusion matrix. For binary classification tasks, the confusion matrix is a 2
x 2 matrix showing the distribution of true negatives, false positives, false
negatives, and true positives. Examining the confusion matrix gives good
insights into the model’s performance, as all metrics are derived from its
quantities. Figure 3 shows the confusion matrix from training on the bonafide
dataset.

Figure 3: Confusion matrix from training on the bonafide dataset. This model
will be referred to as LGBMv0.1. Metrics: TPR = 97.120%, TNR = 95.455%, prec =
99.731%, F-measure = 0.984, weighted g-mean = 0.970  Figure 4: Examples of
LPVs classified as Type II LPVs. Type II LPVs display sinusoidal evolution,
however, additional trends such as a constant rise or deacy, varying
amplitudes of oscillation or presence of several oscillation frequencies are
also present.

Following this training procedure, it is important to understand samples
misclassified by the model. We save the predicted probabilities and labels for
all samples in the training set and compare labels with the true values.
Lightcurves of all misclassified samples are plotted and saved along with the
confusion matrix. Since we perform a randomized search for hyperparameters, we
run multiple training iterations and pick the best-performing model.
Examination of confusion matrix and misclassified samples is particularly
helpful in choosing the optimal model in cases where certain samples may be
mislabelled in the training set. This is particularly unavoidable when we
build the training set in an iterative process and predictions from initial
models form the true labels of the next iteration.

###  3.2 Expanding the training set

The first step in expanding the training dataset is to choose a larger set of
data to obtain predictions on. We extracted features for 35 million
lightcurves in the PGIR lightcurve database. Of these, we choose a subsample
of $\sim$12000 sources by applying cuts on the von Neumann score ($\eta$ <
1.5), peak-to-peak amplitude (ptp > 1) and Lomb-Scargle score (LSscore >
0.75). We run the LGBMv0.1 classifier on all of these sources and choose only
those samples with a prediction probability >0.9, resulting in a larger
training set of 9815 LPVs and 443 non-LPVs.

However, upon visual inspection of the classifications, we find that the
misclassification rate is high in this dataset. We observe that a large
fraction of the classified non-LPVs are LPVs that show variations in
amplitude, show multiple pulsation modes or have a slow background decay or
rise in addition to the few hundred day periodic variations. Upon training on
the extended dataset with corrected labels, we observe that the classifier
struggles to separate these erratic LPVs, generally deviating from a clean
sinusoidal evolution, from non-LPVs. Hence we create a third category – “Type
II LPVs", based purely on lightcurve morphology comprising of such non-
sinusoidal LPVs. We note that this classification does not necessarily hold a
physical meaning, as many AGB stars typically do display such evolution and
can have multiple simultaneous periodic excitations. However, separating the
LPVs into cleanly sinusoidal and erratic evolution better enables the
classifier to pick up the general features of each class. We combine the
sinusoidal LPVs and type II LPVs into the LPV class following classification.
Representative examples of type II LPVs are shown in Figure 4.

We repeat the training procedure on this larger sample of labeled data. In
this iteration, we observe that the classifier has improved performance and is
better able to distinguish the type II LPVs and non-LPVs, as shown in the
confusion matrix in Figure 6. We assemble a clean set of non-vars by visually
inspecting a subset of lightcurves chosen to have von Neumann score less than
1.5, amplitude less than 1 mag and LPV probability less than 0.3 from the
classifier trained on bonafide variables. Finally, we build a master training
set consisting of 2335 LPVs, 444 Type-II LPVs, 1332 non-vars and 166 non-LPVs.
The distribution of select features for LPVs and non-LPVs in the bonafide set
is shown in Figure 5. Several training iterations of the model on this
training set resulted in a best-performing model with a g-mean score of 0.95.

Figure 5: Distributions of select features for variables in the expanded
training set. We observe that the features occupy similar values for the
classes shown here as the classifier has to separate low amplitude variables
whose feature phase space is not too different from non-variables.

The model shows high confusion between LPVs and Type II LPVs, but is able to
minimize confusion between the combined set of these two and the combined set
of non-vars and non-LPVs. We sum the probability score of LPV and type II LPV
is to calculate the net probability of a source being an LPV. This allows us
to separate LPVs from non-LPVs and non-vars in an efficient manner. We group
all of the type II LPVs into the LPV class and the non-vars into the non-LPV
class to obtain the final confusion matrix and metrics shown in Figure 7. We
also list the optimal parameters of the gradient-boosted classifier in Table
2.

Figure 6: Confusion matrix from training on the extended labeled set. This
model will be referred to as LGBMv0.2. Since we classify between three
classes, we calculate metrics for each class. Metrics (listed for LPVs, Type
II LPVs, non-vars and non-LPVs respectively): TPR = 70.80%, 65.81%, 93.86%,
95.35%, TNR = 94.71%, 87.40%, 91.72%, 98.02%, precision = 94.12% 35.0%,
84.51%, 63.08%, F-measure = 0.81, 0.46, 0.89, 0.76 and weighted g-mean = 0.80,
0.69, 0.93, 0.95  Figure 7: Confusion matrix from training on the extended
labeled set with Type II LPVs grouped with LPVs and non-LPVs grouped with non-
vars. Metrics: TPR = 94.89%, TNR = 88.78%, precision = 82.56%, F-measure =
0.88, weighted g-mean = 0.95

###  3.3 Feature Importance

Figure 8: Feature importance of various features used in training. The feature
importance quantifies how effective a particular feature is in separating an
LPV and a non-LPV. We notice features such as the von Neumann score and period
ratio are very effective, while the Lomb-Scargle score is not as useful for
the classifier.

Tree-based classifiers provide convenient methods to quantify the relative
importance of different features used in training. The LightGBM framework, in
particular, calculates feature importance in two ways, using “split" and
“gain". Split-wise importance is calculated as the number of times a
particular feature is used in the model to partition data. Since the model
tries to split data to minimize a loss function, the number of times a feature
is used in this process is a good measure of its effectiveness. Gain-wise
importance calculates the impurity reduction in each split; impurity reduction
is the difference in loss of parent tree node and the weighted sum of loss of
child tree nodes after the split. This quantifies how much a split improves
the model’s performance. Impurity reduction for each feature is accumulated
for all trees in the ensemble and normalized to obtain the final relative
importance. We adopt gain-based importance as the preferred metric as it
provides a direct measure of improvement in accuracy contributed by a feature.

Figure 8 shows the calculated feature importance. We note that the period
ratio is among the most important features for the classifier through both the
gain and split importance calculations. This is consistent with other variable
star or transient alert classifiers (e.g. Sánchez-Sáez et al., 2021, Masci et
al., 2017, Kim, Dae-Won & Bailer-Jones, Coryn A. L., 2016). The $\chi^{2}$ of
a sinusoid fit to the phase folded lightcurve and the J-H color also rank
highly in both metrics, while the Stetson L index and peak-to-peak amplitude
perform decently. We also observe that the Lomb-Scargle scores and percentile
difference perform poorly in separating LPVs and non-LPVs.

Table 2: Optimal hyperparameters of the LightGBM classifier.  Hyperparameter | Value
---|---
objective | multiclass
boosting_type | gbdt
colsample_bytree | 0.9
max_depth | 9
min_child_samples | 40
min_split_gain | 1
n_estimators | 460
num_leaves | 197
reg_alpha | 0.035
reg_lambda | 0.03
subsample | 0.90
subsample_freq | 59

##  4 The Catalog

Figure 9: Spatial distribution of long-period variables from Gattini. A large
fraction (152,220) of LPVs from the catalog lie within the galactic latitudes
of $-10^{\circ}$ to $+10^{\circ}$. The missing patch of sky is below a
declination of $-28^{\circ}$, which PGIR cannot observe. The color bar
represents the density of stars, with a high density around the galactic
centre.

We use the LGBMv0.2 classifier to predict the classes of 35 million PGIR
lightcurves, chosen to have a von Neumann score less than 2, and requiring at
least 30 detections. We present a catalog of 159,696 stars having
$\textit{prob\\_lpv}>0.5$ from the PGIR survey in Table 3, containing the
following columns. “Name" refers to the internal designation of the variable,
“RA" and “Dec" give the spatial position, “num_detections" is the number of
detections, “lcdur" is the baseline of observation, given by the difference of
the time of final detection and the time of initial detection, “bestperiod"
and “amplitude" are the period and amplitude, “$\left<J\right>$" is the mean J
band magnitude, “prob_lpv" is the classifier confidence. In addition, all
extracted features and 2MASS and WISE color information are also presented.
Figure 10 shows example lightcurves of 20 LPVs from the catalog.

Figure 10: Example lightcurves of 20 LPVs, with the periods and LPV
probabilites from our classifier. We see a good mix of large amplitude and
small amplitude LPVs, showing a variety of morphologies. The period and
internal designation of each LPV is listed. Lightcurves for all LPVs from the
catalog are publicly available. A few suspicious lightcurves can also be seen
in the above sample, which have a low classifier confidence (typically less
than 0.6).

Name | RA | Dec | num_dets | lcdur | amplitude | bestperiod | $\left<J\right>$ | prob_lpv
---|---|---|---|---|---|---|---|---
PGIRV_1008_0_3_357 | 172.238841 | -7.260807 | 83 | 962.6 | 0.07 | 130 | 9.41 | 0.86
PGIRV_1012_0_1_365 | 190.461586 | -7.169434 | 64 | 915.6 | 0.15 | 378 | 14.91 | 0.73
PGIRV_1012_2_0_355 | 192.904506 | -8.410068 | 68 | 915.6 | 0.18 | 409 | 13.03 | 0.50
PGIRV_1015_0_2_1656 | 205.447468 | -8.229639 | 69 | 892.6 | 0.12 | 155 | 12.87 | 0.96
PGIRV_1017_3_2_1044 | 218.119784 | -5.521623 | 51 | 885.6 | 1.26 | 276 | 13.98 | 0.93
PGIRV_1018_1_3_1424 | 220.276317 | -4.953930 | 66 | 904.6 | 0.05 | 240 | 11.37 | 0.52
PGIRV_1049_2_1_115 | 13.338281 | -11.976007 | 90 | 989.1 | 0.20 | 195 | 11.47 | 0.99
PGIRV_1054_0_1_194 | 35.383878 | -12.128101 | 96 | 838.8 | 0.07 | 173 | 10.12 | 0.70
PGIRV_1084_1_3_1287 | 183.570633 | -9.013661 | 83 | 934.6 | 0.60 | 159 | 13.22 | 0.92
PGIRV_1086_0_3_472 | 194.237602 | -11.771883 | 84 | 943.6 | 0.20 | 236 | 10.73 | 0.94
PGIRV_1091_3_1_1781 | 219.437017 | -9.584581 | 73 | 909.6 | 0.07 | 229 | 10.12 | 0.62
PGIRV_1119_1_0_0 | 354.895573 | -10.027897 | 71 | 975.2 | 0.07 | 188 | 10.13 | 0.52
PGIRV_1119_3_1_1517 | 357.465153 | -9.885638 | 90 | 991.1 | 0.15 | 353 | 11.44 | 0.51
PGIRV_1120_2_2_851 | 4.066668 | -17.603277 | 85 | 981.2 | 0.59 | 132 | 11.52 | 0.88
PGIRV_1121_3_0_1348 | 7.516286 | -15.770552 | 17 | 642.2 | 1.06 | 102 | 15.75 | 0.54
PGIRV_1122_3_1_364 | 13.373073 | -14.509416 | 88 | 989.2 | 0.27 | 175 | 10.64 | 0.99
PGIRV_1123_0_0_73 | 15.983818 | -17.514027 | 97 | 993.1 | 0.14 | 158 | 11.05 | 0.93
PGIRV_1125_0_3_507 | 26.950635 | -16.721510 | 42 | 784.8 | 0.22 | 241 | 16.53 | 0.55
PGIRV_1130_2_0_1225 | 52.649587 | -17.941429 | 88 | 846.7 | 0.16 | 103 | 15.87 | 0.68
PGIRV_1157_1_1_652 | 185.468635 | -13.886072 | 81 | 893.6 | 0.69 | 222 | 8.659 | 0.92
PGIRV_1190_1_0_597 | 350.546129 | -15.304388 | 86 | 989.2 | 0.37 | 159 | 10.86 | 0.90
PGIRV_1191_1_1_948 | 355.128686 | -13.872989 | 85 | 970.1 | 0.05 | 182 | 9.748 | 0.57
PGIRV_1192_2_0_573 | 3.142701 | -22.921270 | 89 | 988.2 | 0.78 | 170 | 10.39 | 0.97
PGIRV_1198_1_3_1166 | 32.186023 | -19.215745 | 10 | 751.9 | 0.80 | 256 | 15.05 | 0.55
PGIRV_1199_1_3_621 | 37.821288 | -19.131971 | 104 | 838.7 | 0.19 | 236 | 10.40 | 0.98
PGIRV_1199_2_3_1245 | 39.979909 | -21.072354 | 10 | 364.0 | 0.42 | 125 | 15.76 | 0.57
PGIRV_1256_2_3_445 | 333.906426 | -22.015062 | 13 | 994.2 | 1.09 | 361 | 16.20 | 0.84
PGIRV_1258_1_2_436 | 341.508786 | -20.476088 | 76 | 991.2 | 0.37 | 135 | 12.57 | 0.68
PGIRV_1260_1_1_910 | 349.897992 | -18.939772 | 78 | 974.2 | 0.81 | 151 | 11.50 | 0.94
PGIRV_1260_2_1_184 | 353.313861 | -21.525274 | 82 | 974.2 | 0.23 | 114 | 11.18 | 0.75
PGIRV_1263_3_3_824 | 9.703086 | -23.926492 | 76 | 786.8 | 0.59 | 177 | 11.21 | 0.99
PGIRV_1266_2_3_1063 | 25.300664 | -26.278947 | 41 | 788.8 | 0.24 | 232 | 16.48 | 0.58
PGIRV_126_1_1_1490 | 179.680170 | 63.354066 | 144 | 977.6 | 0.25 | 139 | 10.43 | 0.64
PGIRV_1296_1_0_1458 | 180.104200 | -24.686672 | 29 | 837.7 | 0.59 | 175 | 10.18 | 0.58
PGIRV_1296_1_1_58 | 181.156892 | -24.160596 | 80 | 876.6 | 0.27 | 139 | 11.24 | 0.90
PGIRV_1299_1_1_373 | 196.753528 | -24.091946 | 73 | 877.6 | 0.45 | 162 | 9.50 | 0.99
PGIRV_129_1_1_808 | 209.470908 | 64.111660 | 141 | 930.7 | 0.69 | 138 | 11.77 | 0.85
PGIRV_130_2_2_1558 | 225.541293 | 59.522531 | 145 | 932.7 | 0.32 | 373 | 8.530 | 0.71
PGIRV_1313_0_3_2444 | 271.707078 | -26.502290 | 34 | 735.9 | 0.25 | 149 | 11.52 | 0.87
PGIRV_1324_1_3_418 | 330.479227 | -23.974536 | 66 | 994.2 | 0.06 | 162 | 9.349 | 0.73
PGIRV_1325_0_1_761 | 334.291335 | -26.117588 | 68 | 987.2 | 0.88 | 341 | 11.04 | 0.93
PGIRV_1325_3_3_1070 | 337.938061 | -23.546005 | 69 | 987.2 | 0.75 | 209 | 11.48 | 0.92
PGIRV_1326_0_1_412 | 339.714107 | -26.934976 | 60 | 988.2 | 0.31 | 130 | 11.57 | 0.92

  * •

Description of columns :
$Name:$ PGIR name for variable; $RA$ and $Dec:$ RA and Dec (in degrees for
equinox 2000); $num\\_detections:$ Number of detections;
$lcdur:$ Observation baseline; $amplitude$ and $bestperiod:$ amplitude and
period; $\left<J\right>$: mean J band magnitude; $prob\\_lpv:$
Classifier confidence;
The complete catalog will be available as a machine-readable table.

Table 3: PGIR Long Period Variables

###  4.1 Catalog Features

Figure 9 shows the spatial distribution of LPVs. We observe that 155,220 of
LPVs are located between galactic latitudes ($b_{gal}$) of $-10^{\circ}$ to
$+10^{\circ}$. We observe that a large number of LPVs have low amplitudes
$(<0.25)$ mag and a low percentile difference, implying that they are likely
to be semi-regular variables, while the LPVs showing amplitudes greater than
0.5 mag are likely to be Miras. Figure 11 shows the distribution of the peak-
to-peak amplitude, the number of detections for each variable, the von Neumann
score ($\eta$), and the reduced $\chi^{2}$ of a sinusoid fit to a lightcurve.

The period histogram of LPVs is shown in Figure 12, with the Miras showing a
doubly peaked distribution and the semiregular variables showing a triply
peaked distribution. We obtain peaks at $\scriptstyle\mathtt{\sim}$200 d and
$\scriptstyle\mathtt{\sim}$340 d for the Miras, consistent with other studies
(Matsunaga et al., 2005, 2009; Glass et al., 2001; Nikzat et al., 2022), which
also find bimodal period distributions for Miras. The semiregular variables
show peaks in their period distribution at $\scriptstyle\mathtt{\sim}$130 d,
$\scriptstyle\mathtt{\sim}$220 d, and $\scriptstyle\mathtt{\sim}$330 d.
Conroy et al. (2015) introduced a period-amplitude relationship for Miras:

$\log_{10}{\Delta{J}}=\alpha\log_{10}P+\beta$ (11)

We examine this relation for the Mira variables from PGIR in the J-band as
shown in Figure 13 and compare it to the relation obtained in Nikzat et al.
(2022). While the general trend holds, there is a significant scatter around
the best-fit curve. We obtain $\alpha=+0.727\pm 0.007$ and $\beta=-1.923\pm
0.019$ from Miras in PGIR data. Our best-fit parameters agree reasonably with
Nikzat et al. (2022), with $\alpha=+0.76\pm 0.07$ and $\beta=-1.91\pm 0.19$.

Figure 11: Distribution of select features in the LPV catalog. The majority of
LPVs in the catalog (118,112) have amplitude less than 0.25 mag and are likely
to be semiregular variables. The von Neumann score shows a bimodal
distribution with peaks corresponding to Miras and semiregular variables.
Figure 12: Distribution of periods in the PGIR LPV catalog, separated by
amplitude, representing Miras (blue histogram) and semiregular variables (pink
histogram). We show the density in each bin. We observe a bimodal period
distribution for Miras with peaks at $\scriptstyle\mathtt{\sim}$200 d and
$\scriptstyle\mathtt{\sim}$340 d, while the semiregular variables show peaks
at $\scriptstyle\mathtt{\sim}$130 d, $\scriptstyle\mathtt{\sim}$220 d, and
$\scriptstyle\mathtt{\sim}$330 d.  Figure 13: Period Amplitude relationship of
Mira variables, chosen to be LPVs from the catalog having amplitude greater
than 0.4 mag. The white curve represents the best-fit parameters from PGIR
Miras, compared to the relation obtained by Nikzat et al. (2022). While the
general trend of the period-amplitude relation holds, there is a significant
scatter around the best-fit curve. The distribution of PGIR Miras in the
log(period)-amplitude space is shown as a 2D histogram with the colorbar
representing the number of objects in each bin.

###  4.2 Catalog Validation

The second catalog of long-period variables (Lebzelter et al., 2023) from Gaia
data release 3 is the largest database of LPVs to date, with 1,720,588 entries
spanning 34 months of Gaia data, with 391,914 LPVs having period information.
We validate our catalog by comparing it with 150,195 LPVs from Gaia, above a
declination of $-28^{\circ}$, that have calculated periods. We query the Gaia
archive to retrieve the positions, periods, colors, and distances (whenever
available) for these LPVs. We crossmatch the Gattini LPVs with Gaia keeping a
maximum offset of 6 arcsecs (corresponding to 2 PGIR pixels), to obtain 37,661
sources. Figure 14 shows the comparison between periods obtained from Gattini
and Gaia. We observe a strong positive correlation with some scatter around
the line of equal periods. We also observe a branch of correlated periods,
where the Gaia period is approximately twice as much as the Gattini period,
like the example shown in Figure 15. Visual inspection of these sources shows
that Gaia overestimates the period of these variables, likely due to sparser
coverage than PGIR.

Figure 14: Comparison between periods of crossmatched LPVs from PGIR and Gaia.
We observe a strong correlation between the two datasets. We also notice a
secondary branch, which approximately represents the line where the period
from Gaia is approximately twice the period from PGIR. The color gradient
represents the number of objects in each bin of the 2D histogram.  Figure 15:
Example of an LPV where period estimated from Gattini lightcurve is more
accurate than the one calculated from Gaia data

Additionally, we also crossmatch our catalog with SIMBAD to understand the
classification accuracy of the machine learning model. The SIMBAD website
provides a collection of almost all variables previously found in the
Catalina, OGLE, Gaia, and ASAS surveys; among others. We consider our LPV
classification to be correct for a star if its previously assigned category
from other surveys falls into one of ’LongPeriodV*’, ’AGB*’, ’C*’, ’Mira’,
’LongPeriodV*_Candidate’, ’C*_Candidate’, ’Variable*’, ’Star’, ’OH/IR*’,
’MidIR’, ’NearIR’, ’Mira_Candidate’ or ’RGB*’. Using a maximum crossmatch
radius of 6 arcsecs, we see that 73,346 (45.9% of the total LPVs) objects in
the PGIR catalog are new. Of the remaining 84,134 objects, 1,867 fall into
categories other than the ones mentioned before, indicating a high purity
(97.78%) of the PGIR catalog. The populations of the different Simbad
categories are listed in Table 4.

Category | Counts
---|---
New | 73,346
LongPeriodV*_Candidate | 40,195
LongPeriodV* | 22,386
Mira | 18,020
Star | 1,845
C* | 805
Variable* | 409
AGB* | 270
OH/IR* | 220
MidIR/NearIR | 127
C*_Candidate | 122
RGB* | 82
Mira_Candidate | 2
Others | 1,867

  * •

Objects with confirmed sub-types are counted into their respective sub-types
(e.g. Mira) and excluded from the broader class (e.g. LongPeriodV*)

Table 4: Simbad types of objects in the PGIR LPV catalog

##  5 Conclusion

We present a catalog of 159,696 Long Period Variables from the PGIR survey
created using a machine-learning classifier. We assemble a training set of
4300 objects spanning LPVs showing clean sinusoidal evolution in their
lightcurves, “Type II LPVs" showing non-sinusoidal but periodic lightcurves,
non-LPVs comprised of R Cor Bor stars, Young Stellar Objects and other classes
that display long timescale erratic variability, and non-variables. We train a
gradient-boosted decision tree classifier on an artificially resampled dataset
using Adaptive synthetic upsampling and allKNN downsampling to class-balance
the training set, achieving a true positive rate of 94.89% for the combined
set of LPVs and Type II LPVs and an overall weighted g-mean score of 0.95,
indicating high accuracy.

Following a feature extraction process for 35 million lightcurves from the
PGIR survey at a rate of 0.1 seconds per lightcurve, we generate the LPV
catalog using classifier predictions. We validate the catalog by crossmatching
it with the Gaia catalog of LPVs from the third data release (Lebzelter et
al., 2023) and the Simbad interface to ascertain purity and completeness. We
find that 73,346 variables are newly discovered, which is 45.9% of the total
catalog size, and find that 1,867 objects in the catalog fall out of desired
Simbad classification categories, indicating a purity of
$\scriptstyle\mathtt{\sim}$97.8%. We also find a good correlation in periods
between LPVs from PGIR and Gaia. We use the PGIR LPV catalog to examine the
period-amplitude relationship for Miras in the J band and verify the
bimodality in period distributions of Miras found in other studies (Matsunaga
et al., 2005, 2009; Glass et al., 2001; Nikzat et al., 2022). The methods
developed in this paper will be useful for more sensitive searches for
periodic variables with the Vera Rubin Observatory (Ivezić et al., 2019).

## Acknowledgements

Palomar Gattini-IR (PGIR) is generously funded by Caltech, Australian National
University, the Mt Cuba Foundation, the Heising Simons Foundation, the
Binational Science Foundation. PGIR is a collaborative project among Caltech,
Australian National University, University of New South Wales, Columbia
University and the Weizmann Institute of Science. MMK acknowledges generous
support from the David and Lucille Packard Foundation. MMK acknowledges the
US-Israel Bi-national Science Foundation Grant 2016227. MMK acknowledges the
Heising-Simons foundation for support via a Scialog fellowship of the Research
Corporation. MMK and AMM acknowledge the Mt Cuba foundation. JS is supported
by an Australian Government Research Training Program (RTP) Scholarship. AS
acknowledges support from the Caltech Summer Undergraduate Research Fellowship
(SURF) funded by the NSF Astronomy and Astrophysics grant number 2206730. KD
was supported by NASA through the NASA Hubble Fellowship grant #HST-
HF2-51477.001 awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. RS acknowledges grant number 12073029
from the National Natural Science Foundation of China (NSFC).

## References

  * Alcock et al. (2001) Alcock, C., Allsman, R. A., Alves, D. R., et al. 2001, ApJ, 554, 298
  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002
  * Bergstra & Bengio (2012) Bergstra, J., & Bengio, Y. 2012, J. Mach. Learn. Res., 13, 281–305
  * Boone (2019) Boone, K. 2019, AJ, 158, 257
  * Breiman (2001) Breiman, L. 2001, Machine Learning, 45, 5
  * Breiman et al. (1984) Breiman, L., Friedman, J. H., Olshen, R. A., & Stone, C. J. 1984, in Classification and Regression Trees
  * Chen et al. (2018) Chen, X., Wang, S., Deng, L., de Grijs, R., & Yang, M. 2018, ApJS, 237, 28
  * Chen et al. (2020) Chen, X., Wang, S., Deng, L., et al. 2020, VizieR Online Data Catalog, J/ApJS/249/18
  * Clayton (2012) Clayton, G. C. 2012, Journal of the American Association of Variable Star Observers (JAAVSO), 40, 539
  * Conroy et al. (2015) Conroy, C., van Dokkum, P. G., & Choi, J. 2015, Nature, 527, 488
  * Cutri & et al. (2012) Cutri, R. M., & et al. 2012, VizieR Online Data Catalog, II/311
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * De et al. (2020) De, K., Hankins, M. J., Kasliwal, M. M., et al. 2020, PASP, 132, 025001
  * Debosscher et al. (2007) Debosscher, J., Sarro, L. M., Aerts, C., et al. 2007, A&A, 475, 1159
  * Dekany et al. (2011) Dekany, I., Catelan, M., Minniti, D., & the VVV Collaboration. 2011, arXiv e-prints, arXiv:1111.0909
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220
  * Glass et al. (2001) Glass, I. S., Matsumoto, S., Carter, B. S., & Sekiguchi, K. 2001, MNRAS, 321, 77
  * He et al. (2008) He, H., Bai, Y., Garcia, E. A., & Li, S. 2008, in 2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence), 1322–1328
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Karambelkar et al. (2019) Karambelkar, V. R., Adams, S. M., Whitelock, P. A., et al. 2019, ApJ, 877, 110
  * Karambelkar et al. (2021) Karambelkar, V. R., Kasliwal, M. M., Tisserand, P., et al. 2021, The Astrophysical Journal, 910, 132
  * Kasliwal et al. (2018) Kasliwal, M., Jencson, J., Lau, R., et al. 2018, SPIRITS: SPitzer InfraRed Intensive Transients Survey, Spitzer Proposal ID #14089
  * Ke et al. (2017) Ke, G., Meng, Q., Finley, T., et al. 2017, in Advances in Neural Information Processing Systems, ed. I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, & R. Garnett, Vol. 30 (Curran Associates, Inc.)
  * Kim, Dae-Won & Bailer-Jones, Coryn A. L. (2016) Kim, Dae-Won, & Bailer-Jones, Coryn A. L. 2016, A&A, 587, A18
  * Lebzelter et al. (2023) Lebzelter, T., Mowlavi, N., Lecoeur-Taibi, I., et al. 2023, A&A, 674, A15
  * Lemaître et al. (2017) Lemaître, G., Nogueira, F., & Aridas, C. K. 2017, Journal of Machine Learning Research, 18, 1
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447
  * Masci et al. (2014) Masci, F. J., Hoffman, D. I., Grillmair, C. J., & Cutri, R. M. 2014, AJ, 148, 21
  * Masci et al. (2017) Masci, F. J., Laher, R. R., Rebbapragada, U. D., et al. 2017, PASP, 129, 014002
  * Matsunaga et al. (2005) Matsunaga, N., Fukushi, H., & Nakada, Y. 2005, MNRAS, 364, 117
  * Matsunaga et al. (2009) Matsunaga, N., Kawadu, T., Nishiyama, S., et al. 2009, MNRAS, 399, 1709
  * Nikzat et al. (2022) Nikzat, F., Ferreira Lopes, C. E., Catelan, M., et al. 2022, A&A, 660, A35
  * Pojmanski et al. (2005) Pojmanski, G., Pilecki, B., & Szczygiel, D. 2005, Acta Astron., 55, 275
  * Qu et al. (2021) Qu, H., Sako, M., Möller, A., & Doux, C. 2021, AJ, 162, 67
  * Samus’ et al. (2017) Samus’, N. N., Kazarovets, E. V., Durlevich, O. V., Kireeva, N. N., & Pastukhova, E. N. 2017, Astronomy Reports, 61, 80
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smak (1966) Smak, J. I. 1966, ARA&A, 4, 19
  * Soszyński et al. (2009) Soszyński, I., Udalski, A., Szymański, M. K., et al. 2009, Acta Astron., 59, 239
  * Stetson (1996) Stetson, P. B. 1996, PASP, 108, 851
  * Sánchez-Sáez et al. (2021) Sánchez-Sáez, P., Reyes, I., Valenzuela, C., et al. 2021, The Astronomical Journal, 161, 141
  * Tomek (1976) Tomek, I. 1976, IEEE Transactions on Systems, Man, and Cybernetics, SMC-6, 448
  * Udalski et al. (1992) Udalski, A., Szymanski, M., Kaluzny, J., Kubiak, M., & Mateo, M. 1992, Acta Astron., 42, 253
  * VanderPlas & Ivezić (2015) VanderPlas, J. T., & Ivezić, Ž. 2015, ApJ, 812, 18
  * von Neumann (1941) von Neumann, J. 1941, The Annals of Mathematical Statistics, 12, 367
  * Williams & Rasmussen (1995) Williams, C., & Rasmussen, C. 1995, in Advances in Neural Information Processing Systems, ed. D. Touretzky, M. Mozer, & M. Hasselmo, Vol. 8 (MIT Press)
  * Wood et al. (1999) Wood, P. R., Alcock, C., Allsman, R. A., et al. 1999, in Asymptotic Giant Branch Stars, ed. T. Le Bertre, A. Lebre, & C. Waelkens, Vol. 191, 151
  * Wray et al. (2004) Wray, J. J., Eyer, L., & Paczyński, B. 2004, MNRAS, 349, 1059
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868


# GPGPU Acceleration of the KAZE Image Feature Extraction Algorithm

Ramkumar B, R. S. Hegde  Department of Electrical Engineering
Indian Institute of Technology
Gandhinagar, Gujarat, India 382355
Email<EMAIL_ADDRESS>Rob Laber, Hristo Bojinov  Innit Inc.
Redwood City, CA 94063, USA

###### Abstract

The recently proposed open-source KAZE image feature detection and description
algorithm [1] offers unprecedented performance in comparison to conventional
ones like SIFT and SURF as it relies on nonlinear scale spaces instead of
Gaussian linear scale spaces. The improved performance, however, comes with a
significant computational cost limiting its use for many applications. We
report a GPGPU implementation of the KAZE algorithm without resorting to
binary descriptors for gaining speedup. For a 1920 by 1200 sized image our
Compute Unified Device Architecture (CUDA) C based GPU version took around 300
milliseconds on a NVIDIA GeForce GTX Titan X (Maxwell Architecture-GM200) card
in comparison to nearly 2400 milliseconds for a multithreaded CPU version (16
threaded Intel(R) Xeon(R) CPU E5-2650 processsor). The CUDA based parallel
implementation is described in detail with fine-grained comparison between the
GPU and CPU implementations. By achieving nearly 8 fold speedup without
performance degradation our work expands the applicability of the KAZE
algorithm. Additionally, the strategies described here can prove useful for
the GPU implementation of other nonlinear scale space based methods.

###### Index Terms:

Nonlinear scale space, Feature detection, Feature description, GPU, KAZE
features.

##  I Introduction

Feature point detection [2] and description [3] is a key tool for many
computer vision based applications, such as visual navigation [4], automatic
target recognition [5], tracking, structure from motion, registration,
calibration, and more. By picking out only the most salient points of an
image, that can be repeatably localized across different images, we can vastly
reduce subsequent data processing. Feature extraction [2], however, still
remains a major bottleneck for many implementations due to the high
computational cost, especially those that are most robust. GPGPUs have become
promising for image processing [6], computer vision [7, 6], and deep learning
tasks [8]. The GPU is an appropriate choice to process image related problems
because these involve large data sizes and exhibit a high arithmetic intensity
(the ratio between arithmetic operations and memory operations).

Scale Invariant Feature Transform (SIFT) [9, 10, 3] is widely considered as
one of the most robust feature descriptors that exhibits distinctiveness and
invariance to common image transformations. Some of the binary descriptors
have been developed to detect and match interest regions that can decrease the
computational cost. FAST [11] key-point detectors with Binary Robust
Independent Elementary Features (BRIEF) feature descriptor [12] have yielded
better results in real time applications. However, BRIEF and other such binary
descriptors are not very robust to image transformations. Binary Robust
Invariant Scalable Keypoints (BRISK) [13] and Oriented FAST and rotated BRIEF
(ORB) [14] methods have been developed by modifying BRIEF and FAST to achieve
scale and rotation invariance to some degree. The recently proposed KAZE [1]
algorithm outperforms SIFT and several other algorithms inspired by it. In the
short span of its introduction, it has already become widely used in image
matching [15], target classification [16], and data mining [17, 18].

Figure 1: Flowchart for KAZE Algorithm

The improved performance of the KAZE algorithm [3] however comes at the cost
of increased computational cost. The original KAZE algorithm provides modified
SURF (Speeded up robust features) like descriptors allowing it to be used as a
replacement for SIFT and SURF. An accelerated version (multithreaded CPU
version) of the original KAZE algorithm termed as A-KAZE [19] was proposed by
the original authors and reduces the computational cost at the expense of
algorithm performance in some scenarios. The A-KAZE uses binary descriptors in
place of SURF-like descriptors and thus may not easily replace SURF in
computer vision pipelines. Subsequently, the authors have reported on a GPU
based implementation of the A-KAZE [20]. Although this has yielded better
performance in terms of computational cost when compared to CPU based KAZE, a
compromise of algorithm performance in certain outcomes is again involved due
to it resorting to binary descriptors.

In this paper, we mainly focus on parallelizing the original version of KAZE
on the GPU (using Compute Unified Device Architecture (CUDA) C) without
resorting to binary descriptors. We primarily take advantage of the
characteristic of the algorithm and find out parallel strategies for each
stage of KAZE algorithm. Our implementation is primarily GPU based; the CPU
just plays a part in controlling the GPU, initializing CUDA kernels and data
allocation and is free to perform other tasks during the bulk of the
computation. The paper is organized as follows, following this introduction
the KAZE algorithm and its parts are described in detail in section II. In
section III, we have discussed the implementation of KAZE algorithm on GPU
using CUDA. The obtained results have been discussed in the section IV.
Section V concludes the paper.

##  II The KAZE algorithm

The KAZE Features [3] algorithm is a novel feature detection and description
method and it belongs to the class of methods which utilize the so-called
scale space. Its novelty arises in that it operates using a nonlinear scale
space whereas previous methods such as SIFT or SURF [21, 3] find features in
the Gaussian scale space (a particular instance of linear diffusion). The
three main steps involved in KAZE feature extraction algorithm are:

  1. 1.

Construct a Nonlinear Scale space pyramid of the original image.

  2. 2.

Determine keypoints using Hessian determinants and multiscale derivatives in
the nonlinear scale space.

  3. 3.

Compute orientation and descriptor vectors for all keypoints.

The detailed overview of these steps in the KAZE algorithm is given in Figure
1 and further explained below.

###  II-A Building the nonlinear scale space

The main idea for constructing scale space is to obtain a separable structures
of image from the original image, such that only fine scale image structures
exist in the multi-scale representation [2, 22] such as scale-space
representation, pyramids, and non-linear diffusion methods. Linear scale space
systems [22] (like the Gaussian) smear the edges such that small-scale
segmentations do not coincide with large-scale segmentations leading to
reduction in localization accuracy. Non linear scale space systems overcome
this by locally blurring the image data so that details remain unaffected and
noise gets blurred. But computing the nonlinear scale spaces are slow and
diffusion function should be chosen carefully .

Consider an image of size $W\times H$ (width and height of image
respectively). Anisotropic diffusion filter has been utilized while smoothing
in order to preserve the edges and is defined as:

$\frac{\mathrm{d}I}{\mathrm{d}t}=div\left(c(x,y,t).\nabla I\right),$ (1)

where $I$ is the image luminance, div is the divergence operator, $c(x,y,t)$
is the conductivity function which plays a key role in controlling the
diffusion of the image, $\nabla$ is the gradient operator and $t$ is the
timing scale parameter. Magnitude of the image gradient controls the
diffusion. The conductivity function is defined as follows:

$c(x,y,t)=g(|\nabla I_{\sigma}((x,y,t)|),$ (2)

where the function $\nabla I$ is the image gradient obtained after applying
Gaussian smoothing on the image. Equation 1 is generally known as Perona and
Malik diffusion equation [23]. The two different forms of conductivity
function ’$g$’ described by Perona and Malik are:

$g1=exp\left(-\frac{\mathrm{|}\nabla
I_{\sigma}|^{2}}{\mathrm{k}^{2}}\right),g2=\frac{\mathrm{1}}{\mathrm{1}+\frac{\mathrm{|}\nabla
I_{\sigma}|^{2}}{\mathrm{k}^{2}}},$ (3)

where ’$k$’ is the contrast factor that decides whether the edges to be
smoothed or filtered out.

We need to approximate the differential equations using numerical methods as
there are no analytical solutions for Equation 1. The simple solution is to
use explicit schemes but these are found to be impractical for feature
detection due to their computational complexity. To accelerate the nonlinear
scale space generation, Additive Operator Fast Explicit Diffusion (FED) [24,
25] scheme was introduced by [1]. The nonlinear scale space can be built
efficiently by means of such Fast Explicit Diffusion (FED) [25] schemes and
they are numerically stable for any step size.

The schemes discretize the Equation 1 as follows:

$\frac{I^{i+1}-I^{i}}{\tau}=\displaystyle\sum_{l=1}^{m}A_{l}(I^{i})I^{i+1},$
(4)

where $A_{l}$ encodes the image conductivities for each dimension, $\tau$ is a
constant time step to maintain stability. The main intention of FED schemes is
to perform $M$ cycles of $n$ explicit diffusion steps with varying step sizes
$\tau_{j}$ and is defined by:

$\tau_{j}=\frac{\tau_{max}}{2\cos^{2}\left(\pi\frac{2j+1}{4n+2}\right)}.$ (5)

Discretize the scale space in logarithmic steps that are arranged in a series
of $O$ octaves and $S$ sub-levels. We always preserve the resolution of
original image at each new octave whereas in SIFT we perform downsampling at
each new octave. Discrete values of octave index $o$ and a sub-level index $s$
are used to identify set of octaves and sub-levels respectively. The octave
and the sub-level indexes are mapped to their corresponding scale $\sigma$ as
defined by the following equation:

$\displaystyle\displaystyle\sigma_{i}(o,s)=\sigma_{0}2^{\frac{o+s}{S}},$ (6)
$\displaystyle o\in[0...O-1],s\in[0...S-1],$ $\displaystyle i\in[0...N],$

where $\sigma_{0}$ is the base scale level and $N$ is the total number of
filtered images.

Pixel units $\sigma_{i}$ in the set of discrete scale levels is converted to
time units because of nonlinear diffusion filtering which is defined in units
of time. In the Gaussian scale space, filtering the image for some time $t$ =
$\frac{\sigma^{2}}{2}$ is equivalent to convolution of the image with Gaussian
of standard deviation $\sigma$ (in pixels). The result of this conversion is
applied to transform the scale space $\sigma_{i}(o,s)$ in units of time and
obtain a set of evolution times by mapping $\sigma_{i}\rightarrow t_{i}$
defined by:

$t_{i}=\frac{1}{2}\sigma_{i}^{2},i={0...N}.$ (7)

Nonlinear scale space is constructed from the obtained set of evolution times
using $\sigma_{i}\rightarrow t_{i}$ mapping. For each filtered image $t_{i}$
in the non linear scale space, the convolution of the original image with a
Gaussian of standard deviation $\sigma_{i}$ does not correspond with the
resulting image.

###  II-B Feature Detection

Feature detection [26, 4] is the identification of interesting image
primitives (e.g. points, lines/curves, and regions) for the purpose of
highlighting salient visual cues in digital images. The primary goal of
feature detection is to effectively extract highly stable visual features.

In order to increase the detection accuracy, Hessian determinant is computed
for each filtered image $L_{i}$ in the nonlinear scale space.
$k_{i,norm}=\frac{k_{i}}{2^{o^{i}}}$ is a scaling factor to normalize the
computed Hessian determinant for the corresponding filtered images.

$L_{i,Hessian}=k_{i,norm}^{2}(L_{i,xx}L_{i,yy}-L_{i,xy}L_{i,xy}).$ (8)

Concatenated Scharr filter of step size $k_{i,norm}$ is used to compute second
order derivatives. When compared to other filters, Scharr filters [27] have
better central differences differentiation and rotation invariance. First, we
search for detector response maxima in spatial location. Preserve the maxima
in $3\times 3$ pixels window with detector responses higher than the threshold
for each step. We cross-check whether the response is maxima with respect to
other keypoints from level $i-1$ to $i+1$ for each of the response directly
below and above in a $\sigma_{i,s}\times\sigma_{i,s}$ pixels window. Hessian
determinant response in $3\times 3$ pixels neighborhood is fitted by the two-
dimensional quadratic function and its maximum value is used to estimate the
two-dimensional position of the keypoint with sub-pixel accuracy.

###  II-C Feature Description

A feature description is a process which takes an image with interest points
and yields feature descriptors (vectors). The obtained feature descriptors
acts as numerical "fingerprint" which differentiates one feature from another
and are ideally invariant to image transformations. There are several methods
to find feature descriptors like SIFT, SURF, Histogram Oriented Gradients
(HOG) [28]. Here we use modified SURF (M-SURF) [29] for computing descriptors.

####  II-C1 Finding the Dominant Orientation

Estimating the dominant orientation with keypoint location as a center in a
local neighborhood is a key step for obtaining rotation invariant descriptors.
We compute the dominant orientation with a sampling step of size $\sigma_{i}$
in a circular area of radius $6\sigma_{i}$. First order derivatives $L_{x}$
and $L_{y}$ are weighted with a Gaussian weighted at the keypoint location as
a center for each of the samples in the circular area. Then, represent the
derivative responses as points in vector space and sum the responses within a
sliding circle segment covering an angle of $\pi/3$ to find the dominant
orientation. Finally, dominant orientation is obtained from the longest
vector.

####  II-C2 Building the Descriptor

Frame work of nonlinear scale space is embedded with M-SURF for computing
descriptors. Compute the first order derivatives $L_{x}$ and $L_{y}$ of size
$\sigma_{i}$ over a $24\sigma_{i}\times 24\sigma_{i}$ rectangular window for a
detected feature at scale $\sigma_{i}$. This rectangular window is divided
into $4\times 4$ subregions with an overlap of $2\sigma_{i}$ and with the size
$9\sigma_{i}\times 9\sigma_{i}$. Gaussian weighted
($\sigma_{1}=2.5\sigma_{i}$) derivative responses in each subregion is
centered on the subregion center and summed into a descriptor vector
$d_{v}=(L_{x},L_{y},|L_{x}|,|L_{y}|)$. Then, each subregion vector is Gaussian
weighted ($\sigma_{2}=1.5\sigma_{i}$ ) over a mask of $4\times 4$ with
keypoint as center. Compute the the derivatives according to the dominant
orientation. In order to achieve the contrast invariance, normalize the
descriptor vector of length 64 into a unit vector.

##  III CUDA Implementation of the KAZE algorithm

While the KAZE algorithm has demonstrated superior performance in terms of
feature robustness in relation to the class of linear scale space based
methods, the computations involving the nonlinear scale space considerably add
to the computational cost. However, a lot of parallelism exists in the
algorithm and it is thus suitable for a GPU based acceleration. We describe
the CUDA based implementation details for all the three main steps of the KAZE
algorithm in this section. As long as they are available we always use
standard OpenCV GPU primitives.

We have used Gpumat object of OpenCV to allocate memory in GPU global memory.
We allocate memory to evolution level with dimensions $O\times S$ (number of
octaves and number of sublevels respectively) for each evolution image of
nonlinear scale space pyramid ($L_{t}$), first order derivatives ($L_{x}$ and
$L_{y}$), second order derivatives ($L_{xx}$, $L_{yy}$ and $L_{xy}$), images
obtained after Gaussian smoothing ( $L_{smooth}$) and hessian determinants
($L_{det}$) of size $W\times H$.

####  III-1 Nonlinear scale space generation

For a given input image, in order to reduce image artifacts and noise, we
convolve the image with a Gaussian kernel of standard deviation $\sigma_{0}$.
The resultant image is then considered as the base image and its image
gradient histogram is computed to obtain the contrast parameter $k$ in an
automatic procedure. For the obtained set of evolution times $t_{i}$ and
contrast parameter, we then build the nonlinear scale space using the FED
scheme. At each evolution level of the nonlinear scale space, several
auxiliary images are generated. The image $L_{smooth}$ is the image obtained
by performing a two-dimensional Gaussian convolution using OpenCV CUDA
primitive kernel on the image generated at the end of the previous evolution
level (or on the base image in the case of the first level). $L_{smooth}$ is
calculated by allocating a block of $8\times 8$ threads and a grid of
$(W+8)/8\times(H+8)/8$ blocks. First order $x$ and $y$ Gaussian derivatives of
$L_{smooth}$ in the form of images $L_{x}$ and $L_{y}$ are then calculated
using Scharr filters [27]. We determine the diffused image $L_{flow}$ at each
level by then using our GPU kernel to compute Perona malik conductivity
equation from the first order Gaussian derivatives $L_{x}$ and $L_{y}$. We
allocate a block of $32\times 32$ threads and a grid of
$(W+32)/32\times(H+32)/32$ blocks and then each thread computes the
conductivity equation from the first order Gaussian derivatives on each pixel
of an image to form the flow image. From the obtained flow images, we need to
determine step images. So, we implemented Fast Explicit Diffusion that
executes over several steps (for defined step size) using a GPU kernel that
yields step image at each evolution level. We have allocated a block
containing $32\times 32$ threads and a grid of $(W+32)/32\times(H+32)/32$
blocks so that each thread computes Fast Explicit Diffusion Scheme (FED) using
diffusivity factor of image over the inner steps at each scale to build
nonlinear scale space. The obtained step image after applying FED scheme is
considered as the image of the corresponding evolution level. The performance
of nonlinear scale space generation kernel has been improved using shared
memory access when appropriate.

####  III-2 Implementation of Feature Points Detection

In order to determine the detector response at each evolution level, we
calculate the multiscale derivatives. To find the first order and second order
derivatives at each evolution level for larger scale size, we use separable
larger linear row filter and the larger linear column filter kernels and
combine them for calculating Scharr derivatives. We have allocated a block
containing $64\times 4$ threads and a grid of $(W+64)/64\times(H+4)/4$ blocks
for both row filter and column filters so that each thread computes the
derivative response of each pixel of the image at each evolution level of
defined kernel window with scale size $\sigma_{i}$. We combine the results of
both the filter to form large and separable linear two dimensional filter
which yields multiscale derivative responses of the image at each evolution
level with scale size $\sigma_{i}$. From the obtained second order derivatives
$L_{xx}$,$L_{yy}$ and $L_{xy}$, we calculate detector response of images at
each evolution level. Although we compute multiscale derivatives for every
pixel in the feature detection step, we also utilize the same set of computed
derivatives in the feature description step in order to reduce computational
cost. We allocate $8\times 8$ threads in a two-dimensional block. First, the
current and two adjacent scales in non-linear space should be bound to three
texture references. Each thread takes out 27 corresponding pixels in three
levels by global fetch. Then, for each of the potential response, we check
that the response is a maxima with respect to other keypoints from level $i-1$
and $i+1$, respectively directly above and directly below in a window of size
$\sigma_{i,s}\times\sigma_{i,s}$ pixels. The central point in these pixels is
compared to others to determine whether it is the extreme point. If a pixel is
an extreme point, it is selected as a candidate keypoint. Last, the edge
responses must be eliminated. To eliminate edge points, the principal
curvature around the surface $D(x,y)$ of a candidate point at (x,y) can be
calculated using the Hessian matrix of candidate keypoints. In practice, the
image pixels are discrete, so the Hessian matrix can be calculated by 8 points
around those points. The Hessian matrix is given by

$H=\begin{pmatrix}D_{xx}&D_{xy}\\\ D_{xy}&D_{yy}\\\ \end{pmatrix}$ (9)

Where $D_{xx}$, $D_{xy}$, and $D_{yy}$ are second order local derivative of
candidate keypoint. Let the largest eigenvalue be $\alpha$, and the smallest
eigenvalue is $\beta$. Then the sum of the eigenvalues from the trace of H and
their product from the determinant can be computed as:

$Tr(H)=D_{xx}+D_{xy}=\alpha+\beta$ (10) $Det(H)=D_{xx}D_{yy}-D_{xy}^{2}$ (11)

Let r be the ratio between the $\alpha$ and $\beta$, so $\alpha=r\beta$. Then,

$\frac{Tr(H)^{2}}{Det(H)}=\frac{(\alpha+\beta)^{2}}{\alpha\beta}=\frac{(r+1)^{2}}{r},$
(12)

Figure 2: Runtimes of multithreaded CPU-KAZE and GPU-KAZE Figure 3: Runtimes
of CPU-KAZE and GPU-KAZE for different scales

###  III-A Implementation of Feature Points Description

####  III-A1 Implementation of Orientation Assignment

To build the orientation, we divide the circle into segments each covering an
angle of $\pi/3$. We have defined two dimensional thread block with 121
threads ($11\times 11$) such that each thread computes the orientation of a
keypoint. Information of all neighborhood pixels around a keypoint in a
circular region with radius of $6\sigma$ ($\sigma$ is 1.5 times as the scale
of the keypoint) is required to compute orientation. The information of
keypoint and their relative pixels that have been computed is transferred to
each thread using shared memory. We compute Gaussian weighted first order
derivatives centered at the keypoint for each sample in the segments of a
circular area. We represent the corresponding derivative responses as points
in vector space and is stored in GPU memory space. Finally, the computed
responses are summed over the sliding segments of the circle over an angle
$\pi/3$. The maximum response in the longest vector is the dominant
orientation.

####  III-A2 Implementation of Feature Points Descriptor

Computation of MSURF-descriptors with the non-linear space is parallelized. We
use shared memory to accelerate the computation of descriptors and store the
derivatives responses over the squared regions. First order derivatives
$L_{x}$ and $L_{y}$ are computed for detected feature at each scale
$\sigma_{i}$ over square region of $24\sigma_{i}\times 24\sigma_{i}$ around
keypoint. We divide each region into $4\times 4$ square sub-regions whose edge
are $9\sigma_{i}$. Compute the Gaussian weighted ($\sigma$ is 2.5 times as the
scale of the keypoint) derivative responses in each subregion around its
center and summed into descriptor vector. Compute Gaussian weighted ($\sigma$
is 1.5 times as the scale of the keypoint) subregion vector around the
keypoint over the mask of $4\times 4$ region. A 4-orientation derivative
response is generated by calculating the contribution of the orientation of
each pixel to the orientation in a sub-region. So we can obtain $4\times
4\times 4$ responses to form a 64-dimensional vector. We use one thread to
process a sub-region. In our allocation strategy, there are 16 threads in a
block to process 4 keypoints. A thread computes the weight of all the pixels
in the sub-region and transfer responses to shared memory. Then the processing
results of the 16 threads could generate a 64-dimensional feature vector.

The number of keypoints decreases rapidly as the scale size $\sigma_{i}$
gradually increases so does the number of threads which are used to process
and compute information of keypoints. A preprocessing before orientation
assignment and keypoints descriptor is put forward. In this preprocess, some
information of images such as size, scale etc., and the memory address of
feature vectors are calculated by CPU and stored in GPU global memory. Because
these data would not be accessed frequently and never modified, they are
stored in the constant memory which saves bandwidth and accelerates the
accessing speed. Moreover, all the images and information of keypoints should
be bound to several texture references. After pre-processing, the kernel will
be initialized.

To improve load balance, we adapted the following process.

  1. 1.

For a given keypoint, the location of the keypoint in image pyramid is
calculated by each thread.

  2. 2.

To compute keypoint descriptor vectors and their corresponding orientations,
we should reuse the information stored in constant memory.

  3. 3.

the computed feature vectors should be stored in the appointed GPU memory
addresses.

In contrast to methods by which the scale space levels are processed one by
one, all the image scales are computed at the same time to make full use of
threads that are potentially starving. We reallocate the blocks and threads
which process the keypoints in different scale images concurrently.
Consequently, this allocation strategy ameliorates the load imbalance.

##  IV Results and Discussion

Figure 4: Individual runtimes of key steps in CPU-KAZE and GPU-KAZE Figure 5:
Runtimes of multithreaded CPU-KAZE and GPU-KAZE for key tasks with respect to
the image dimensions of $1920\times 1200$ and $480\times 640$ respectively

In this paper, we implement the program on the GPU NVIDIA Geforce GTX TITAN X
(MAXWELL). It has 12 GB of global memory and 64 KB of constant memory. This
GPU has 48KB shared memory and at most 65536 registers on each Streaming
Multiprocessors (SM). Moreover, This GPU also has 32 blocks and 64 warps with
32 threads per warp on each SM. The CPU version of KAZE is implemented on
Intel(R) Xeon(R) CPU E5-2650 processor configured with 8 cores and 16 threads.

In Figure 2(a), we plot the runtimes for both the multithreaded CPU-KAZE with
$m=16$ (where ’$m$’ denotes the number of threads) and our GPU-KAZE
implementation for various image sizes. While the runtimes increase as a
function of the image size, the rate of increase is significantly larger for
the CPU version. The speedup of the GPU version is thus seen to improve with
increasing image size. GPU-KAZE is around 10 times faster than the CPU-KAZE
for the image of size $1920\times 1200$. In Figure 2(b), we have plotted the
runtimes of CPU-KAZE with varying number of threads ($(m)=4,8,12,16$) with
respect to increasing image dimensions. It is seen that the speedup does not
scale well with increased thread count. The GPU version, however, will scale
well with the number of CUDA cores in the GPU chip.

For a given image size the number of keypoints will change depending on the
complexity of the image and thus runtimes for images of the same size may
change. In order to study the effect of image complexity, we created an image
dataset of 8 different image dimensions ranging from $240\times 400$ (lower
image dimension) to $1920\times 1200$ (Higher image dimension) with 6 images
of varying complexity in each dimension. In Figure 3, we have plotted the
average time taken by CPU-KAZE and GPU-KAZE for each image in this dataset
grouped by their dimension. Although there is now variance in the runtime for
an image of a given size, the speedup factors observed earlier still seem to
hold.

In Figure 4, we assess the runtimes with respect to the individual key steps
involved in the KAZE algorithm. We can notice almost 10 times faster
performance of GPU-KAZE than CPU-KAZE for each step as the image dimension
increases. We can also observe that building non-linear scale space consumes
most of the time when compared to feature detection (computes multiscale
derivatives and keypoints) and feature description (computes descriptors).

Within these three major steps, construction of the Nonlinear scale space
pyramid, computation of the keypoints orientation and descriptor calculation
are the substeps with the largest computational burden. The runtimes of GPU-
KAZE and multithreaded CPU-KAZE (with variation in the number of threads
$(m)=1,4,16$) for these subtasks are plotted in Figure 5. Two different image
dimensions were considered. It is clearly seen that better speedups are
obtained especially for the nonlinear pyramid construction step for the image
with larger dimensions. The speedup factor for the smaller and large dimension
is listed in Table I and Table II respectively.

Figure 6: Memory foot print of CPU-KAZE and GPU-KAZE for different scales
Figure 7: Keypoints obtained from CPU-KAZE and GPU-KAZE

Efficient memory optimization is also a key factor that affects the
performance capabilities of any algorithm. The impact on performance by
effective utilization, optimization, and reuse of heap space can be observed
as image dimensions increases or when a huge number of keypoints are
extracted. In Figure 6, we plot the memory utilization of CPU-KAZE that has
been analyzed using Valgrind tool whereas GPU-KAZE is analyzed using NVIDIA
profiler. Memory utilized by CPU-KAZE is shown in Figure 6(a). In Figure 6(b),
we also plot memory utilized by GPU-KAZE and memory utilized by scratch images
say $L_{step}$ (represented by the green colored bar chart) with varying image
dimensions. We can notice that the memory utilized by CPU-KAZE for the larger
image dimension ($1200\times 1920$) is around twenty times higher than the
memory utilized by CPU-KAZE for the lower image dimension ($480\times 640$)
whereas memory utilized by GPU-KAZE for the larger image dimension
($1200\times 1920$) is around three times higher than the memory utilized by
GPU-KAZE for the lower image dimension ($480\times 640$) implying better
memory optimization in the case of GPU-KAZE when compared with CPU-KAZE.

Table I: GPU-KAZE over multithreaded CPU- ($m=$ denotes number of CPU threads) speedup factor for the image size $480\times 640$ | $m=1$ | $m=4$ | $m=16$
---|---|---|---
Construct nonlinear scale space | 5 | 3.5 | 3.5
Compute feature descriptors | 130 | 60 | 20
Compute keypoint orientations | 200 | 75 | 20
Table II: GPU-KAZE over multithreaded CPU-KAZE ($m=$ denotes number of CPU threads) speedup factor for the image size $1200\times 1920$ | $m=1$ | $m=4$ | $m=16$
---|---|---|---
Construct nonlinear scale space | 18 | 16 | 13
Compute feature descriptors | 120 | 45 | 18
Compute keypoint orientations | 200 | 65 | 20

In Figure 7, we compare the extracted keypoints of CPU-KAZE and GPU-KAZE
implementations. The CPU-KAZE implements an approximate serial procedure by
first checking for maxima on a $3\times 3$ window instead of the
$\sigma\times\sigma$ window. We have implemented both this approximate
procedure and also the exact procedure. In the GPU version, no runtime penalty
was observed by implementing the exact procedure instead of the approximate
one. We have considered an image with dimension $480\times 640$ to check the
accuracy of the keypoint extraction for the GPU-KAZE as compared to the CPU-
KAZE. We have obtained a similar number of keypoints (i.e., around 2200
keypoints are extracted) for both CPU-KAZE and GPU-KAZE. We have noticed
relatively less number of keypoints when using the exact procedure in
comparison to the number of keypoints obtained using CPU-KAZE approximate
procedure.

##  V Conclusion

The computational cost for building non-linear scale space, feature
extraction, and feature description of GPU-KAZE are in the ratio $4:3:1$
respectively. From the keen analysis of the results obtained from GPU-KAZE, we
can notice that building nonlinear scale space and feature detection are the
key tasks that are computationally intensive. To improve the nonlinear scale
pyramid construction, new approaches to speed the solution of the Perona Malik
PDE are needed. We are investigating schemes that have been utilized to
perform nonlinear diffusion filtering [30, 31, 32] and can solve Partial
Differential Equations (PDE) [33] more efficiently. The feature detection
computations rely on spatial derivative calculation in a three dimensional
space. Utilizing the spatial locality features in GPU texture memory is one
way to improve these.

## Acknowledgment

RB and RSH acknowledge funding support from Innit Inc. Consultancy grant
CNS/INNIT/EE/P0210/1617/007 and High Performance Computing Lab support from
Mr. Sudeep Banerjee.

## References

  * [1] P. Alcantarilla, A. J.Davison, and A. Bartoli, “Kaze Features,” _Procedings of the British Machine Vision Conference_ , vol. 7577 LNCS, no. 6, pp. 13.1–13.11, 2013.
  * [2] T. Lindeberg, “Feature Detection with Automatic Scale Selection,” _International Journal of Computer Vision_ , vol. 30, no. 2, pp. 79 – 116, 1998.
  * [3] S. Bianco, D. Mazzini, D. P. Pau, and R. Schettini, “Local detectors and compact descriptors for visual search: A quantitative comparison,” _Digital Signal Processing: A Review Journal_ , vol. 44, no. 1, pp. 1–13, 2015.
  * [4] S. Gauglitz, T. Hollerer, and M. Turk, “Evaluation of interest point detectors and feature descriptors for visual tracking,” _International Journal of Computer Vision_ , vol. 94, no. 3, pp. 335–360, 2011.
  * [5] A. Sanna and F. Lamberti, “Advances in target detection and tracking in forward-looking infrared (FLIR) imagery,” _Sensors (Switzerland)_ , vol. 14, no. 11, pp. 20 297–20 303, 2014.
  * [6] J. Fung and S. Mann, “Using graphics devices in reverse: GPU-based Image Processing and Computer Vision,” _2008 IEEE International Conference on Multimedia and Expo, ICME 2008 - Proceedings_ , pp. 9–12, 2008.
  * [7] I. P. Seung, S. P. Ponce, J. Huang, Y. Cao, and F. Quek, “Low-cost, high-speed computer vision using NVIDIA’s CUDA architecture,” _Proceedings - Applied Imagery Pattern Recognition Workshop_ , 2008.
  * [8] A. Coates, B. Huval, T. Wang, D. Wu, and A. Y. Ng, “Deep learning with COTS HPC systems,” _Proceedings of The 30th International Conference on Machine Learning_ , pp. 1337–1345, 2013.
  * [9] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International Journal of Computer Vision_ , vol. 60, no. 2, pp. 91–110, 2004\.
  * [10] K. Mikolajczyk, K. Mikolajczyk, C. Schmid, and C. Schmid, “A performance evaluation of local descriptors,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 27, no. 10, pp. 1615–1630, 2005.
  * [11] E. Rosten and T. Drummond, “Machine learning for high-speed corner detection,” _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , vol. 3951 LNCS, pp. 430–443, 2006.
  * [12] M. Calonder, V. Lepetit, C. Strecha, and P. Fua, “BRIEF: Binary robust independent elementary features,” _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , vol. 6314 LNCS, no. PART 4, pp. 778–792, 2010.
  * [13] S. Leutenegger, M. Chli, and R. Y. Siegwart, “BRISK: Binary Robust invariant scalable keypoints,” _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 2548–2555, 2011.
  * [14] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “ORB: An efficient alternative to SIFT or SURF,” _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 2564–2571, 2011.
  * [15] Y. Lehiani, M. Preda, M. Maidi, and F. Ghorbel, “Object identification and tracking for steady registration in mobile augmented reality,” _IEEE 2015 International Conference on Signal and Image Processing Applications, ICSIPA 2015 - Proceedings_ , pp. 54–59, 2016.
  * [16] X. Gao, W. Li, M. Loomes, and L. Wang, “A fused deep learning architecture for viewpoint classification of echocardiography,” _Information Fusion_ , vol. 36, pp. 103–113, 2017. [Online]. Available: http://dx.doi.org/10.1016/j.inffus.2016.11.007
  * [17] A. Camargo, D. Papadopoulou, Z. Spyropoulou, K. Vlachonasios, J. H. Doonan, and A. P. Gay, “Objective definition of rosette shape variation using a combined computer vision and data mining approach,” _PLoS ONE_ , vol. 9, no. 5, 2014.
  * [18] Y. Zhai, Y. S. Ong, and I. W. Tsang, “The emerging Big dimensionality,” _IEEE Computational Intelligence Magazine_ , vol. 9, no. 3, pp. 14–26, 2014\.
  * [19] P. Alcantarilla, J. Nuevo, and A. Bartoli, “Fast Explicit Diffusion for Accelerated Features in Nonlinear Scale Spaces,” _Procedings of the British Machine Vision Conference 2013_ , pp. 13.1–13.11, 2013.
  * [20] A. Pieropan, M. Björkman, N. Bergström, and D. Kragic, “Feature Descriptors for Tracking by Detection: a Benchmark,” 2016. [Online]. Available: http://arxiv.org/abs/1607.06178
  * [21] H. Bay, A. Ess, T. Tuytelaars, and L. Van Gool, “Speeded-Up Robust Features (SURF),” _Computer Vision and Image Understanding_ , vol. 110, no. 3, pp. 346–359, 2008.
  * [22] R. W. Harvey, A. Bosson, and J. A. Bangham, “A Comparison of Linear and Non-Linear Scale-Space Filters in Noise,” _Signal Processing VIII_ , vol. 1, no. 6, pp. 1777–1781, 1996.
  * [23] “Scale-Space and Edge Detection Using Anisotropic Diffusion,” pp. 629–639, 1990\.
  * [24] L. Feng, Z. Wu, and X. Long, “Fast Image Diffusion for Feature Detection and Description,” _International Journal of Computer Theory and Engineering_ , vol. 8, no. 1, pp. 58–62, 2016.
  * [25] S. Grewenig, J. Weickert, and A. Bruhn, “From box filtering to fast explicit diffusion,” _In: Goesele M., Roth S., Kuijper A., Schiele B., Schindler K. (eds) Pattern Recognition. DAGM 2010. Lecture Notes in Computer Science_ , vol. 6376 LNCS, pp. 533–542, 2010.
  * [26] Y. Li, S. Wang, Q. Tian, and X. Ding, “A survey of recent advances in visual feature detection,” _Neurocomputing_ , vol. 149, no. PB, pp. 736–751, 2015\.
  * [27] J. Weickert and H. Scharr, “A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance,” _Journal of Visual Communication and Image Representation_ , vol. 13, no. 1-2, pp. 103–118, 2002\. [Online]. Available: http://linkinghub.elsevier.com/retrieve/pii/S104732030190495X
  * [28] N. Dalal and B. Triggs, “Histograms of oriented gradients for human detection,” _Proceedings - 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, CVPR 2005_ , vol. I, pp. 886–893, 2005\.
  * [29] H. Huang, L. Lu, B. Yan, and J. Chen, “A new scale invariant feature detector and modified SURF descriptor,” _Proceedings - 2010 6th International Conference on Natural Computation, ICNC 2010_ , vol. 7, pp. 3734–3738, 2010.
  * [30] Y. Chen and T. Pock, “Trainable Nonlinear Reaction Diffusion: A Flexible Framework for Fast and Effective Image Restoration,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 39, no. 6, pp. 1256–1272, 2015.
  * [31] W. Hu, R. Hu, N. Xie, H. Ling, and S. Maybank, “Image Classification Using Multiscale Information Fusion Based on Saliency Driven Nonlinear Diffusion Filtering,” _Image Processing, IEEE Transactions on_ , vol. 23, no. 4, pp. 1513–1526, 2014.
  * [32] B. Chen, X.-H. Zhou, L.-W. Zhang, J. Wang, W.-Q. Zhang, and C. Zhang, “A new nonlinear diffusion equation model for noisy image segmentation,” _Advances in Mathematical Physics_ , vol. 2016, 2016.
  * [33] N. Katzourakis, “Generalised solutions for fully nonlinear PDE systems and existence–uniqueness theorems,” _Journal of Differential Equations_ , vol. 263, no. 1, pp. 641–686, 2017. [Online]. Available: http://dx.doi.org/10.1016/j.jde.2017.02.048