cross section? What is the magnetic field strength? There is a strong need for
realistic and detailed models for the amplification and twisting of the field
lines, plus some physics criterion for knowing when the field lines will
reconnect. (The need for better micro-physics theory actually applies to
reconnection events with solar flares, ordinary flare stars, and for
Superflare stars all across the HR diagram.) And we need to know why the many
other CVs, especially the ‘sister-RN’ U Sco, do not have any known Superflares
despite very intensive observations over many decades. And we need some sort
of understanding as to why V2487 Oph is so extreme amongst all other
Superflare stars.

##  8 Acknowledgements

We thank the allocation committees at McDonald Observatory, CTIO, and the
Kepler K2 mission for granting us time to observe V2487 Oph. We thank the many
AAVSO observers who tirelessly produced professional quality photometry,
despite no pay and no specific knowledge of how their photometry will be used.
We thank the PTF, Pan-STARRS, and ZTF teams for producing their light curves,
with the ZTF flare providing a strong proof that the Superflares are intrinsic
to the old nova system.

This paper includes data collected by the Kepler mission and obtained from the
MAST data archive at the Space Telescope Science Institute (STScI). Funding
for the Kepler mission is provided by the NASA Science Mission Directorate.
STScI is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5-26555. This research made use of
Astropy (Astropy Collaboration et al., 2018), Lightkurve, a Python package for
Kepler and TESS data analysis (Lightkurve Collaboration et al., 2018), and
EVEREST (Luger et al., 2016, 2018).

Partial financial support for BES and AP’s early work was provided by the
National Science Foundation, under grant AST-0708079. Partial financial
support for AP was provided by the Oak Ridge National Laboratory’s Ralph E.
Powe Junior Faculty Enhancement Award. Financial support for SZ was provided
by the College of Charleston’s School of Science and Mathematics Summer
Research Award. AP would additionally like to thank the LWS Slack Channel and
MS in particular for many useful co-writing sessions.

##  9 Data Availability

All of the photometry data used in this paper are explicitly given in Tables
2, 3, and 4, and are publicly available from the references and links in Table
1.

## References

  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Balona (2012) Balona L. A., 2012, MNRAS, 423, 3420
  * Balona et al. (2015) Balona L. A. et al., 2015, MNRAS, 450, 956
  * Bellm et al. (2019) Bellm E. C. et al., 2019, PASP, 131, 18002
  * Carrington (1859) Carrington R. C., 1859, MNRAS, 20, 13
  * Chambers ert al. (2016) Chambers K. C. et al., 2016, arXiv:1612.05560
  * de Laverny et al. (1998) de Laverny P., Mennessier M., Mignard F., Mattei J., 1998, A&A, 330, 169
  * Filippenko et al. (1998) Filippenko A. V., Leonard D. C., Modjaz M., 1998, IAU Circ., 6943
  * Frank et al. (2002) Frank J., King A., Raine D., 2002, Accretion Power in Astrophysics. Cambridge Univ. Press, Cambridge
  * Fryxell and Woosley (1982) Fryxell B. A., Woosley S. E., 1982, ApJ, 261, 332
  * Gunn et al. (1980) Gunn A. G., Migenes V., Doyle J. G., Spencer R. E., Mathioudakis M., 1997, MNRAS, 287, 199
  * Hachisu and Kato (2018) Hachisu I., Kato M., 2018, ApJS, 237, 4
  * Hachisu et al. (2002) Hachisu I., Kato M., Kato T., Matsumoto K., 2002, in Gänsicke B. T., Beuermann K., Reinsch K., eds, ASP Conf. Ser. Vol. 261, The Physics of Cataclysmic Variables and Related Objects. Astron. Soc. Pac., San Francisco, p. 629
  * Haisch et al. (1991) Haisch B. M., Strong K. T., Rodono M., 1991, ARA&A, 29, 275
  * Hernanz and Sala (2002) Hernanz M., Sala, G., 2002, Science, 298, 393
  * Hesse and Cassak (2020) Hesse M., Cassak P. A., 2020, JGR Space Physics, 125, e2018JA025935
  * Howard et al. (2018) Howard W. S. et al., 2018, ApJ, 860, L30
  * Howell et al. (2014) Howell S. B. et al., 2014, PASP, 126, 398
  * Janvier (2017) Janvier M., 2017, J. Plasma Phys., 83, 535830101, see arXiv:1612.06513
  * Kinemuchi et al. (2012) Kinemuchi K., Barclay T., Fanelli M., Pepper J., Still M., Howell S. B., 2012, PASP, 124, 963
  * Kretzschmar (2011) Kretzschmar M., 2011, A&A, 530, A84
  * Kunieda et al. (1984) Kunieda H. et al., 1984, PASJ, 36, 215
  * Lewin et al. (1993) Lewin W. H. G., van Paradijs J., Taam R. E., 1993, Space Sci. Rev., 62, 223
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration et al., 2018, ascl.soft. ascl:1812.013
  * Lingam and Loeb (2017) Lingam M., Loeb A., 2017, ApJ, 848, L41
  * Luger et al. (2016) Luger R., Agol E., Kruse E., Barnes R., Becker A., Foreman-Mackey D., Deming D., 2016, AJ, 152, 100
  * Luger et al. (2018) Luger R., Kruse E., Foreman-Mackey D., Agol E., Saunders N., 2018, AJ, 156, 99
  * Luri et al. (2018) Luri, X. et al., 2018, A&A, 616, A9
  * Lynch et al. (2000) Lynch D. K., Rudy R. J., Mazuk S., Puetter R. C., 2000, ApJ, 541, 791
  * Maehara et al. (2012) Maehara H. et al., 2012, Nature, 485, 478
  * Maehara et al. (2015) Maehara H. et al., 2015, Earth Planets & Space, 67, 59
  * Nomoto (1982) Nomoto K., 1982, ApJ, 253, 798
  * Pagnotta and Schaefer (2014) Pagnotta A., Schaefer B. E., 2014, ApJ, 788, 164
  * Pagnotta et al. (2009) Pagnotta A., Schaefer B. E., Xiao L., Collazzi A. C., Kroll P., 2009, AJ, 138, 1230
  * Patterson et al. (2021) Patterson J. et al., 2021, ApJ, in press, see arXiv:2010.07812
  * Paudel et al. (2020) Paudel R. R., Gizis J. E., Mullan D. J., Schmidt S. J., Burgasser A. J., Williams P. K. G. 2020, MNRAS, 494, 5751
  * Ragan (2009) Rau A. et al., 2009, PASP, 121, 1395
  * Rodono (1992) Rodono M., 1992, in Kondo Y., Sistero R. F., Polidan R. S., eds, IAU Symposium 151, Evolutionary Processes in Interacting Binary Stars. Kluwer, Boston, p. 71
  * Rubenstein and Schaefer (2000) Rubenstein E. P., Schaefer B. E., 2000, ApJ, 529, 1031
  * Schaefer-normals (1989) Schaefer B. E., 1989, ApJ, 337, 927
  * Schaefer-faints (1990) Schaefer B. E., 1990, ApJ, 353, L25
  * Schaefer-miras (1991) Schaefer B. E., 1991, ApJ, 366, L39
  * Schaefer (2005) Schaefer B. E., 2005, ApJ, 621, L53
  * Schaefer (2010) Schaefer B. E., 2010, ApJS, 187, 275
  * Schaefer (2018) Schaefer B. E., 2018, MNRAS, 481, 3033
  * Schaefer et al. (1992) Schaefer B. E. et al., 1992, ApJS, 81, 321
  * Schaefer et al. (2000) Schaefer B. E., King J. R., Deliyannis C. P., 2000, ApJ, 529, 1026
  * Schlafly and Finkbeiner (2011) Schlafly E.F., Finkbeiner D.P., 2011, ApJ 737, 103
  * Shara et al. (2018) Shara M. M., Prialnik D., Hillman Y., Kovetz A., 2018, ApJ, 860, 110
  * Shen and Bildsten (2009) Shen K. J., Bildsten L., 2009, ApJ, 699, 324
  * Shibata et al. (2013) Shibata K. et al., 2013, PASJ, 65, 49
  * Shibayama et al. (2013) Shibayama T., Maehara H., Notsu S., Notsu Y., 2013, ApJS, 209, 5
  * Simon et al. (1997) Simon T., Linsky J. L., Schiffer F. H. III, 1980, ApJ, 239, 911
  * Starrfield et al. (2008) Starrfield S., Iliadis C., Hix W., 2008, in Bode M. F., Evans A., eds, Classical Novae, second edition. Cambridge University Press, Cambridge, p. 77
  * Stella et al. (1988) Stella L., Haberl F., Lewin W. H. G., Parmar A. N., van Paradijs J., White N. E., 1988, ApJ, 324, 379
  * Tu et al. (2020) Tu Z.-L., Yang M, Zhang Z. J., Wang F. Y., 2020, ApJ, 890, 46
  * Vanderburg & Johnson (2014) Vanderburg A., Johnson J. A., 2014, PASP, 126, 948
  * Warner (2008) Warner B., 2008, in Bode M. F., Evans A., eds, Classical Novae, second edition. Cambridge University Press, Cambridge, p. 16
  * Yang and Liu (2019) Yang H., Liu J., 2019, ApJS, 241, 29


Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2017.DOI

This work was supported by the Russian Science Foundation (Project No.
20-61-47089)

Corresponding author: Lev Teplyakov (e-mail: teplyakov@iitp.ru).

# LSDNet: Trainable Modification of LSD Algorithm for Real-Time Line Segment
Detection

LEV TEPLYAKOV1  LEONID ERLYGIN1,2, EVGENY SHVETS1  Institute for Information
Transmission Problems, Russian Academy of Sciences, 119991 Moscow, Russia
Moscow Institute of Physics and Technology, 117303 Moscow, Russia

###### Abstract

As of today, the best accuracy in line segment detection (LSD) is achieved by
algorithms based on convolutional neural networks – CNNs. Unfortunately, these
methods utilize deep, heavy networks and are slower than traditional model-
based detectors. In this paper we build an accurate yet fast CNN-based
detector, LSDNet, by incorporating a lightweight CNN into a classical LSD
detector. Specifically, we replace the first step of the original LSD
algorithm – construction of line segments heatmap and tangent field from raw
image gradients – with a lightweight CNN, which is able to calculate more
complex and rich features. The second part of the LSD algorithm is used with
only minor modifications. Compared with several modern line segment detectors
on standard Wireframe dataset, the proposed LSDNet provides the highest speed
(among CNN-based detectors) of $214$ FPS with a competitive accuracy of $78$
$F^{H}$. Although the best-reported accuracy is $83$ $F^{H}$ at $33$ FPS, we
speculate that the observed accuracy gap is caused by errors in annotations
and the actual gap is significantly lower. We point out systematic
inconsistencies in the annotations of popular line detection benchmarks –
Wireframe and York Urban, carefully reannotate a subset of images and show
that (i) existing detectors have improved quality on updated annotations
without retraining, suggesting that new annotations correlate better with the
notion of correct line segment detection; (ii) the gap between accuracies of
our detector and others diminishes to negligible $0.2$ $F^{H}$, with our
method being the fastest.

###### Index Terms:

Convolutional neural networks, edge detection, line segment detection, U-net,
LSD

=-15pt

Figure 1: Overview of the proposed approach. A lightweight neural network
predicts line segment mask and tangent vector field, which are then clustered
and each cluster is represented as a line segment.

##  I Introduction

Automatic general-purpose line segment detection is a long-standing computer
vision problem of high practical importance. Line segment detectors are
exploited to construct an intermediate representation of image contents in
visual recognition systems over a wide range of applications, such as
autonomous vehicle localization [1, 2, 3], infrastructure maintenance with an
UAV [4, 5], document recognition [6, 7].

Traditionally, the problem of line segment detection was approached with so-
called model-based algorithms [8, 9, 10, 11]. These algorithms operate by
searching an image for elements that satisfy an explicit definition of a
salient line segment, for example, ``line segment is a strip-like set of image
pixels with similar gradients'' [8, 9] or ``an image region is a line segment
if its contour map triggers a peak in Hough space'' [10, 11]. These algorithms
typically have the benefits of being fast and having interpretable parameters.
However, they may miss the segments that are salient for human, but for some
reason don't match the exact explicit implemented definition. They are also
prone to over-segmentation (splitting a single segment into parts) and
sometimes demand nontrivial problem-specific postprocessing [12, 13].

The troubling problem of formulating an explicit criterion that matches the
human expectation of what exactly constitutes a ``salient line segment'' can
be avoided by manually annotating images and training a CNN, which then learns
an implicit algorithm from data samples. This is an approach that yields the
best accuracy in line segment detection task today [14, 15, 16, 17, 18, 19,
20]. Skipping ahead, let us note that such annotation is not a simple task
either - existing datasets on line segment detection have numerous and
sometimes extreme internal inconsistencies - probably caused by the inherent
ambiguity of the task, lack of clear labeling instructions and the tediousness
of the task, leading to missed segments.

From a technical perspective, there also is a challenge in designing a CNN-
based line detector. The typical solution is that CNN constructs an
intermediate representation - encoding - which is then converted into a set of
answers by some hand-crafted algorithm. Object detection networks have solved
this problem by using so-called anchors [21] and, later, more simple
anchorless detectors such as FCOS [22]. But line segment detectors need
alternative encodings, suitable not for bounding boxes, but for line segments.
The encodings should effectively deal with the fact that the line segments to
be detected in a typical image intersect with each other a lot, while
encodings for bounding boxes are effective only when ``overlapping mostly
happens between objects with considerably different sizes'' [22], making them
hard to exploit for line segment detection. While object detectors encondings
- after many iterations of refinement - have become fast and elegant, we
believe that intermediate representations of most line detector used today are
still either imprecise, slow or unintuitive.

So whether it is the complexity of the interpreter or the sheer weight of the
CNN backbone, CNN-based detectors that outperform the traditional ones in
accuracy are also computationally harder [20]. Their complexity limits the
scope of application of such algorithms in cases where speed, energy
consumption, or hardware price are critical.

In this work we propose a fast yet accurate CNN-based line segment detection
algorithm, LSDNet, built on the basis of a widely used model-based detector,
LSD [8]. LSDNet overview is presented in Figure 1. The first step of LSD is
the calculation of image gradient's orientation and magnitude. We view this
step as an estimation of an intermediate representation, composed of line
segments' heatmap, estimated as gradient's magnitude, and tangent field,
estimated as gradient's orientation. We substitute this step with a
lightweight CNN to generate the heatmap and tangent field from a more diverse
and complex set of features than a simple gradient; the second step,
conversion the intermediate representation into a set of line segments, is
taken from the LSD almost as-is. This substitution boosts the LSD accuracy and
simplifies the postprocessing due to more accurate heatmap and tangent field.
The heatmap and tangent field generation is a relatively simple task - the
answer can be correctly inferred from the local context - which allows to use
a lightweight CNN.

We compare LSDNet and a selection of existing detectors: LSD [8], L-CNN[18],
HAWP[15], TP-LSD[17], M-LSD and M-LSD-tiny[20] and show that LSDNet provides
competitive accuracy of $78$ $F^{H}$ score on Wireframe dataset while being
the fastest: $214$ FPS for $288\times 288$ input on conventional hardware.
LSDNet ourperforms the second fastest approach, M-LSD-tiny[20], both in
accuracy and speed.

We discuss considerable inconsistencies (section IV) in ground truth labeling
of current benchmarks for line segment detection accuracy - datasets Wireframe
[23] and YorkUrban [24]. The labeling in these datasets is inconsistent not
only between images: within the very same image many segments, similar in
appearance, are often marked up differently - some as a positives, others as
negatives. We speculate that these datasets in their current state are flawed
for assessing the accuracy of general purpose line segment detectors.

So, in addition to measuring the accuracy of our detector on these datasets,
we select and reannotate a part of Wireframe dataset - consisting mainly of
simple images with less ambiguity - and show that the reannotated subset
correlates better with conventional notion of correct line segment detection.
Specifically, it narrows the accuracy gap between LSD and CNN-based
approaches, and makes the gap between the proposed LSDNet and L-CNN [18], one
the most accurate CNN-based approaches, negligible.

##  II Related Work

In this section we cover the model-based LSD algorithm [8] and the existing
CNN-based approaches to the problem of line segment detection.

###  II-A LSD Detector

LSD is one of the most popular general-purpose line segment detectors and
serves as a common baseline for CNN-based algorithms both in accuracy and
speed [15, 16, 17].

The first step of LSD detector is gradient calculation; then the algorithm
builds so-called line support regions (LSRs) [25]. LSRs are image segments
(not to be confused with line segments), spanning actual line segments in an
image. They are built by iteratively grouping neighbouring pixels with high
gradients' magnitudes and similar orientations. After the formation of initial
LSRs, they undergo several steps of filtering and refinement. These include,
among others, splitting LSRs of ``hockey stick'' shape - effectively
decoupling distinct but merged regions; removal of LSRs with a high deviation
of gradients' orientations - possibly false-positives. Typical resulting LSRs
are long, straight and a few pixels thick. Finally, each such region is
individually encoded as a pair of points - effectively, a line segment.

LSD detector is fast, reaching 185 FPS for $320\times 320$ images on a
conventional CPU, and provides an accuracy of $63.3$ $F^{H}$ on the Wireframe
dataset.

###  II-B CNN-based approaches

The CNN-based approaches are typically composed of two modules: the CNN itself
predicts an intermediate representation, then a postprocessing module
reconstructs line segments from this representation.

We consider the design of the intermediate representation to be the key growth
point of CNN-based line segment detectors since the desired detector's output
- a unknown-size set of line segments with potentially high overlap - is hard
to represent as a ``CNN-friendly'' fixed-shape tensor [26]. The survey below
covers most popular intermediate representations used in existing CNN-based
approaches.

The first CNN-based detectors represented line segments in an image as a set
of endpoints and their connectivity graph [23, 18]. The endpoints were
detected as local maxima of a CNN-produced heatmap. The connectivity of the
endpoints was deduced either with the help of edge map heuristics [23] or with
a trainable classifier [18]. To generate the classifier's input, a fixed
number of uniform spaced points was sampled from the feature map between the
endpoints. This operation was called LoI pooling [18].

In [14] a representation called attraction field was proposed (distance field
[19] is a similar concept). Line segments were represented as a 2D vector
field of translations to a nearest line segment. The postprocessing step for
such a representation required nontrivial line segments extraction from
heatmap-like prediction. Interestingly, the postprocessing of even a perfect
attraction field generated from dataset annotations did not provide absolute
detection accuracy [14] \- in other words, such representation is inherently
ambiguous.

The covered representations of endpoints and attraction field were combined
and modified in [15]. The proposed CNN predicted both the endpoints and the
attraction field, which was enriched to encode the translations to the both
segment's endpoints as 4D vector field. Then the endpoints and the attraction
field were used to refine each other. The refined line segments proposals
underwent final verification with the help of LoI pooling and a trainable
classifier. The detector provides state-of-the-art quality of $83$ $F^{H}$ on
Wireframe [23] dataset to date, but with low speed of $33$ FPS on GPU.

The recently proposed ``tri-points'' representation [17] is focused on
speeding up the detector. A line segment was represented by its center point
and two vectors to its endpoints. It allowed to significantly boost the speed
up to $50$ FPS, driven by a much faster postprocessing requiring trivial
conversion from ``tri-point'' to a line segment and non-maximum supression.
The CNN itself remained comparably slow. In [20] some further enhancements
were proposed. A lightweight CNN was designed and the training procedure was
improved by augmentations and more sophisticated loss function. It resulted in
the fastest CNN-based detector to date with $200$ FPS overall and $241$ FPS
for standalone CNN.

##  III Proposed Approach

A perfect line segments representation should make it possible to design a CNN
and a postprocessing module both being fast and accurate. We argue that in the
search for such a representation there is no need to develop a brand new one
from scratch; instead, the representation used implicitly by LSD detector -
line segments heatmap and tangent field - already possesses all the desired
properties. Indeed, the heatmap and tangent field could be inferred from local
image context and does not require the reasoning of complex abstract features,
which allows to use a lightweight CNN. On the other hand, as proven by LSD,
the representation could be efficiently postprocessed to actual line segments.
In the next sections we cover each step of the algorithm in detail.

###  III-A CNN

Line Segment Representation. We represent a set of line segments in an image
as a 2-channel feature map of the same height and width as the image. The
first channel denoted by $M$ contains a line segment mask. The second channel
denoted as $F$ contains tangent vector field of line segments. Since the
values of $F$ are unit vectors, we encode $F$ as one channel feature map.

Let $l=(x_{1},y_{1},x_{2},y_{2})$ be a line segment, $\varphi_{l}$ \-
segment's level line angle in range $[0,\pi)$, $p=(x,y)$ \- an image pixel,
$L_{p}=l_{1},l_{2},...$ \- set of line segments crossing $p$. Then if
$L_{p}=\emptyset$, then $M(p)=0$ and $F(p)$ is arbitrary; otherwise $M(p)=1$,
$F(p)=\varphi_{l_{1}}$.

Note that in case of overlapping segments ($||L_{p}||>1$) $F(p)$ is defined by
one arbitrary segment of overlap. This ambiguity probably can be ignored since
only about $1\%$ of line segments' pixels lie on overlaps - as measured on
Wireframe dataset [23].

Loss Function. The loss function $L=L_{mask}+\alpha L_{field}$ used to fit the
network is composed of two independent weighted terms, one responsible for
mask $M$, other - for the tangent vector field $F$.

While the prediction of mask $M$ is a straightforward segmentation problem
with $L_{mask}$ being a conventional cross enthropy loss, to correctly
estimate error in prediction of the tangent field $F$ we should account for
the following property of line segments' level line angles - the distance
between angles $0^{0}$, $10^{0}$ and $0^{0}$, $170^{0}$ should be equal.

This problem can be approached by computing several distances between the
original angles and the angles, shifted by $\pm\pi$, and picking the minimum
distance [27]. The calculation can be done simpler: let $\varphi_{1}$,
$\varphi_{2}$ \- angles between which the distance is to be computed,
$z_{1}=e^{i\varphi_{1}},z_{2}=e^{i\varphi_{2}}\in\mathbf{C}$ \- the
representation of the angles as complex numbers with unit length and phases
$\varphi_{1}$, $\varphi_{2}$, then

$\rho(\varphi_{1},\varphi_{2})=||z_{1}^{2}-z_{2}^{2}||_{2}$ (1)

This distance function has a geometrical interpretation, illustrated in Fig.
2a. It equals to $\ell^{2}$-norm of vector difference between unit vectors
with phases $2\varphi_{1}$, $2\varphi_{2}$. Turning to the aforementioned
example, the doubled phase makes vectors, corresponding to angles $10^{0}$ and
$170^{0}$, equally close to the horizontal, corresponding to angle $0^{0}$.

Given the angle distance function $\rho$, the predicted field $F_{p}$, the
reference mask $M_{t}$ and reference field $F_{t}$, the tangent vector field
loss $L_{field}$ is defined as

$L_{field}=\frac{1}{\sum_{p}M_{t}(p)}\sum_{p:M_{t}(p)=1}\rho^{2}(F_{t}(p),F_{p}(p)),$
(2)

where $M_{t}$ is the reference line segment mask - essentially, loss is the
average tangent angle discrepancy over the pixels that correspond to the
ground truth line segments.

(a)

(b)

Figure 2: (a) A visualization of the distance function between angles
$\varphi_{1}$, $\varphi_{2}$. It accounts for the periodicity of angles. (b)
Line segment extraction from a region.

CNN Architecture. To predict the proposed feature map, we use a CNN of U-Net
[28] family. The architecture we use differs from the original one in the
following simplifications. We exploit padded convolutions providing the same
input and output spatial sizes of convolutional layers, which allows not to
crop the feature maps feeded to skip connections. Instead of transposed
convolutions we use bilinear upsampling. We reduce the depth of encoder-
decoder branches up to 3 maxpooling and 3 upsampling layers, correspondingly,
and use fewer filters in convolutional blocks - 16, 32, 64, 128 filters per
block (the number of blocks is greater than the number if maxpooling layers by
one). The resulting CNN has $\approx 0.5$M trainable parameters and can run at
48 FPS on CPU and at 695 FPS on GPU (refer to section V for benchmarking
details).

###  III-B Line Segments Reconstruction

Let us consider how the predicted segments mask $M$ and tangent field $F$ are
converted into the desired output - a set of line segment
$((x_{1},y_{1}),(x_{2},y_{2})),...$. This process has three steps: firstly,
the predicted features are coarsely segmented into lines and background
(section ``Foreground segmentation''). Secondly, the lines are finely
segmented into several line support regions (LSRs) (section ``Region
grouping''). Finally, a line segment and its confidence is extracted from each
LSR (section ``Line segments extraction'').

Foreground segmentation. The first step of line segments reconstruction is to
segment foreground (lines) from background (not lines).

We use a coarse-to-fine binarization approach by multiplying the masks of
global thresholding $M(p)>\tau$ and local thresholding

$M(p)>\sum_{d\in K(p)}W_{p}(d)M(d)-\theta$ (3)

where $\theta$ \- threshold, $K_{p}$ \- a window centered at pixel $p$,
$W_{p}(d)$ \- Gaussian averaging weight.

Global thresholding with a small threshold gives a coarse extraction of line
segments mask, but often incorrectly joins close - but separate - line
segments. On the contrary, local thresholding provides much finer local
distinction of line segments, but can produce clumps of false positive
detections in low intensity areas (Fig. 3). The combination of these
binarizations by simple multiplication of the resulting masks allows to filter
out false positives of both types and achieve better accuracy.

(a) Prediction

(b) Global

(c) Local

(d) Combined

Figure 3: Binarization of CNN prediction. Global binarization joins close
regions, while local binarization is noisy for the regions of near-zero
intensity. The combined binarization is free of both drawbacks.

Region grouping. The goal of this step (being a modification of a similar step
of LSD algorithm [8]) is to split the foreground, binarized at the previous
step, into narrow strip-like LSRs, one per true line segment.

Informally, we want neighbouring pixels $p_{1},p_{2}$ to be assigned to one
LSR, if the values of $M(p_{1})$, $F(p_{1})$ and $M(p_{2})$, $F(p_{2})$ are
similar. The algorithm grows LSRs iteratively, starting from pixels with
highest $M$ value and adding new pixels to the existing LSR which are
geometrically close to the pixels and have similar features.

Let us formally introduce the similarity measure used to decide whether pixel
$g$ is fit to be joined into a LSR. Let $R=p_{1},p_{2},...,p_{n}$ be the set
of pixels of this LSR, $I_{R}=1/n\cdot\sum_{p\in R}M(p)$ be the mean line
segments' mask over it, and $\phi_{R}=\angle\sqrt{\sum_{p\in R}e^{2iF(p)}}$ \-
the average tangent field (here $\angle
z=\operatorname{atan2}(\operatorname{Im}(z)/\operatorname{Re}(z))$ is the
phase of a complex number). Then the similarity function is given by

$d_{M,F}(g,R)=\rho^{2}(F(g),\phi_{R})+\alpha(M(g)-I_{R})^{2}$ (4)

The first term defines similarity of tangent field orientation (refer to Eq.
(1) for details), the second - the similarity of line mask, $\alpha$ \-
weighting coefficient. Given the distance function, LSRs are built with an
iterative growing algorithm, presented in Algorithm 1.

input : Line heatmap $M\in[0,1]^{hw}$,
tangent field $F\in[0,\pi)^{hw}$

output : Set of line support regions $\mathcal{R}=\\{R\\}$

param : $\tau\in\mathbf{R}$ \- distance threshold

param : $n\in\mathbf{N}$ \- minimum region size

$B\in\\{0,1\\}^{hw}$ \- foreground segmentation of $M$;

$P=(p_{i})$, such that $B(p)=1,M(p_{i})\geq M(p_{i+1})$;

$P_{used}\leftarrow\varnothing$;

for  _$p\in P$_ do

if  _$p\in P_{used}$_ then

continue;

end if

$R\leftarrow p$;

do

$region\\_updated\leftarrow False$;

for  _$g\in N_{8}(R)\setminus P_{used}$_ do

if  _$d_{M,F}(g,R) <\tau$_ then

Add $R\leftarrow g$;

Add $P_{used}\leftarrow g$;

$region\\_updated\leftarrow True$;

end if

end for

while  _$region\\_updated$_ ;

Add $\mathcal{R}\leftarrow R$;

end for

for  _$R\in\mathcal{R}$_ do

if  _$|R| <n$_ then

remove $R$ from $\mathcal{R}$;

end if

end for

Algorithm 1 Region grouping algorithm. $d_{M,F}(p,R)$ \- distance function
between pixel $p$ and region $R$, $N_{8}(R)$ \- 8-neighbours of region $R$.

Line segments extraction. Each LSR $R=p_{1},p_{2},...,p_{n}$, should be
converted into a line segment satisfying the following criteria.

  * •

The segment goes through LSR's center of mass $p_{\mu}$

$p_{\mu}=\frac{1}{\sum_{p\in R}M(p)}\sum_{p\in R}M(p)p.$ (5)

  * •

The segment is collinear with minor eigenvector $a$ of region's inertia tensor
$I$ defines as follows

$I=\sum_{p\in R}M(p)I^{*}(p-p_{\mu})$ (6)
$I^{*}(x,y)=\begin{pmatrix}y^{2}&-xy\\\ -xy&x^{2}\\\ \end{pmatrix}$ (7)

  * •

The segment spans the furthest LSR's points, projected onto axis $a$. The
segment's confidence is mean value of $M$ over the region $R$.

Line segment extraction is visualized in Figure 2b.

##  IV Dataset

In this section we analyze the issues of the existing line segment detection
datasets and propose a dataset Wireframe-tiny++, a subset of Wireframe dataset
[23] with refined annotations.

###  IV-A The existing datasets

To the best of our knowledge, there are two widely-used public line segment
detection datasets: Wireframe [23] and YorkUrban [24]. The former is composed
of 5.000 train and 462 test images, the latter is composed of 120 test images.
The datasets contain both indoor and outdoor colour images of various man-made
environments. Some samples from the datasets are presented in Figures 4 and 5.
The datasets are annotated with a list of point pairs, representing line
segments.

York dataset was annotated under so-called Manhattan world assumption [29],
which means that the annotated line segments are those aligned with the basis
of some Cartesian coordinate system (specifically, with axes parallel to image
sides), while the others are ignored. Wireframe dataset did not follow
Manhattan world assumption and was annotated with line segments, from which
``meaningful geometric information of the scene can be extracted'' [23], which
also resulted in some salient line segments being not annotated.

So, the ground truth labeling in these datasets is explicitly limited to some
category of line segments - which means other categories of line segments are
viewed as negatives. CNNs that are trained on these datasets (and/or with high
accuracy on them) will have to systematically classify these categories of
line segments as negatives and therefore - by design! - can't be viewed as
general-purpose line detectors.

While such annotations can be useful to train and test some specific niche
line segment detectors (e.g. for indoor robot navigation [30]) we believe they
are flawed as datasets for general purpose line segment detection.
Specifically, we would like to highlight the following problems (also
illustrated in Fig. 4b).

(a) Original image

(b) Wireframe annotation

(c) Wireframe-tiny++ annotation

Figure 4: An image from Wireframe dataset (a) and two versions of its
annotation - original (b) and proposed (c). Best viewed in color and zoom-in.

One problem is the inconsistency of the annotations - it is easily noticeable
on strip-like objects having two side line segments looking almost exactly
alike - however one of them is annotated while the other is not.

Some categories of salient line segments are systematically not annotated -
e.g. shadows and reflections. We believe the fact that these segments are not
``real'' physical objects should not be considered in the context of general
purpose line segment detection and these segments should be annotated as well.

Finally, some line segments lying on the same straight line are falsely merged
(vertical segments on the bed canopy's frame in Fig. 4b). It happens when a
long line segment is intercepted by another object. Although for some
applications it could be desirable to avoid such a splitting and there are
approaches to achieve that [12, 13], we believe, that for general-purpose
detector splitting is the desired detector's behaviour.

###  IV-B Wireframe-tiny++

To approach the covered issues with the existing datasets, we selected 20
random images from Wireframe test subset and reannotated them to make the
annotations more accurate and consistent. We call the selected subset of
images with the original markup Wireframe-tiny, and the resulting dataset with
enhanced annotations - Wireframe-tiny++.

Comparing to the original annotations, we mainly added unannotated segments, 9
per image on average. Some segments are removed as undetectable. Some segments
are divided into several smaller segments due to occlusion. The refined
annotations are presented in Fig. 4c.

##  V Experimental Setting

Datasets. The proposed algorithm is trained and evaluated with the following
datasets. Wireframe dataset [23] consisting of 5000 training and 462 test
images is used both to train and evaluate LSDNet. Datasets YorkUrban [24],
Wireframe-tiny and its reannotated version Wireframe-tiny++ (refer to the
previous section for details), composed of 120 and 20 images correspondingly,
are used solely for evaluation.

Accuracy. To evaluate LSDNet accuracy we use standard [15, 18] quality score
$F^{H}=100\cdot 2\cdot p\cdot r/(p+r)$, where $p,r$ stand for precision and
recall. Multiplier $100$ is added for readability, making $F^{H}$ fall in
$[0,100]$ range. The score is evaluated pixel-wise by rasterizing both the
predicted and the reference line segments. A pixel of a predicted line segment
is considered true positive, if its distance to a pixel of a reference segment
does not exceed $1\%$ of image diagonal. For evaluation we use $F^{H}$
implementation provided with L-CNN [18].

Quality score $F^{H}$ was criticised [18] for being not sensitive towards
overlapping and splitted line segments. We consider such an insensibility is
not critical: LSDNet can not produce overlapping segments by design, since
line support regions (LSRs) can not overlap, and we did not observe a notable
amount of splitted line segments for any CNN-based algorithm.

Speed. CNN is benchmarked on Quadro GV100 GPU for comparison with other
detectors. Reconstruction algorithm is benchmarked on Core i5 9300hf CPU. CNN
and the reconstruction algorithm are benchmarked independently. The speed of
the latter depends on its input, we report the average speed over the dataset
given the trained preprocessing network.

Baselines. We compare the proposed LSDNet with classical algorithm LSD [8] and
several state-of-the-art CNN-based detectors L-CNN[18], HAWP[15], TP-LSD[17],
M-LSD and M-LSD-tiny[20].

The reported FPS for all CNN-based methods is cited as in [20], where
benchmarking was performed on Tesla V100 GPU with practically the same
characteristics as the GPU used in our experiments.

The reported $F^{H}$ is also cited as in [20] for all methods except LSD and
L-CNN [18], for which it was reproduced by our means. We tried to reproduce
the stated quality measurements for other approaches with the help of their
open-source implementations, but they appeared notably lower than the reported
ones. Therefore on datasets Wireframe-tiny and Wireframe-tiny++ we compare
LSDNet only to L-CNN and LSD.

Preprocessing. For LSDNet, all images are resized to $288\times 288$, which
appeared to be the optimal input shape in terms of speed-accuracy tradeoff.
Pixel intensities are simply converted from 8-bit unsigned integer to 32-bit
floating-point with $1/255$ scaling coefficient. During training, random
horizontal and vertical flips and gamma correction are applied.

For baseline methods, the preprocessing from the corresponding paper is
applied. For LSD, we use $320\times 320$ image shape.

Hyperparameters. LSDNet is initialized with He uniform initialization[31] and
trained by Adam optimizer [32] with $10^{-4}$ weight decay and $8$ images per
batch for $180$ epochs. The initial learning rate is $10^{-3}$ and is reduced
by half if the value of loss function does not improve for $15$ epochs.

Implementation. CNN training and inference is implemented in TensorFlow [33]
and ONNX Runtime [34], correspondingly. Reconstruction algorithm is
implemented in C++ with the help of OpenCV [35].

##  VI Results and Analysis

In this section we quantitatively and qualitatively analyze LSDNet performance
and compare it to a wide range of state-of-the-art line segment detectors.
Please refer to Section V for evaluation and comparison details.

Figure 5: Speed (FPS) and accuracy ($F^{H}$) (accuracy) comparisons of different line segment detectors. Size of circles indicates the number of trainable parameters. LSD is marked as cross since it has no trainable parameters. The proposed LSDNet outperforms all previous detectors in speed with $214$ FPS and outperforms the nearest fastest counterpart, M-LSD-tiny, both in speed and $F^{H}$ quality. Method | accuracy, $F^{H}$ | speed, FPS
---|---|---
Wireframe | York | CNN | postproc | total
LSD [8] | 63.3 | 59.0 | - | 185 | 185
L-CNN [18] | 81.2 | 65.4 | 55 | 24 | 17
HAWP [15] | 83.1 | 66.3 | 55 | 82 | 33
TP-LSD [17] | 82.0 | 67.3 | 65 | 201 | 49
M-LSD [20] | 80.0 | 64.2 | 132 | 883 | 115
M-LSD-tiny [20] | 76.8 | 61.9 | 241 | 1203 | 200
LSDNet (proposed) | 77.5 | 64.6 | 695 | 301 | 214
TABLE I: Detection accuracy on datasets Wireframe and York Urban. Bold and underlined values stand for top-1 and top-2, correspondingly. Column ``postproc'' shows the speed of postprocessing the CNN's prediction. Method | detection accuracy, $F^{H}$
---|---
Wireframe | Wireframe-tiny | Wireframe-tiny++
LSD[8] | 63.3 | 69.0 | 73.2
L-CNN [18] | 81.2 | 83.6 | 85.5
LSDNet (proposed) | 77.5 | 82.0 | 85.3
TABLE II: Detection accuracy on datasets Wireframe, Wireframe-tiny, Wireframe-
tiny++ datasets. Bold and underlined values stand for top-1 and top-2,
correspondingly.

Figure 6: Qualitative evaluation. From left to right - images from Wireframe
dataset of varying number of salient line segments. From top to bottom -
original image, detection results from LSD [8], the model-based detector; HAWP
[15], the most accurate CNN-based detector; M-LSD-tiny [20], the fastest CNN-
based detector except for the proposed one; the proposed LSDNet. Best viewed
in color and zoom-in.

Public datasets. Table I and Figure 5 summarize the results on Wireframe and
York Urban datasets. It shows that LSDNet achieves state-of-the-art inference
speed of $695$ FPS for standalone CNN and $214$ FPS for overall detector
alongside competitive accuracy with $77.5$ $F^{H}$ and $64.6$ $F^{H}$ on
Wireframe and York Urban datasets, correspondingly.

In comparison with the nearest fastest counterpart, M-LSD-tiny, LSDNet
outperforms it both in quality and speed with the absolute increase in
accuracy of $+1.0$ $F^{H}$ on Wireframe, $+2.7$ $F^{H}$ on York Urban and
$+14$ FPS speed-up.

Compared to the fastest detector outperformimg LSDNet in accuracy on Wireframe
dataset, M-LSD, the proposed approach is approximately two times faster with
$214$ FPS against $115$ FPS. The fastest algorithm to outperform LSDNet on
both datasets is TP-LSD, which is approximately four times slower with $49$
FPS.

Custom datasets. Table II summarizes the results on Wireframe dataset, its
subset Wireframe-tiny and its reannotated version Wireframe-tiny++. Please
refer to section IV-B for details.

On Wireframe-tiny++, all the detectors demonstrate higher accuracy than those
on Wireframe-tiny. Since these datasets are composed of the same images and
differ only in annotations, such a consistent accuracy growth indicates that
Wireframe-tiny++ annotation is more suitable for the problem of general
purpose line segment detection.

All the approaches, being arranged by $F^{H}$, show the same relative order on
all the datasets, but the absolute differences change significantly. The gap
between LSDNet and L-CNN has shrinked from $3.7$ $F^{H}$ on Wireframe to
negligible $0.2$ $F^{H}$ on Wireframe-tiny++. We believe it could be explained
by different learning capacity of detectors' CNNs. Expressive L-CNN with 9.8M
parameters managed to learn the subtle notion of a line segment implied by
Wireframe train dataset annotation (discussed in Sec. IV-A); whereas
lightweight LSDNet with only 0.5M parameters learned the general line segment
detection with no capacity to learn the subtle details. It made L-CNN good for
wireframe-like detection problems with the goal to detect line segments, from
which ``meaningful geometric information of the scene can be extracted''. But
it could possess confusing properties in terms of general-purpose line segment
detection, making LSDNet a better choice in such a case.

Qualitative results. Qualitative comparison of LSDNet to other line segment
detectors is illustrated in Figure 6. In this section we refer to LSD and
LSDNet as LSR-based and to HAWP and M-LSD-tiny as endpoint-based detectors,
since the methods within these groups demonstrate similar behaviour.

Endpoint-based detectors demonstrate the selectivity of line segments, which
could not be attributed to overall segments' saliency. This effect is mostly
notable in the foreground in the right column in Figure 6. LSR-based LSD
detects all the shadows and the carpet in front of the sofa, while it can't
``see'' floor tiles due to their low contrast. LSDNet detects all the shadows,
the carpet and the floor tiles. Whereas endpoint-based detectors HAWP and
M-LSD-tiny detect these objects poorly, but at the same time they detect way
less salient segments in the background. We believe such a selectivity could
be attributed to the combination of high expressive power of the underlying
CNN and annotation inconsistencies of train dataset, discussed in Section
IV-A. This effect could be undesirable in an application requiring that very
class of segments, which is missed by endpoint-based detectors.

Another interesting difference between the detectors' groups occurs due to the
types of misdetections. In terms of a quality measure, LSR-based LSD and
LSDNet can detect false positives, typically corresponding to line segment-
like patterns on highly structured image regions, and miss some annotated
segments (false negatives), which are usually poorly visible. These errors
could be, at least partially, attributed to the ill-posed nature of the task.
The endpoints-based methods are also prone to miss poorly visible segments,
and possess an advantage of not detecting false positives on structured image
regions. However, an potential drawback of endpoints-based detectors is that
they can produce hard false positives — line segments of high confidence score
with no evidence of a true line segment in an image. We believe the reason for
hard false positives is a classification error of line segment verification
module. An example could be seen in the middle column in Figure 6, please note
the salient diagonal segments in the bedhead (fourth row) and right part of
the carpet (third row). This issue is to be approached prior to successful
exploitation of an endpoint-based detector.

##  VII Conclusion

In this study we introduce a fast and accurate line segment detector LSDNet.
The detector is composed of a lightweight encoder-decoder CNN, which predicts
line segment heatmap and tangent field, and a postprocessing module – a
modification of the famous LSD algorithm. When benchmarked on the traditional
Wireframe dataset against several SOTA methods, LSDNet shows the highest FPS
of $214$ – though it achieves detection accuracy of $78$ $F^{H}$ – lower than
the best methods (82 and 83.1 $F^{H}$).

However, we speculate that this gap in detection accuracy is primarily caused
by the imperfections of the dataset rather than the network itself. We analyze
the commonly used line segment detection datasets – Wireframe and York Urban –
and point out numerous and significant inconsistencies in their annotation. By
carefully reannotating a part of the Wireframe test dataset, we show that (i)
all detectors demonstrate better quality on improved annotations (without any
re-training), which indicates that the refined annotations correlate better
with the notion of correct line segment detection, (ii) the gap between
accuracies of our detector and others is reduced to almost non-existent - with
our method being the fastest.

##  VIII Acknowledgements

We would like to thank Marina Tepliakova for making the illustrations; Alexey
Savchik and Veniamin Blinov for reviewing the early versions of the
manuscript; Dmitry Nikolaev for his strong belief in the power of fusing
model-based algorithms with light-weight neural networks, which inspired this
work.

## References

  * [1] D. C. Andrade, F. Bueno, F. R. Franco, R. A. Silva, J. H. Z. Neme, E. Margraf, W. T. Omoto, F. A. Farinelli, A. M. Tusset, S. Okida _et al._ , ``A novel strategy for road lane detection and tracking based on a vehicle’s forward monocular camera,'' _IEEE Transactions on Intelligent Transportation Systems_ , vol. 20, no. 4, pp. 1497–1507, 2018.
  * [2] T. M. Hoang, H. G. Hong, H. Vokhidov, and K. R. Park, ``Road lane detection by discriminating dashed and solid road lanes using a visible light camera sensor,'' _Sensors_ , vol. 16, no. 8, p. 1313, 2016.
  * [3] O. Shipitko, V. Kibalov, and M. Abramov, ``Linear features observation model for autonomous vehicle localization,'' 2020.
  * [4] A. Ceron, F. Prieto _et al._ , ``Power line detection using a circle based search with uav images,'' in _2014 International Conference on Unmanned Aircraft Systems (ICUAS)_. IEEE, 2014, pp. 632–639.
  * [5] Y. Cao and L. Yan, ``Automatic road network extraction from uav image in mountain area,'' in _2012 5th International Congress on Image and Signal Processing_. IEEE, 2012, pp. 1024–1028.
  * [6] J. Shemiakina, A. Zhukovsky, and D. Nikolaev, ``The method for homography estimation between two planes based on lines and points,'' in _Tenth International Conference on Machine Vision (ICMV 2017)_ , vol. 10696. International Society for Optics and Photonics, 2018, p. 106961G.
  * [7] A. Zhukovsky, D. Nikolaev, V. Arlazarov, V. Postnikov, D. Polevoy, N. Skoryukina, T. Chernov, J. Shemiakina, A. Mukovozov, I. Konovalenko _et al._ , ``Segments graph-based approach for document capture in a smartphone video stream,'' in _2017 14th IAPR International Conference on Document Analysis and Recognition (ICDAR)_ , vol. 1. IEEE, 2017, pp. 337–342.
  * [8] R. G. Von Gioi, J. Jakubowicz, J.-M. Morel, and G. Randall, ``Lsd: A fast line segment detector with a false detection control,'' _IEEE transactions on pattern analysis and machine intelligence_ , vol. 32, no. 4, pp. 722–732, 2008\.
  * [9] C. Akinlar and C. Topal, ``Edlines: Real-time line segment detection by edge drawing (ed),'' in _2011 18th IEEE International Conference on Image Processing_. IEEE, 2011, pp. 2837–2840.
  * [10] R. F. Guerreiro and P. M. Aguiar, ``Connectivity-enforcing hough transform for the robust extraction of line segments,'' _IEEE Transactions on Image Processing_ , vol. 21, no. 12, pp. 4819–4829, 2012.
  * [11] Z. Xu, B.-S. Shin, and R. Klette, ``Accurate and robust line segment extraction using minimum entropy with hough transform,'' _IEEE Transactions on Image Processing_ , vol. 24, no. 3, pp. 813–822, 2014.
  * [12] N. Hamid and N. Khan, ``Lsm: perceptually accurate line segment merging,'' _Journal of Electronic Imaging_ , vol. 25, no. 6, p. 061620, 2016.
  * [13] Q. Yu, G. Xu, Y. Cheng, and Z. H. Zhu, ``Plsd: A perceptually accurate line segment detection approach,'' _IEEE Access_ , vol. 8, pp. 42 595–42 607, 2020.
  * [14] N. Xue, S. Bai, F. Wang, G.-S. Xia, T. Wu, and L. Zhang, ``Learning attraction field representation for robust line segment detection,'' in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1595–1603.
  * [15] N. Xue, T. Wu, S. Bai, F. Wang, G.-S. Xia, L. Zhang, and P. H. Torr, ``Holistically-attracted wireframe parsing,'' in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 2788–2797.
  * [16] Y. Lin, S. L. Pintea, and J. C. van Gemert, ``Deep hough-transform line priors,'' _arXiv preprint arXiv:2007.09493_ , 2020.
  * [17] S. Huang, F. Qin, P. Xiong, N. Ding, Y. He, and X. Liu, ``Tp-lsd: Tri-points based line segment detector,'' _arXiv preprint arXiv:2009.05505_ , 2020.
  * [18] Y. Zhou, H. Qi, and Y. Ma, ``End-to-end wireframe parsing,'' in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 962–971.
  * [19] K. Huang and S. Gao, ``Wireframe parsing with guidance of distance map,'' _IEEE Access_ , vol. 7, pp. 141 036–141 044, 2019.
  * [20] G. Gu, B. Ko, S. Go, S.-H. Lee, J. Lee, and M. Shin, ``Towards real-time and light-weight line segment detection,'' _arXiv preprint arXiv:2106.00186_ , 2021.
  * [21] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao, ``Yolov4: Optimal speed and accuracy of object detection,'' _arXiv preprint arXiv:2004.10934_ , 2020.
  * [22] Z. Tian, C. Shen, H. Chen, and T. He, ``Fcos: Fully convolutional one-stage object detection,'' in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 9627–9636.
  * [23] K. Huang, Y. Wang, Z. Zhou, T. Ding, S. Gao, and Y. Ma, ``Learning to parse wireframes in images of man-made environments,'' in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 626–635.
  * [24] P. Denis, J. H. Elder, and F. J. Estrada, ``Efficient edge-based methods for estimating manhattan frames in urban imagery,'' in _European conference on computer vision_. Springer, Berlin, Heidelberg, 2008, pp. 197–210.
  * [25] J. B. Burns, A. R. Hanson, and E. M. Riseman, ``Extracting straight lines,'' _IEEE transactions on pattern analysis and machine intelligence_ , no. 4, pp. 425–455, 1986.
  * [26] L. Teplyakov, K. Kaymakov, E. Shvets, and D. Nikolaev, ``Line detection via a lightweight cnn with a hough layer,'' _arXiv preprint arXiv:2008.08884_ , 2020\.
  * [27] Z. Li, F. Wang, and N. Wang, ``Lidar r-cnn: An efficient and universal 3d object detector,'' _arXiv preprint arXiv:2103.15297_ , 2021.
  * [28] O. Ronneberger, P. Fischer, and T. Brox, ``U-net: Convolutional networks for biomedical image segmentation,'' in _International Conference on Medical image computing and computer-assisted intervention_. Springer, Cham, 2015, pp. 234–241.
  * [29] J. M. Coughlan and A. L. Yuille, ``Manhattan world: Orientation and outlier detection by bayesian inference,'' _Neural computation_ , vol. 15, no. 5, pp. 1063–1088, 2003.
  * [30] E. Delage, H. Lee, and A. Y. Ng, ``Automatic single-image 3d reconstructions of indoor manhattan world scenes,'' in _Robotics Research_. Springer, Berlin, Heidelberg, 2007, pp. 305–321.
  * [31] K. He, X. Zhang, S. Ren, and J. Sun, ``Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,'' in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1026–1034.
  * [32] D. P. Kingma and J. Ba, ``Adam: A method for stochastic optimization,'' _arXiv preprint arXiv:1412.6980_ , 2014.
  * [33] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , ``Tensorflow: A system for large-scale machine learning,'' in _12th symposium on operating systems design and implementation_ , 2016, pp. 265–283.
  * [34] ``Optimize and accelerate machine learning inferencing and training.'' [Online]. Available: https://www.onnxruntime.ai/
  * [35] G. Bradski, ``Opencv library,'' Accessed on: Apr. 22, 2022. [Online]. Available: https://opencv.org/

|  Lev Teplyakov received the B.S. and M.S. degrees at Moscow Institute for
Physics and Technology, Moscow, Russia in 2017 and 2019, correspondingly. We
is currently pursuing the Ph.D. degree at Institute for Information
Transmission Problems, Moscow, Russia. His master`s thesis considered the
problem of training cascade CNN classifiers, which process an image with fast
or slow branches, depending on the complexity of an image. His research
interests include video analytics, action recognition, deep learning for
object detection, fast neural networks.
---|---
|  Leonyd Erlygin graduated from Moscow Institute for Physics and Technology
and received his bachelor degree in 2021. Now he is a Master`s student at
Institute for Information Transmission Problems, Moscow, Russia. His research
interests include CNN-based object detection and classification, with the
focus on generation of synthetic training datasets with the help of data
augmentation and 3d modelling.
---|---
|  Evgeny Shvets graduated from Moscow Institute for Physics and Technology,
Moscow, Russia and received his Ph.D. degree in Technology in 2017 in
Institute for Information Transmission Problems, Moscow, Russia. His Ph.D.
thesis focused on distributed control of a swarm for area surveillance. His
research interests include image processing, image registration and deep
learning, including generation of synthetic datasets for neural network
training.
---|---


# Efficient Calculations of the Mode-Resolved ab-initio Thermal Conductivity
in Nanostructures

Giuseppe Romano<EMAIL_ADDRESS>Institute for Soldier Nanotechnologies,
Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge,
Massachusetts 02139, USA

###### Abstract

First-principles calculations of thermal transport in homogeneous materials
have reached remarkable predicting power. Modeling deterministically phonon
transport in nanostructures, however, poses novel challenges; notably, it
entails solving as many algebraic equations as the number of combinations of
wave vectors in the discretized Brillouen and polarizations. We show that,
within the relaxation time approximation of the Boltzmann transport equation
(BTE), this issue is resolved by interpolating the phonon distributions in the
vectorial phonon mean free paths (MFP) space. The coupling between structure
and mode-resolved heat transport is investigated in terms of angular-resolved
bulk thermal conductivity and phonon suppression function, the latter being
associated primarily to the material’s geometry. Our method, termed the
anisotropic MFP-BTE (aMFP-BTE), allows for fast and accurate thermal
conductivity calculations in nanomaterials regardless of the number of phonon
branches and wave vectors. Furthermore, it naturally blends with first-
principles thermal transport calculations, therefore allowing for multiscale,
parameter-free simulations. We apply the aMFP-BTE to compute the mode-resolved
effective thermal conductivity of porous Si membranes, achieving up to 50x
speed with respect to the case with no interpolation. The proposed approach
unlocks the engineering of novel nanostructures, with applications to
thermoelectrics and heat management.

##  I Introduction

Tuning thermal transport via nanostructuring is beneficial to several
applications, including thermal energy harvesting Vineis  _et al._ (2010);
Lorenzi  _et al._ (2018); Lee  _et al._ (2010), heat management Kim  _et al._
(2007); Moore and Shi (2014) and routing Anufriev  _et al._ (2017); Zhang  _et
al._ (2013). The key advantage of nanomaterials is the possibility of
exploiting ballistic phonon effects, in contrast to macroscopic systems, where
transport is mostly diffusive Chen (2005); Murthy  _et al._ (2005). Modeling
this transport regime, however, is challenging because it entails solving, at
the very least, the Boltzmann Transport Equation (BTE), a much more
complicated model than the standard Fourier equation Ziman (2001). The steady-
state, mode-resolved BTE, in fact, requires tracking phonons in both real and
momentum space, and, when solved deterministically, amounts to inverting as
many matrices as the number of phonon modes in our system Murthy and Mathur
(1998), which can easily reach hundreds of thousands in realistic materials.
To overcome this limitation, several approaches have been proposed; for
example in the frequency-dependent BTE (FD-BTE), the BZ is made isotropic from
a slice taken along a high-symmetry path, and then sampled in frequency space
Minnich  _et al._ (2011); Singh  _et al._ (2011); Loy  _et al._ (2013). In a
previous work, we developed a formalism, called the MFP-BTE, where the
sampling is carried out in the MFP-space Romano and Grossman (2015); Romano
and Kolpak (2017). Recently, a deterministic approach based on the self-
adjoint version of the BTE has been proposed Harter  _et al._ (2019). While
these methods initially considered isotropic dispersions, they can be in
principle extended to arbitrary anisotropic dispersions, and unnravel the rich
physics of the coupling between the crystal structure and nanoscale
geometries, as captured by several MonteCarlo simulations Vermeersch  _et al._
(2016); Wu  _et al._ (2016); Mei  _et al._ (2014); Landon and
Hadjiconstantinou (2014). Furthermore, these approaches, unlike single-MFP
methods, have allowed for incorporating first-principle data Romano and
Grossman (2014); Hao  _et al._ (2009); Romano  _et al._ (2012). Building upon
the MFP-BTE formalism, we propose a technique to solve efficiently the BTE
over the whole BZ, basing on the fact that, in the relaxation-time-
approximation (RTA), the nonequilibrium phonon populations are a smooth
function of the their vectorial MFPs; we exploit this property by
interpolating them onto a regular spherical grid, reducing dramatically
computing times while not compromising on the accuracy. Crucially, our method
computes heat transport in constant time with respect to the number of phonon
branches and wave vectors, opening up the possibility of simulating complex
unit-cell materials, such as Bi2Te3 and SnSe- based nanostructures. In
deriving our approach, we also introduce the spherical and polar
representations of the thermal conductivity and phonon suppression function,
both pivotal to unraveling the coupling between the geometry and the
underlying material. As an example, we first apply our method to Si membrane
with infinite thickness, obtaining a speed up of about 50x with respect to the
case with no interpolation. A potential speed up of about 3 is found for a
three-dimensioanal membrane, which is bound to increase for more complex
materials. The proposed method in practice extends the MFP-BTE to anistropic
materials, thus we call it the anisotropic MFP-BTE, or aMFP-BTE. Taken
together, we expect our approach to significantly expands the engineering
space of materials for thermoelectrics and heat management applications, while
keeping the computational effort amenable to most laptops.

##  II The Mode-Resolved Boltzmann transport equation

In absence of perturbation, e.g. an applied temperature gradient, the phonon
populations are in equilibrium with temperature $T_{0}$ and are given by the
Bose-Einstein distribution

$\bar{n}_{p}(\mathbf{q})=\frac{1}{e^{\frac{\hbar\omega_{p}(\mathbf{q})}{k_{B}T_{0}}}-1},$
(1)

where $p$ is the phonon branch and $\mathbf{q}$ the wave vector. When a
temperature gradient is applied, the phonon distributions depart from the
equilibrium ones, $n_{p}(\mathbf{r},\mathbf{q})=\bar{n}_{p}(\mathbf{q})+\Delta
n_{p}(\mathbf{r},\mathbf{q})$, and non-zero current develops. At the steady
state, the deviational distributions are obtained by the linearized, time-
independent BTE Ziman (2001),

$\displaystyle-$ $\displaystyle\mathbf{v}_{p}(\mathbf{q})\cdot\nabla\Delta
n_{p}(\mathbf{r},\mathbf{q})=$ (2) $\displaystyle=$
$\displaystyle\sum_{p^{\prime}}\int_{\mathrm{BZ}}\frac{d\mathbf{q}^{\prime}}{\Omega_{BZ}}A_{pp^{\prime}}(\mathbf{q},\mathbf{q}^{\prime})\Delta
n_{p^{\prime}}(\mathbf{r},\mathbf{q}^{\prime}),$

where $A_{pp^{\prime}}(\mathbf{q},\mathbf{q}^{\prime})$ is the scattering
operator and $\mathbf{v}_{p}(\mathbf{q})$ is the group velocity. Equation 2
has six unknowns, three for space and three for momentum. The total number of
branches is $N_{p}$. Let us assume that the BZ has been discretized uniformnly
into $N_{q}$ control volumes $\Delta\mathbf{q}_{k}$; we can then discretize
the momentum space by simply integrating both sides over
$\Delta\mathbf{q}_{k}$ and $\Delta\mathbf{q}_{k^{\prime}}$. Assuming that
distributions and group velocities are constant within the control volume, Eq.
2 becomes

$-\mathbf{v}_{\mu}\cdot\nabla\Delta
n_{\mu}(\mathbf{r})=\sum_{\mu^{\prime}}A_{\mu\mu^{\prime}}(\mathbf{r})\Delta
n_{\mu^{\prime}}(\mathbf{r}),$ (3)

where $\mu$ collectively indicates phonon branches and wave vectors. We
conveniently choose to work in the temperature formulation of the BTE Romano
and Grossman (2015); Romano (2020)

$-\mathbf{S}_{\mu}\cdot\nabla\Delta T_{\mu}=\sum_{\nu}W_{\mu\nu}\Delta
T_{\nu},$ (4)

where $\Delta T_{\mu}=\Delta n_{\mu}\hbar\omega_{\mu}/C_{\mu}$ are the phonon
pseudotemperatures (or simply “temperatures” hereafter),
$W_{\mu\mu^{\prime}}=A_{\mu\mu^{\prime}}\omega_{\mu}C_{\mu^{\prime}}C_{\mu}^{-1}\omega_{\mu^{\prime}}^{-1}$,
and $\mathbf{S}=\mathbf{v}_{\mu}C_{\mu}$; the term
$C_{\mu}=k_{B}(\eta\sinh{\eta})^{-2}$ is the mode-resolved heat capacity, with
$\eta=k_{B}\hbar\omega_{\mu}/(2k_{B}T_{0})$. Within the temperature
formulation, the heat flux is $\mathbf{J}=\sum_{\mu}C_{\mu}\Delta
T_{\mu}\mathbf{v}_{\mu}$.

In this work, we employ the relaxation time approximation (RTA), which is
accurate for many relevant materials, including silicon Esfarjani  _et al._
(2011). Within RTA, the scattering operator is $W_{\mu\nu}\approx
C_{\nu}\tau_{\nu}^{-1}\left(\delta_{\mu\nu}\Delta T_{\nu}-\Delta
T^{L}\right)$, where $\Delta T^{L}$ is a local pseudotemperature. This
quantity is computed by setting $\nabla\cdot\mathbf{J}=0$ in Eq. 4, leading to
the scattering operator
$W_{\mu\nu}=C_{\nu}/\tau_{\nu}\left(\delta_{\mu\nu}-a_{\nu}\right)$, where
$a_{\nu}=C_{\nu}/\tau_{\nu}\left[\sum_{k}C_{k}/\tau_{k}\right]^{-1}$. Equation
4 then simplifies to

$\mathbf{F}_{\mu}\cdot\nabla\Delta T_{\mu}^{(n)}+\Delta
T_{\mu}^{(n)}=\sum_{\nu}a_{\nu}\Delta T_{\nu}^{(n-1)},$ (5)

where $\mathbf{F}_{\mu}=\tau_{\mu}\mathbf{v}_{\mu}$ is the vectorial MFP.
Equation 5 is solved iteratively, with $\Delta T_{\nu}^{(0)}$ given by the
diffusive equation Romano and Grossman (2015).

We consider periodic materials along both $x$\- and $y$\- axis, whereas heat
flux is enforced by applying difference of temperature $\Delta
T_{\mathrm{ext}}=\Delta T_{\mu}^{\mathrm{R}}-\Delta T_{\mu}^{\mathrm{L}}$
along the $x$-axis. With $L$ being the distance between the right (R) and left
(L) contact, the effective thermal conductivity is computed by

$\kappa^{\mathrm{eff}}=-\frac{L}{\Delta
T_{\mathrm{ext}}A_{\mathrm{R}}}\int_{\mathrm{R}}\mathbf{J}\cdot\mathbf{\hat{n}}dS,$
(6)

where $A_{\mathrm{R}}$ is the area of the right contact. Along the wall we
apply the total diffuse scattering boundary condition, e.g. phonons outgoing
from a pore are fixed to the boundary temperature Landon and Hadjiconstantinou
(2014)

$T_{B}=\sum_{\mu}\Delta
T_{\mu}\mathrm{ReLu}(\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})\left[\sum_{\nu}\mathrm{ReLu}(\mathbf{S}_{\nu}\cdot\mathbf{\hat{n}})\right]^{-1},$
(7)

where $\mathrm{ReLu}(x)$ is

$\mathrm{ReLu}(x)=\begin{cases}x,&\text{if }x>0\\\
0,&\text{otherwise}.\end{cases}$ (8)

It’s straightforward to show that Eq. 7 satisfies the condition
$\sum_{\mu}\mathbf{J}_{\mu}\cdot\mathbf{\hat{n}}=0$. Equation 5 along with
boundary conditions is the mode-resolved BTE under RTA (MR-BTE), and can be
discretized using the finite-volume techniques, with $\tau_{\mu}$, $C_{\mu}$
and $\mathbf{v}_{\mu}$ obtained from first-principles calculations. In
practice, however, solving directly Eq. 5 is challenging: It requires solving
an algebraic equation as many times as all the combinations of phonon branch
and phonon vectors ($N_{p}N_{q}$), which can easily lead to prohibitive
computational load. To overcome this issue, we present the anisotropic MFP-BTE
(aMFP-BTE), described in the next section.

##  III The anisotropic MFP-BTE

At each iteration, the mode-resolved RTA-BTE, encoded in Eq. 5, can be
computed more efficiently if we exploit the fact that the temperatures are a
smooth function in the vectorial MFPs, $\mathbf{F}_{\mu}$; in fact, we have
equations of the form
$f(\mathbf{F},\mathbf{r})=g(\mathbf{r})-\mathbf{F}\cdot\nabla
f(\mathbf{F},\mathbf{r})$, with $g(\mathbf{r})$ being associated to the local
temperature. We can, therefore, solve Eq 5 for vectorial MFPs located on a
uniform grid and then retrieve $\Delta T_{\mu}$ by interpolation. We choose a
spherical grid with nodes $\mathbf{F}_{mlk}$, where $m$, $l$ and $k$ label the
magnitude ($\Lambda_{m}$), azimuthal ($\theta_{l}$) and polar angles
($\phi_{k}$), respectively. The polar and azimuthal angles span uniformly in
linear scale, while the magnitudes spread on log scale.

The generic vectorial MFP is

$\displaystyle\mathbf{F}_{mkl}=\Lambda_{m}\bigg{[}\sin{\phi_{k}}\sin{\theta_{l}}\mathbf{\hat{x}}$
$\displaystyle+$
$\displaystyle\cos{\phi_{k}}\sin{\theta_{l}}\mathbf{\hat{y}}+$ (9)
$\displaystyle+$
$\displaystyle\cos{\theta_{l}}\mathbf{\hat{z}}\bigg{]}=\Lambda_{m}\mathbf{S}_{kl},$

where $m$, $l$ and $k$ run up to $N_{\Lambda}$, $N_{\theta}$ and $N_{\phi}$,
respectively. The generic mode-resolved phonon temperature is

$\Delta T_{\mu}=\sum_{mlk}c_{mlk}^{\mu}\Delta T_{mlk},$ (10)

where $c_{mlk}^{\mu}$ are linear coefficients such that

$\mathbf{F}_{\mu}=\sum_{mlk}c_{mlk}^{\mu}\mathbf{F}_{mlk}.$ (11)

Within this formalism, the temperature is

$\Delta T^{L}=\sum_{\nu}a_{\nu}\sum_{mlk}a_{mlk}^{\mu}\Delta T_{mlk};$ (12)

finally, Eq. 5 becomes

$\displaystyle~{}\Lambda_{m}\mathbf{S}_{kl}\cdot\nabla\Delta T_{mkl}^{(n)}$
$\displaystyle+$ $\displaystyle\Delta T_{mkl}^{(n)}=$ (13) $\displaystyle=$
$\displaystyle\sum_{m^{\prime}k^{\prime}l^{\prime}}a_{m^{\prime}k^{\prime}l^{\prime}}\Delta
T_{m^{\prime}k^{\prime}l^{\prime}}^{(n-1)},$

where $a_{mkl}=\sum_{u}c_{mkl}^{\mu}$. Similarly, heat flux is given by
$\mathbf{J}=\sum_{mkl}\Delta T_{mkl}\mathbf{G}_{mkl}$, where
$\mathbf{G}_{mkl}=\sum_{\mu}c_{mkl}^{\mu}C_{\mu}\tau_{\mu}^{-1}\mathbf{F}_{\mu}$.
Upon convergence, the effective thermal conductivity is now provided in terms
of $\Delta T_{mkl}$, i.e.

$\kappa_{\mathrm{eff}}=\sum_{\mu}\frac{C_{\mu}}{\tau_{\mu}}F_{\mu,x}\langle\Delta
T_{\mu}\rangle=\sum_{mkl}G_{mkl,x}\langle\Delta T_{mkl}\rangle,$ (14)

where $\langle f\rangle=-L\Delta
T_{\mathrm{ext}}^{-1}A_{\mathrm{R}}^{-1}\int_{\mathrm{R}}fdS$. To better
understand the effect of the geometry on $\kappa_{\mathrm{eff}}$, it is
convenient to define a suppression function, which is a measure on how much
heat is carried in the nanomaterial compared to that from the bulk.
Originally, this tool was conceived as a MFP- or frequency-dependent function
Minnich (2012) and later was generalized to include directionality Romano and
Kolpak (2017). Here we define the mode-resolved suppression function as
$S_{\mu}=\langle\Delta T_{\mu}\rangle\left(v_{\mu,x}\tau_{\mu}\right)^{-1}$,
which is calculated with respect to the mode-specific vectorial MFP projected
onto the direction of the applied temperature. Using $S_{\mu}$, the effective
thermal conductivity is
$\kappa^{\mathrm{eff}}=\sum_{\mu}\kappa_{\mu}^{xx}S_{\mu}=\sum_{mkl}\kappa_{mkl}S_{mkl}$,
where

$\kappa_{mkl}=\sum_{\mu}\frac{C_{\mu}}{\tau_{\mu}}F_{\mu,x}a_{mkl}^{\mu}F_{mkl,x}$
(15) $S_{mkl}=\frac{\langle\Delta T_{mkl}\rangle}{F_{mkl,x}},$ (16)

are defined as the spherical bulk thermal conductivity and spherical
suppression function, respectively. Equation 13 along with the definitions
15-16 constitute the aMFP-BTE model. The corresponding workflow is summerized
in algorithm 1.

Algorithm 1 aMFP-BTE. Note that the $mkl$ labels have been grouped into $g$.

1:Solve $\Delta T^{L,(0)}\leftarrow\nabla^{2}\Delta T^{L,(0)}=0$

2:$n=1$

3:while error $>1e^{-3}$ do

4: $\Delta T^{L,(n+1)}=0$

5: $\kappa^{\mathrm{eff},(n+1)}=0$

6: for g = 1:$N_{\Lambda}N_{\theta}N_{\phi}$ do

7: $\Delta T_{g}^{(n+1)}\leftarrow\mathbf{F}_{g}\cdot\nabla\Delta
T_{g}^{(n+1)}+\Delta T_{g}^{(n+1)}=\Delta T^{(0)}$

8: $\Delta T^{L,(n+1)}\leftarrow T^{L,(n+1)}+a_{g}\Delta T_{g}^{(n+1)}$

9: $S_{g}\leftarrow\langle\Delta T_{g}^{(n+1)}\rangle F_{g,x}^{-1}$

10:
$\kappa^{\mathrm{eff},(n+1)}\leftarrow\kappa^{\mathrm{eff},(n+1)}+\kappa_{g}S_{g}$

11: end for

12: error
$\leftarrow|(\kappa^{\mathrm{eff},(n+1)}-\kappa^{\mathrm{eff},(n)})|/\kappa^{\mathrm{eff},(n+1)}$

13: $n\leftarrow n+1$

14:end while

##  IV Membranes with Infinite Thickness

We first apply the aMFP-BTE to porous Si membrane with circular pores and
infinite thickness. The simulation domain comprises a square unit cell of size
$L$, containing one circular pore, to which a difference of temperature
$\Delta T_{\mathrm{ext}}$ is applied. The chosen porosity is 0.2. The
scattering times at 300 K, group velocities and phonon frequencies are
computed with density functional theory and supercell approach, implemented in
AlmaBTE Carrete  _et al._ (2017), which also provides pre-computed second- and
third- order force constants. Convergence was found with a uniform wave vector
grid of 32x32x32 points, to which it corresponds a bulk thermal conductivity,
$\kappa_{\mathrm{bulk}}$, of 160 Wm-1K-1. Naturally occurring isotope disorder
was included in the calculation. At the macroscopic level, reduction in
thermal transport is given by Fourier’s law, which, in our case, gives
$\kappa_{f}$=107.4 W m-1K-1, in line with the prediction from Eucken theory
Hasselman and Johnson (1987), i.e.
$\kappa_{f}\approx\kappa_{\mathrm{bulk}}\frac{1-\phi}{1+\phi}$=106.7 W m-1K-1.

Figure 1: The effective thermal conductivity versus the number of iterations
for the case with L = 200 nm. Results are shown for the MR-BTE and the aMFP-
BTE, for different spherical discretizations. In this work, we choose
$N_{\Lambda}$ = 40 and $N_{\phi}$ = 96.

The value for $\kappa_{f}$ does not vary with $L$, as long as the porosity is
kept constant. On the other side, phonon size effects, computed by the BTE,
depend on the size of the unit cell, as shown below.

((a))

((b))

((c))

((d))

Figure 2: a) The effective thermal conductivity of porous Si with porosity 0.2
and periodicity L = 10 nm, 200 nm and 30 $\mu$m. In the inset, a map of the
magnitude of thermal flux (normalized to 1) b) The average phonon suppression
function for different periodicities. In the inset, a 3D represention of
$S_{tk}$ for the case with L = 10 nm. c) Mode-resolved
$\kappa^{\mathrm{eff}}_{\mu}$ versus the projected MFP
$\mathbf{F}_{\mu}\cdot\mathbf{\hat{x}}$ for different periodicities. d) Upper
panel: polar phonon suppression function, $S_{tk}$. A difference of
temperature is imposed along the $x$-axis. Bottom panel: polar effective
thermal conductivity, $\kappa_{\mu}^{\mathrm{eff}}$.

In the case with infinite thickness, the BTE is a two-dimensional system but
still with the azimuthal angle explicitly taken into account, i.e.

$\displaystyle\Lambda_{m}$
$\displaystyle\sin(\theta_{l})\bigg{[}\sin(\phi_{k})\frac{\partial\Delta
T_{mkl}}{\partial x}+$ (17) $\displaystyle+$
$\displaystyle\cos(\phi_{k})\frac{\partial\Delta T_{mkl}}{\partial
y}\bigg{]}+\Delta T_{mkl}=\Delta T_{L}.$

Unlike with the 3D case, the phonon temperatures are a smooth function of
$\Lambda_{m}\sin(\theta_{l})$, which then can be reformulated as an
“effective” MFP $\xi_{t}$. Leveraging this result, the interpolation is
performed between the bulk MFPs projected onto the xy plane, and

$\mathbf{F}_{tk}=\xi_{t}\left[\sin{(\phi_{k})}\mathbf{\hat{x}}+\cos{(\phi_{k})}\mathbf{\hat{y}}\right],$
(18)

with the corresponding interpolation coefficients being $c_{tk}^{\mu}$. The
polar representations of the bulk thermal conductivity and suppression
function are
$\kappa_{tk}=\sum_{\mu}C_{\mu}/\tau_{\mu}F_{\mu,x}c_{tk}^{\mu}F_{\mu,x}$ and
$S_{mt}=\langle\Delta T_{tk}\rangle/F_{tk,x}$, respectively. Since only the
discretization in the polar angle and MFPs is needed, this treatment greatly
enhances the computational efficiency, unlocking realistic phonon transport
simulations on common laptops.

To assess the efficiency of the proposed model, we solve the MR-BTE, i.e. Eq.
4, and the aMFP-BTE, defined by Eq. 13. The mode resolution is the same as the
q-grid used for bulk thermal conductivity calculations, $N_{p}N_{q}$ = 196608.
On the other hand, we choose a variable spherical grid for the aMFP-BTE,
obtaining an agreement between the two models within $0.5\%$ with
$N_{\Lambda}$ = 40 and $N_{\phi}$ = 96, as shown in Fig. 1. We note that for
the aMFP-BTE case, the “interpolated modes” were only 3840, obtaining roughly
a 50x speed up with respect to the MR-BTE.

Figure 2(a) shows $\kappa^{\mathrm{eff}}$ for L = 50 nm, 200 nm and 30 $\mu$m.
The last case employs a multiscale approach, which will be discussed
elsewhere. As expected, strong suppressions are achieved for smaller
periodicities since the distance between the pores’ walls shrinks and phonon
transport is suppressed over a wider range of the MFP distribution Romano and
Grossman (2015). For large L, we approach the diffusive limit, described
above. As shown in the inset of 2(a), thermal flux is concentrated in the
space between the pores, a well-known signature of ballistic transport
Anufriev and Nomura (2020).

The influence of the structure on thermal transport is best described by the
suppression function, as detailed in the previous section. Its angular
average, given by $S_{t}=\sum_{k}S_{tk}$ and illustrated in Fig. 2(b), shows
three regimes Minnich (2012): (i) the large-MFP regime, namely where transport
is mostly ballistic, $S_{t}\approx\xi_{t}^{-1}$; (ii) the intermediate-MFP
regime, where both ballistic and diffusive transport are present and (iii) the
small-MFP regime, where transport is mainly diffusive. Although all the
configurations have the same porosity, their small-MFP limits do not match the
prediction from Fourier’s law. In fact, in nongray materials, ballistic
phonons may also effect those who travel diffusively via the definition of
$\Delta T^{L}$ Romano and Kolpak (2019). For structures whose characteristic
length is larger than most heat-carrying phonons, however,
$S_{0}\approx\kappa_{\mathrm{fourier}}/\kappa_{\mathrm{bulk}}$. Once $\Delta
T_{tk}$ are computed, it is possible to compute the mode-resolved temperature
$\Delta T_{\mu}=\sum_{tk}c_{tk}^{\mu}\Delta T_{tk}$; from the mode-resolved
temperature, we can then obtain $S_{\mu}$ and thus the mode resolved effective
thermal conductivity
$\kappa^{\mathrm{eff}}_{\mu}=C_{\mu}\mathbf{v}_{\mu}F_{\mu,x}S_{\mu}$. In Fig.
2(c), we plot $\kappa^{\mathrm{eff}}_{\mu}$ versus $\mathbf{F}_{\mu,x}$ for
different L; we note that, as L becomes smaller, a wider spectrum of the MFP
is suppressed.

The aMFP-BTE allows to explore the coupling between the structure and the
material at the angular level. In the case with infinite thickness, only polar
discretization is needed thus we can plot $S_{tk}$ as a surface in 3D, as
shown in the inset of Fig. 2(b) for the case with L = 50 nm; we note there are
two main lobes, corresponding to the forward and backward direct paths Romano
and Kolpak (2017). In Fig. 2(d), we show the top-view of the polar suppression
function for all L as well the polar representation of
$\kappa^{\mathrm{eff}}$. For L = 10 nm, $S_{tk}$ is strongly anisotropic
showing four additional lobes. For L = 200, however, this anisotropy becomes
less pronounced, and for L = 30 $\mu$m, $S_{tk}$ is mostly isotropic. This
last case can be regarded as the diffusive limit. Note that for
$\mathbf{F}_{tk}$ such that
$|\mathbf{F}_{tk}\cdot\mathbf{\hat{y}}|/|\mathbf{F}_{tk}|$ is small,
numerically instabilities regarding the calculation of $S_{tk}$ occur. To
understand this issue, it is convenient to rewrite the temperatures as $\Delta
T_{tk}=\Delta T^{L}-\mathbf{F}_{tk}\cdot\nabla T_{tk}$. The suppression
function then reads

$S_{tk}=-\frac{\partial\langle\Delta T_{tk}\rangle}{\partial
x}-\frac{1}{\tan{\phi_{k}}}\frac{\partial\langle\Delta T_{tk}\rangle}{\partial
y};$ (19)

the second term of Eq. 19, for $\phi_{k}\approx 0$ and $\phi_{k}\approx\pi$,
becomes numerically challenging. For this reason, in Fig. 2(c), we plot
$S_{tk}$ only for regions far away from those two critical cones. We note
that, however, these inaccuracies do not alter $\kappa^{\mathrm{eff}}$ since
$\kappa^{\mathrm{bulk}}_{\mu}\approx 0$ for $\mathbf{F}_{\mu}$ aligned with
$\mathbf{\hat{y}}$ and -$\mathbf{\hat{y}}$. Furthermore, $S_{tk}$ is computed
only for analysis purposes, while in the actual calculation for
$\kappa^{\mathrm{eff}}$ we use Eq. 14, which overcomes these inaccuracies.

Similarly to $S_{tk}$, we explore the polar representation of the effective
thermal conductivity, i.e.
$\kappa^{\mathrm{eff}}_{tk}=\kappa^{\mathrm{bulk}}_{tk}S_{tk}$, shown in Fig.
2(d). For all L, we note that most heat is carried by phonons along high-
symmetry axis, and that the anisotropy of $\kappa^{\mathrm{eff}}_{tk}$ is
influenced by $L$. In fact, for the case with L = 200 nm, we can also observe
a circular band in addition to peaks around high-symmetry axis. Furthermore,
we note that the case with L = 10 nm, the extra lobes in $S_{tk}$ are aligned
with the 100 axis of the crystal, thus capturing heat transport in that region
of the polar space. Lastly, for L = 30 $\mu$m,
$\kappa^{\mathrm{eff}}_{tk}\approx\kappa_{tk}$; this filtering effect is made
accessible thanks to the aMFP-BTE and can be exploited to enhance thermal
transport tunability.

##  V Membranes with Finite Thickness

For membranes with finite thickness, thermal transport is further reduced due
to the top and bottom surfaces scattering phonons. In this case, the
characteristic size is determined by both pore-pore distance and thickness. To
show this effect, we consider a periodic membrane with circular pores,
porosity 0.2, $L$ = 50 nm and thickness t = 10 nm. Convergence was found with
$N_{\Lambda}$=30, $N_{\phi}$ = 48 and $N_{\theta}$ = 48, at which it
corresponds $\kappa^{\mathrm{eff}}\approx$ 11.3 W m-1K-1. The case with
infinite thickness (t=$\infty$), modeled as described in the previous section,
gives 23.7W m-1K-1. The mode-resolved $\kappa^{\mathrm{eff}}$, reported in
Fig. 3, reveals that most modes have stronger suppression with respect to the
case with t=$\infty$; these modes are those who scatter with the top and
bottom surfaces. On the other side, there still a small fraction of phonons
that travel without being strongly influenced by the top and bottom surfaces,
i.e. those with vectorial MFPs that are aligned with the $x$-axis. In the
inset of Fig. 3, we report a 3D map of the magnitude of thermal flux, which,
analogously to the case with $t=\infty$, shows high values in the space
between pores.

Figure 3: Mode-resolved effective thermal conductivity,
$\kappa_{\mu}^{\mathrm{eff}}$, versus the projected MFP
$\mathbf{F}_{\mu}\cdot\mathbf{\hat{x}}$, for bulk and for thicknesses t = 10
nm and $\infty$. In the inset, a 3D map of the magnitude of thermal flux,
normalized to 1.

Lastly, we note that in this case the speed gain with respect to the MR-BTE is
around 3, assuming that convergences occur with comparable number of steps.
However, with more complex materials, the number of branches as well as the
number of wave vectors increases, with little or no increase in the
computational demand from aMFP-BTE; therefore, the aMFP-BTE may introduce
large speed ups for 3D systems, as well.

##  VI Conclusion

We introduce the aMFP-BTE, a method that reduces dramatically the
computational demand of thermal conductivity calculations in nanostructures,
leveraging interpolation of phonon distributions in the vectorial MFP space.
The application of the aMFP-BTE to Si porous membranes with infinite thickness
reveals a speed-up of about 50x. Taking into account the whole BZ, our method
is able to elucidate on the interplay between the anisotropy of the bulk
thermal conductivity and the material’s geometry, highlighting novel
mechanisms for thermal transport tuning. The aMFP-BTE solves for thermal
transport in constant time with respect to the number of phonon branches and
wave vectors, opening up the possibility of simulating phonon size effects in
promising thermoelectric materials, such as SnSe and Bi2Te3. The code used for
this work will be released with the package OpenBTE Romano (2019).

###### Acknowledgements.

Research was partially supported by the Solid-State SolarThermal Energy
Conversion Center (S3TEC), an Energy Frontier Research Center funded by the
U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences
(BES), under Award No. DESC0001. The author thanks Steven G. Johnson for
helpful discussions.

##  Appendix A Finite-Volume Formulation of the BTE

We discretize Eqs. 2-13 using the finite-volume techniques. For simplicity, we
label the phonon temperatures with $\mu$, but the same procedure holds for
both the MR-BTE and the aMFP-BTE. The mesh, generated by GMSH Geuzaine and
Remacle (2009), is unstructured, with elements at the border of domain
matching their periodic counterparts. The BTE for a generic $\mathbf{F}_{\mu}$
reads

$\mathbf{F}_{\mu}\cdot\nabla\Delta T_{\mu}(\mathbf{r})+\Delta
T_{\mu}(\mathbf{r})=\Delta T^{L}(\mathbf{r}),$ (20)

which, after integrating both side over the control volume $V_{c}$, becomes

$\frac{1}{V_{c}}\int_{\partial V_{c}}\Delta
T_{\mu}\mathbf{F}\cdot\mathbf{\hat{n}}dS+\Delta T_{c\mu}=\Delta T^{L}_{c},$
(21)

where we uses Gauss’ theorem. The RHS of Eq. 21 is simply given by $\Delta
T^{L}_{c}=\sum_{\nu}a_{\nu}\Delta T_{c\nu}$. The term associated to the
surface integral is discretized according to the upwind method Murthy and
Mathur (1998); Romano and Di Carlo (2011), yielding

$\frac{1}{V_{c}}\int_{\partial V_{c}}\Delta
T\mathbf{F}_{\mu}\cdot\mathbf{\hat{n}}dS=I_{c\mu}^{\mathrm{in}}+I_{c\mu}^{\mathrm{out}}+I_{c\mu}^{\mathrm{B}}+I_{c\mu}^{\mathrm{P}}.$
(22)

The contribution $I_{c\mu}^{\mathrm{in}}$ accounts for all the flux incoming
to a given element, and is given by

$I_{c\mu}^{\mathrm{in}}=-\sum_{c^{\prime}}\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{R}_{cc^{\prime}})\Delta
T_{c^{\prime}\mu},$ (23)

where

$\mathbf{R}_{cc^{\prime}}=\begin{cases}\frac{1}{V_{c}}A_{cc^{\prime}}\mathbf{\hat{n}}_{cc^{\prime}}^{\prime},&\text{if
}c\text{ and }c^{\prime}\text{ are neighbors}\\\
0,&\text{otherwise}.\end{cases}$ (24)

The terms $\mathbf{\hat{n}}_{cc^{\prime}}$ and $A_{cc^{\prime}}$ are the
normal (pointing toward the volume $c^{\prime}$) and the area of the side
between the volume $c$ and $c^{\prime}$. The ReLu function filters only the
incoming contributions to the element $c$, while $\mathbf{R}_{cc^{\prime}}$
represents the connections between the volumes. The contribution
$I_{c\mu}^{\mathrm{out}}$ is the flux leaving the volume, and reads

$I_{c\mu}^{\mathrm{out}}=\sum_{c^{\prime\prime}}\mathrm{ReLu}(\mathbf{F}_{\mu}\cdot\mathbf{R}_{cc^{\prime\prime}})\Delta
T_{c^{\prime\prime}}.$ (25)

In the upwind scheme, this term amounts to the diagonal of the stiffness
matrix. The third term of Eq. 22 arises from the flux bouncing back from an
adiabatic boundary, and is given by

$I_{c\mu}^{\mathrm{B}}=-\sum_{s}\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{\hat{n}}_{s})g_{sc}\Delta
T_{s}.$ (26)

where $\mathbf{\hat{n}_{s}}$ is the normal of the surface pointing outward
with respect to the computational domain, and $T_{s}$ is the boundary
temperature, computed with Eq. 7. The term $g_{sc}$ is

$g_{sc}=\begin{cases}\frac{1}{V_{c}}A_{s},&\text{if }s\text{ is a side of
}c\\\ 0,&\text{otherwise},\end{cases}$ (27)

where $A_{s}$ is the area of the side $s$. Using Eq. 7, we can rewrite the
boundary contribution as $I_{c\mu}^{\mathrm{B}}=\sum_{\nu}H_{\mu\nu c}\Delta
T_{c\nu}$, where

$H_{c\mu\nu}=-\sum_{s}\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{\hat{n}}_{s})g_{sc}\frac{\mathrm{ReLu}(\mathbf{S}_{\nu}\cdot\mathbf{\hat{n}}_{s})}{\sum_{k}\mathrm{ReLu}(\mathbf{S}_{k}\cdot\mathbf{\hat{n}}_{s})}.$
(28)

The fourth and last contribution is due to the periodic boundary conditions,

$I_{c\nu}^{\mathrm{P}}=-\sum_{c^{\prime}}\mathrm{ReLu}(-\mathbf{F}_{\nu}\cdot\mathbf{R}_{cc^{\prime}}^{P})\Delta
T_{\mathrm{ext}}$ (29)

where

$\mathbf{R}_{cc^{\prime}}^{P}=\begin{cases}\frac{1}{V_{c}}A_{cc^{\prime}}\mathbf{\hat{n}}_{cc^{\prime}},&\text{if
}c\text{ and }c^{\prime}\text{ are periodic}\\\
0,&\text{otherwise}.\end{cases}$ (30)

Putting these contributions together, along with Eq. 20, we obtain the
iterative linear system

$\sum_{c^{\prime}}A_{cc^{\prime}\mu}\Delta
T_{c^{\prime}\mu}^{(n)}=\sum_{\nu}(a_{\nu}+H_{c\mu\nu})\Delta
T_{c\nu}^{(n-1)}+I_{c\mu}^{P}$ (31)

where

$\displaystyle A_{cc^{\prime}\mu}$
$\displaystyle=\delta_{cc^{\prime}}+\sum_{c^{\prime\prime}}\mathrm{ReLu}(\mathbf{F}_{\mu}\cdot\mathbf{R}_{cc^{\prime\prime}})\delta_{cc^{\prime}}-$
(32)
$\displaystyle-\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{R}_{cc^{\prime\prime}})\delta_{c^{\prime}c^{\prime\prime}}.$

Finally, we can write Eq. 31 in vector notation

$\mathbf{A}_{\mu}\mathbf{\Delta
T_{\mu}}^{(n)}=\sum_{\nu}\mathbf{S}_{\mu\nu}\mathbf{\Delta
T}_{\nu}^{(n-1)}+\mathbf{I}_{\mu},$ (33)

where $S_{c\mu\nu}=\alpha_{\nu}+H_{c\mu\nu}$.

## References

  * Vineis  _et al._ (2010) C. J. Vineis, A. Shakouri, A. Majumdar, and M. G. Kanatzidis, Adv. Mater. 22, 3970 (2010).
  * Lorenzi  _et al._ (2018) B. Lorenzi, R. Dettori, M. T. Dunham, C. Melis, R. Tonini, L. Colombo, A. Sood, K. E. Goodson, and D. Narducci, J. Electron. Mater. 47, 5148 (2018).
  * Lee  _et al._ (2010) H. Lee, D. Vashaee, D. Wang, M. S. Dresselhaus, Z. Ren, and G. Chen, J. Appl. Phys. 107, 094308 (2010).
  * Kim  _et al._ (2007) W. Kim, R. Wang, and A. Majumdar, Nano Today 2, 40 (2007).
  * Moore and Shi (2014) A. L. Moore and L. Shi, Mater. Today 17, 163 (2014).
  * Anufriev  _et al._ (2017) R. Anufriev, A. Ramiere, J. Maire, and M. Nomura, Nat. Commun. 8, 15505 (2017).
  * Zhang  _et al._ (2013) K. Zhang, Y. Zhang, and S. Wang, Scientific reports 3, 1 (2013).
  * Chen (2005) G. Chen,  _Nanoscale energy transport and conversion: a parallel treatment of electrons, molecules, phonons, and photons_ (Oxford University Press, USA, 2005).
  * Murthy  _et al._ (2005) J. Y. Murthy, S. V. J. Narumanchi, P.-G. Jose’A, T. Wang, C. Ni, and S. R. Mathur, Int. J. Multiscale Com. 3, 5 (2005).
  * Ziman (2001) J. M. Ziman,  _Electrons and Phonons_ (Oxford University Press, 2001).
  * Murthy and Mathur (1998) J. Murthy and S. Mathur, J. Thermophys. Heat Trans. 12, 313 (1998).
  * Minnich  _et al._ (2011) A. J. Minnich, G. Chen, S. Mansoor, and B. S. Yilbas, Phys. Rev. B 84, 235207 (2011).
  * Singh  _et al._ (2011) D. Singh, J. Y. Murthy, and T. S. Fisher, J. Heat Transfer 133 (2011).
  * Loy  _et al._ (2013) J. M. Loy, J. Y. Murthy, and D. Singh, J. Heat Transfer 135 (2013), 10.1115/1.4007654.
  * Romano and Grossman (2015) G. Romano and J. C. Grossman, J. Heat Transf. 137, 71302 (2015).
  * Romano and Kolpak (2017) G. Romano and A. M. Kolpak, Sci. Rep. 7, 44379 (2017).
  * Harter  _et al._ (2019) J. R. Harter, S. A. Hosseini, T. S. Palmer, and P. A. Greaney, Int. J. Heat Mass Transf. 144, 118595 (2019).
  * Vermeersch  _et al._ (2016) B. Vermeersch, J. Carrete, and N. Mingo, Appl. Phys. Lett. 108, 193104 (2016).
  * Wu  _et al._ (2016) R. Wu, R. Hu, and X. Luo, J. Appl. Phys. 119, 145706 (2016).
  * Mei  _et al._ (2014) S. Mei, L. N. Maurer, Z. Aksamija, and I. Knezevic, J. Appl. Phys. 116, 164307 (2014).
  * Landon and Hadjiconstantinou (2014) C. D. Landon and N. G. Hadjiconstantinou, J. Appl. Phys. 116, 163502 (2014).
  * Romano and Grossman (2014) G. Romano and J. C. Grossman, Appl. Phys. Lett. 105, 33116 (2014).
  * Hao  _et al._ (2009) Q. Hao, G. Chen, and M.-S. Jeng, J. Appl. Phys. 106, 114321 (2009).
  * Romano  _et al._ (2012) G. Romano, A. Di Carlo, and J. C. Grossman, J. Comput. Electron. 11, 8 (2012).
  * Romano (2020) G. Romano, arXiv preprint arXiv:2002.08940 (2020).
  * Esfarjani  _et al._ (2011) K. Esfarjani, G. Chen, and H. T. Stokes, Phys. Rev. B 84, 085204 (2011).
  * Minnich (2012) A. J. Minnich, Phys. Rev. Lett. 109, 205901 (2012).
  * Carrete  _et al._ (2017) J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, T. Wang, G. K. Madsen, and N. Mingo, Comput. Phys. Commun. 220, 351 (2017).
  * Hasselman and Johnson (1987) D. Hasselman and L. F. Johnson, J. Compos. Mater. 21, 508 (1987).
  * Anufriev and Nomura (2020) R. Anufriev and M. Nomura, Mater. Today Phys. 15, 100272 (2020).
  * Romano and Kolpak (2019) G. Romano and A. M. Kolpak, J. Heat Transfer 141, 12401 (2019).
  * Romano (2019) G. Romano, in  _APS March Meeting Abstracts_, Vol. 2019 (2019) pp. T70–321.
  * Geuzaine and Remacle (2009) C. Geuzaine and J.-F. Remacle, Int. J. Numer. Meth. Eng. 79, 1309 (2009).
  * Romano and Di Carlo (2011) G. Romano and A. Di Carlo, IEEE Trans. Nanotechnol. 10, 1285 (2011).