# PaintNeSF: Artistic Creation of Stylized Scenes with Vectorized 3D Strokes

Hao-Bin Duan1 Miao Wang1,2* Yan-Xun Li1 Yong-Liang Yang3
1State Key Laboratory of Virtual Reality Technology and Systems, Beihang
University
2Zhongguanchun Laboratory 3Department of Computer Science, University of Bath
https://buaavrcg.github.io/PaintNeSF

###### Abstract

We present Paint Neural Stroke Field (PaintNeSF), a novel technique to
generate stylized images of a 3D scene at arbitrary novel views from multi-
view 2D images. Different from existing methods which apply stylization to
trained neural radiance fields at the voxel level, our approach draws
inspiration from image-to-painting methods, simulating the progressive
painting process of human artwork with vector strokes. We develop a palette of
stylized 3D strokes from basic primitives and splines, and consider the 3D
scene stylization task as a multi-view reconstruction process based on these
3D stroke primitives. Instead of directly searching for the parameters of
these 3D strokes, which would be too costly, we introduce a differentiable
renderer that allows optimizing stroke parameters using gradient descent, and
propose a training scheme to alleviate the vanishing gradient issue. The
extensive evaluation demonstrates that our approach effectively synthesizes 3D
scenes with significant geometric and aesthetic stylization while maintaining
a consistent appearance across different views. Our method can be further
integrated with style loss and image-text contrastive models to extend its
applications, including color transfer and text-driven 3D scene drawing.

Figure 1: We propose a method to stylize a 3D scene from multi-view 2D images
using vectorized 3D strokes based on geometric primitives and splines. The
four scenes from left to right are drawn with axis-aligned box, oriented box,
ellipsoid, and cubic Bézier curve, respectively.

††* Corresponding author.

##  1 Introduction

Artistic style image creation, historically a domain requiring significant
skill and time, has been revolutionized by neural network-based Style Transfer
techniques [7, 6]. These methods usually separate, manipulate, and merge the
content and style of images to create an artistic effect. Extending this to 3D
scenes, however, is challenging due to complex geometries and appearance
traits [27], and traditional convolutional neural network (CNN) methods [24]
are not readily adaptable to 3D spaces. Recently, 3D implicit scene
representations, particularly Neural Radiance Fields (NeRF), have become
popular, since they are fully differentiable and easy to optimize. Efforts to
apply artistic styles to NeRF [14, 21, 10] often involve separate NeRF
training and style transfer, followed by integration style transfer loss from
multiple views [5, 30, 3]. The fundamental approach continues to be the
independent manipulation of content and style. Nevertheless, existing NeRF
stylization techniques only alter the color aspects, leaving the density-and
consequently the geometry—unchanged, leading to stylizations that lack
geometric changes.

Style transfer methods in 2D and 3D primarily manipulate pixels or voxels to
achieve artistic effects, differing from traditional art where artists use
brushstrokes and a degree of randomness in brushes, materials, and colors.
This traditional approach is time-consuming and skill-intensive. To automate
it, research has explored image generation with vectorized strokes [9],
defined by position, color, size, and direction, and optimized to match target
images. 3D painting tools like Google’s Tilt Brush [4] use VR/AR to create art
scenes with various brushes in virtual spaces. However, the field of automated
stroke-based 3D art generation remains relatively unexplored.

In this work, we present a novel technique for transforming 2D images with
known poses into stylized 3D scenes using vectorized 3D strokes. Our approach
recreates 3D scenes that exhibit distinct geometric and appearance styles,
emulating the stroke-by-stroke painting process employed by human artists. The
vector-based nature of our system enables rendering at any desired resolution.
It accommodates a range of 3D stroke styles and supports more intricate
stylizations and generations through the use of style and semantic losses. The
core of our method is a differentiable, vectorized 3D scene representation,
which diverges from NeRF by using parameterized 3D strokes to synthesize 2D
images at different viewpoints. This differentiable representation allows for
the direct optimization of stroke parameters via gradient descent,
sidestepping the need for previous greedy-search-based or reinforcement
learning-based prediction methods in stroke-based image synthesis.
Furthermore, we analyze the gradient behaviors of our stroke-based 3D
representation and introduce a stroke initialization and update scheme that
avoids sub-optimal initialization where optimization is difficult due to local
minima.

Our method was evaluated using the multi-view datasets of real-world and
synthetic images. Our experiments demonstrate that it effectively creates
high-quality artistic 3D scene renderings by maintaining both global visual
fidelity and local color accuracy while ensuring perspective consistency. Our
contribution is summarized as follows:

  * •

We propose a novel method to translate multi-view 2D images into stylized 3D
scenes using 3D strokes based on basic primitives and spline curves whose
parameters can be learned through gradient descent.

  * •

We present a novel initialization scheme for the parameters in our stroke-
based representation to address issues of flat gradients during optimization.

  * •

Our method supports various 3D stroke styles and can be applied for
sophisticated geometry and appearance stylizations similar to those created by
human painters.

##  2 Related Work

###  2.1 Image Painting

Image painting, evolving over a long history, uses vectorized strokes in
various colors and styles for 2D image creation. Early methods like ”Paint by
Numbers” [8] introduced brush parameters such as position, color, size, and
direction, optimized through techniques like stroke zooming, noise adding, and
color enhancement. Subsequently, more expressive forms emerged, employing
slender spline curves [9] and rectangular strokes [28] to better capture
source image details. The development of brush libraries [25, 31] allowed for
more diverse styles in paintings, notably in oil painting effects. With the
advent of deep learning, techniques like generative neural networks [20] have
further refined style representation. Advanced loss based on optimal
transportation [32] and feed-forward paint transformer [15] are used to
improve fidelity and speed up the painting process. While these 2D image
synthesis techniques are relatively mature, applying styles directly to 2D
images from 3D scenes lacks consistency across views, which motivates our
present work on using 3D strokes.

###  2.2 Stylization of 3D Scenes

The process of 3D scene stylization involves applying the visual style of a
chosen reference to a target 3D scene while preserving its inherent 3D
structure. Traditional methods [26] faced challenges like limited control over
perceptual factors and the inability to selectively target stylization
effects. The advent of Neural Radiance Fields (NeRF) has provided a more
flexible representation for 3D scene stylization. Techniques like ARF [30]
transfer artistic features from 2D style images to 3D scenes, treating
stylization as an optimization problem within the NeRF framework. Chiang et
al. [3] combines NeRF’s implicit representation with a hypernetwork for style
transfer. StylizedNeRF [10] jointly learns 2D and 3D stylization, and StyleRF
[14] learns high-level features in 3D space for fast zero-shot style transfer.
Our approach diverges from existing methods, which rely on additional style
references and fundamentally contrast with the processes human artists use to
create artworks. We focus on generating artistic images through 3D brush
strokes, eliminating the need for style images. Our method more closely
resembles image painting techniques than voxel-level manipulations.

##  3 Methodology

In this section, we introduce our 3D scene stylization framework with stroke-
based representation. Our framework mainly consists of three parts: (1) a
collection of 3D strokes based on primitives and spline curves, (2)
differentiable rendering and composition of strokes, (3) a stroke
initialization and update scheme that stabilizes the training process.

###  3.1 Overview of Stroke Field

NeRF [18] represents the scene as the density
$\sigma(\mathbf{x})\in\mathbb{R}^{+}$ and RGB radiance
$\mathbf{c}(\mathbf{x},\mathbf{d})$, modeled by an MLP that takes the spatial
coordinates $\mathbf{x}\in\mathbb{R}^{3}$ and view directions
$\mathbf{d}\in\mathbb{R}^{2}$ as input. Given a camera pose, a batch of rays
$\mathbf{r}(t)=\mathbf{o}+t\mathbf{d}$ is sampled and cast from the camera’s
center $o\in\mathbb{R}^{3}$ along the direction $\mathbf{d}\in\mathbb{R}^{3}$
passing through each pixel of the image. The color of a ray $\mathbf{r}$ is
given by differentiable volume rendering:

$C(\mathbf{r})=\int_{t_{n}}^{t_{f}}\sigma(\mathbf{r}(t))\mathbf{c}(\mathbf{r}(t),\mathbf{d})\exp(\int_{t_{n}}^{t}\sigma(\mathbf{r}(s))ds)dt$
(1)

As shown in Fig. 2, while NeRF models the scene at per-voxel level with an
implicit field, our method represents the scene as a set of vectorized 3D
strokes. Similar to NeRF, the stroke field is defined by two spatially varying
functions for density $\sigma(\mathbf{x})\in\mathbb{R}^{+}$ and RGB color
$\mathbf{c}(\mathbf{x},\mathbf{d})$ and rendered into 2D images at a given
camera pose, using the same differentiable volume rendering formula as in
NeRF.

While the stroke field shares NeRF’s field definition, their core formulations
are fundamentally different. Each point’s density and color in the stroke
field are set by 3D strokes, shaped and styled to resemble brush traces in
human drawings. This is analogous to the difference between a rasterized image
and a vectorized image but in 3D space. Specifically, the geometry,
appearance, and opacity of 3D strokes are given by three parameters, including
shape $\mathbf{\theta_{s}}\in\mathbb{R}^{N\times d_{s}}$, color
$\mathbf{\theta_{c}}\in\mathbb{R}^{N\times d_{c}}$ and density
$\mathbf{\theta_{\sigma}}\in(\mathbb{R}^{+})^{N}$, where $N$ is the total
number of strokes in the field, $d_{s}$ and $d_{c}$ are the number of
parameters of the specific stroke shape and color, respectively. We combine
these strokes in a differentiable way to acquire the density and color field:

$(\sigma,\mathbf{c})=StrokeField(\mathbf{x},\mathbf{d};\,\mathbf{\theta_{s}},\mathbf{\theta_{c}},\mathbf{\theta_{\sigma}})$
(2)

Figure 2: Our method learns a vectorized stroke field instead of MLP-based
implicit representation to represent a 3D scene.

###  3.2 3D Strokes

We first define the shape of our 3D strokes in order to paint them into the 3D
space. The shape of a 3D stroke is essentially a volume region formed by a
closed two-dimensional surface. To define such volume, we use the Signed
Distance Field (SDF), a scalar value function $\mathrm{sdf}(\mathbf{x})\to
s\in\mathbb{R}$ that gives the signed distance from a point to the closed
surface to describe the stroke shape in 3D. The sub-space
$\\{\mathbf{p}\in\mathbb{R}^{3}\;|\;\mathrm{sdf}(\mathbf{p})\leq 0\\}$
enclosed by the zero level-set is the stroke volume. We construct two types of
strokes based on basic geometric primitives and spline curves respectively.

####  3.2.1 Basic Primitives

Our first category of 3D strokes is based on common geometric primitives such
as spheres, cubes, etc. For simplicity, we normalize these geometries at the
origin of the unit space, apply a transformation matrix to the unit
geometries, and acquire their SDFs in the scene space. The unit geometry can
optionally contain a shape parameter
$\theta_{s}^{\text{basic}}\in\mathbb{R}^{d_{\text{basic}}}$ describing its
transformation-agnostic shape feature, where $d_{\text{basic}}$ is the number
of basic shape parameters. The SDF of unit primitive is defined as
$\mathrm{sdf}_{\text{unit}}:(\mathbf{\hat{p}},\theta_{s}^{\text{basic}})\to
s\in\mathbb{R}$, where $\mathbf{\hat{p}}$ is a 3D point in the unit space. We
list the common unit geometric primitives in Tab. 1.

Primitive | Params |  SDF formula
---|---|---
Sphere | None |  $\|\mathbf{p}\|_{2}-1$
Cube | None |  $\min(\max(\mathbf{q}_{x},\mathbf{q}_{y},\mathbf{q}_{z}),0)+\|\max(\mathbf{q},\mathbf{0})\|_{2}$, where $\mathbf{q}=|\mathbf{p}|-\mathbf{1}$
Tetrahedron | None |  $(\max(|\mathbf{p}_{x}+\mathbf{p}_{y}|-\mathbf{p}_{z},|\mathbf{p}_{x}+\mathbf{p}_{y}|+\mathbf{p}_{z})-1)/\sqrt{3}$
Octahedron | None |  $(\|\mathbf{p}\|_{1}-1)/\sqrt{3}$
Round Cube | $r$ |  $\min(\max(\mathbf{p}_{x},\mathbf{p}_{y},\mathbf{p}_{z}),0)+\|\max(\mathbf{p},\mathbf{0})\|_{2}-r$
Triprism | $h$ |  $\max(|\mathbf{p}_{y}|-h,\max(|\mathbf{p}_{x}|*\sqrt{3}/2+\mathbf{p}_{z}/2,-\mathbf{p}_{z})-0.5)$
Capsule Line | $h,r_{\delta}$ |  $\|\mathbf{p}-[0,\min(\max(\mathbf{p}_{y},-h),h),0]\|_{2}-r_{\delta}\min(\max((0.5*(\mathbf{p}_{y}+h)/h,0),1)-1$

Table 1: SDF formula for common unit geometry primitives.

We apply a transformation matrix $\mathbf{T}\in\mathbb{R}^{4\times 4}$ to
transform primitives from the unit coordinates $\mathbf{\hat{p}}$ to the
shared scene coordinates $\mathbf{p}=\mathbf{T}\mathbf{\hat{p}}$. The
parameters of transformation are composed of a translation
$\mathbf{t}\in\mathbb{R}^{3}$, a rotation described by Euler angle
$\mathbf{r}\in\mathbb{R}^{3}$, and a uniform or anisotropic scale
$\mathbf{s}\in\mathbb{R}^{3}$. To acquire the SDF of primitives in the scene
space, we inversely transform the position to the unit space and query the
unit SDF:

$\mathrm{sdf}(\mathbf{p};\theta_{s})=\mathrm{sdf}_{\text{unit}}(\mathbf{T}(\mathbf{t},\mathbf{r},\mathbf{s})^{-1}\mathbf{p},\theta_{s}^{\text{basic}}),$
(3)

where
$\theta_{s}=\\{\theta_{s}^{\text{basic}},\mathbf{t},\mathbf{r},\mathbf{s}\\}$
is the final shape parameters for primitive-based 3D strokes. In practice, we
may use a subset of translation, rotation, and scale to combine different
styles of primitive strokes. We list all the primitive-based 3D strokes in the
supplementary. More primitive can be easily added given a defined unit SDF.

####  3.2.2 Spline Curves

In addition to basic primitives, we use volumetric 3D curves with a given
radius to simulate the trace of strokes in human paintings. The curves are
defined by parametric 3D splines in the scene space:
$C:(t,\theta_{s}^{\text{curve}})\to\mathbf{x}\in\mathbb{R}^{3}$, where
$t\in[0,1]$ is the interpolation parameter, and
$\theta_{s}^{\text{curve}}\in\mathbb{R}^{d_{c}}$ is the parameter of the 3D
spline. To simulate the brushstroke effects, we define two different radii
$r_{a},r_{b}\in\mathbb{R}^{+}$ at the two endpoints of the spline curve
respectively, and interpolate the radius at any position on the curve smoothly
based on $t$ as $r(t;r_{a},r_{b})=r_{a}(1-t)+r_{b}t$.

We utilize common polynomial splines, including quadratic and cubic Bézier and
Catmull-Rom splines, to define 3D curves. To compute the SDF of a 3D curve, we
need to solve for the value of $t$ corresponding to the nearest point on the
spline curve for any point $\mathbf{p}\in\mathbb{R}^{3}$ in space. While there
exists an analytic solution for quadratic splines, it is difficult to directly
solve for the nearest point’s $t$ value for cubic or more complex splines in a
differentiable way.

Therefore, we use a general approximation method to compute the SDF: we
uniformly sample $K+1$ positions on the curve to get $K$ line segments, and
calculate the distance from these segments to the query point respectively,
thereby obtaining an approximate nearest point. This algorithm can be applied
to any parametric curve spline and allows simple control over the trade-off
between computational complexity and accuracy. With the found $t$ value of the
nearest point on the 3D spline to query position $\mathbf{p}$ as $t^{*}$, the
SDF of the 3D curve is defined as

$\mathrm{sdf}(\mathbf{p};\theta_{s})=\|\mathbf{p}-C(t^{*},\theta_{s}^{\text{curve}})\|_{2}-r(t^{*};r_{a},r_{b})$
(4)

where $\theta_{s}=\\{\theta_{s}^{\text{curve}},r_{a},r_{b}\\}$ is the final
shape parameters for spline-based 3D strokes. We leave the details of the
nearest point finding algorithm in the supplementary.

###  3.3 Differentiable Rendering of 3D Strokes

With the SDF of 3D strokes defined in Sec. 3.2, we now convert it into the
density field $\sigma(\mathbf{x})$ and color field
$\mathbf{c}(\mathbf{x},\mathbf{d})$ for differentiable rendering.

##### Deriving Density Field from SDF.

Theoretically, we consider the inner region of a 3D stroke with an SDF value
less than or equal to zero. Therefore, we can define an indicator region
function $\alpha(\mathbf{x})\in[0,1]$ based on whether a point is inside the
SDF:

$\alpha(\mathbf{x})=\left\\{\begin{aligned} &1,\;\mathrm{sdf}(\mathbf{x})\leq
0\\\ &0,\;otherwise\end{aligned}\right.$ (5)

Assuming uniform density inside one stroke, we define the density field of
each 3D stroke as

$\sigma(\mathbf{x})=\theta_{\sigma}\alpha(\mathbf{x})$ (6)

where $\theta_{\sigma}\in\mathbb{R}^{+}$ is the density parameter of a 3D
stroke. However, due to the discontinuity of the indicator function, the
gradient of the density field w.r.t. the shape parameters of the SDF is zero,
thus hindering the use of losses to optimize the 3D stroke. To render the 3D
stroke in a differentiable way, we aim to derive a smooth approximation of the
region function to gain non-zero gradients regarding the shape parameters. We
use the Laplace cumulative distribution function (CDF) to approximate the
discrete region function, similar to VolSDF [29]:

$\alpha(\mathbf{x})=\left\\{\begin{aligned}
&1-\frac{1}{2}\exp(\mathrm{sdf}(\mathbf{x})/\delta),\;\mathrm{sdf}(\mathbf{x})\leq
0\\\
&\frac{1}{2}\exp(-\mathrm{sdf}(\mathbf{x})/\delta),\;otherwise\end{aligned}\right.$
(7)

where $\delta$ controls the width of the smooth transitional interval, as
shown in Fig. 3. This definition ensures that all points in space have non-
zero derivatives w.r.t. the shape parameters, although the gradient would be
too small to achieve meaningful optimization excluding the near region besides
the zero level-set boundary of the SDF. This approximation approaches the
discrete indicator function when $\delta\to 0$, while larger $\delta$ allows
smooth gradients flow in a larger area near the boundary, at the cost of
making the shape boundary more blurred.

Figure 3: Differential region function by approximating step function with the
CDF of Laplace distribution. The parameter $\delta$ controls the width of the
transitional area with respect to the SDF of 3D strokes, where a larger
$\delta$ leads to smoother gradient at the cost of making the shape boundary
blurrier.

##### Adaptive Choice of $\delta$ based on Cone Tracing.

Since the pixel corresponding to each ray has a specific planner size during
rendering, the pinhole camera model casts a frustum with non-zero volume
instead of a ray. The frustum far away from the camera will cover a large
region of the scene, while the near frustum will only affect a small region.
Using a uniform $\delta$ across the entire scene to render the 3D strokes will
make the near objects over-blurry, and far objects lack enough gradients to
optimize the stroke shape. To solve this issue, we adopt a similar cone
tracing strategy as in MipNeRF [1], and adjust the $\delta$ value of the
region function adaptively based on the cone size. For simplicity, we use
isotropic spheres instead of anisotropic multivariate Gaussians to approximate
the volume of the cone. Specifically, assuming the radius of each pixel is
$\dot{r}$, we calculate the radius of a sphere at sample points as
$r=\dot{r}tf^{-1}$, where $t$ is the un-normalized ray distance, and $f$ is
the focal length. We then compute the adaptive $\delta$ as
$\delta=(k_{\delta}r)^{-1}$, where $k_{\delta}$ is a hyper-parameter
controlling how much to “dilate” the boundary. Moreover, as the scale of
$\delta$ is defined in scene space, for basic primitives with a
transformation, we further adjust $\delta$ based on the scaling factor of the
inverse transform.

##### Color Field.

We mainly use constant color for each stroke, which is irrelevant to the
position and view direction and can be expressed as
$\mathbf{c}(\mathbf{x},\mathbf{d})=\theta_{\mathbf{c}}$, where
$\theta_{\mathbf{c}}\in\mathbb{R}^{3}$ is the RGB color parameter. More
diverse visual effects can be achieved by spatially varying color fields and
joint modification of the density and color. The color field can be readily
expanded to support view-dependent effects by replacing the RGB color with
spherical harmonics, however, we have found that the view dependency does not
contribute much to the visual effect of stroke-based scenes in practice.

###  3.4 Composition of 3D Strokes

Given the density and color field of individual 3D strokes, we now compose
them into a stroke field for rendering the full 3D scene. In the procedure of
image painting, different strokes are painted onto the canvas sequentially,
and new strokes are overlayed on the old strokes. This behavior is simulated
using the alpha composition. Similarly, we use an “overlay” method to compose
multiple 3D strokes.

Suppose we have the density fields $\\{\sigma_{i}(\mathbf{x})\\}$ and color
fields $\\{\mathbf{c}_{i}(\mathbf{x})\\}$ of $N$ strokes, where
$i\in\\{1,2,\cdots,N\\}$. We use the region function defined in Sec. 3.3 as
the blending weight for overlay. Given $\\{\alpha_{i}(\mathbf{x})\\}$ as the
region function of each 3D stroke, the density and color of the final stroke
field are expressed as

$\begin{split}\sigma(\mathbf{x})&=\sum_{i}^{n}\sigma_{i}(\mathbf{x})\alpha_{i}(\mathbf{x})T_{i}(\mathbf{x}),\\\
\mathbf{c}(\mathbf{x})&=\frac{\sum_{i}^{n}\mathbf{c}_{i}(\mathbf{x})\alpha_{i}(\mathbf{x})T_{i}(\mathbf{x})}{1-T_{n}(\mathbf{x})},\end{split}$
(8)

where

$T_{i}(\mathbf{x})=\prod_{j}^{i}(1-\alpha_{j}(\mathbf{x})).$ (9)

Note that we normalize $\mathbf{c}(\mathbf{x})$ using the total weight of all
region functions, making the final color unaffected by the number of strokes.

The “overlay” composition is sensitive to the painting order of multiple
strokes, i.e., the 3D strokes painted later will have larger blending weights
than the previously painted strokes. Normally this is the behavior we desire,
however, we sometimes would like the final outcome irrelevant to the painting
order. The simplest way is choosing the 3D stroke with the maximum
$\alpha(\mathbf{x})$ as:
$\sigma(\mathbf{x})=\sigma_{i}(\mathbf{x}),\;\mathbf{c}(\mathbf{x})=\mathbf{c}_{i}(\mathbf{x})$,
where $i=\arg\max_{j}\\{\alpha_{j}(\mathbf{x})\\}$. This “max” composition
will only attribute the gradient of the density field to the nearest 3D
stroke. A smoother way is to use Softmax for computing weights as:
$\sigma(\mathbf{x})=\sigma_{i}(\mathbf{x})\omega_{i}(\mathbf{x}),\;\mathbf{c}(\mathbf{x})=\mathbf{c}_{i}(\mathbf{x})\omega_{i}(\mathbf{x})$,
where
$\omega_{i}(\mathbf{x})=\frac{\exp(\alpha_{i}(\mathbf{x})/\tau)}{\sum_{j}^{n}\exp(\alpha_{j}(\mathbf{x})/\tau)}$,
and $\tau$ is a hyper-parameter that controls the blending smoothness. We
analyze the effect of different compositions in Sec. 4.3.

Figure 4: Novel view synthesis results of various 3D strokes on the synthetic
scenes. Each scene contains 500 strokes. Our vectorized stroke representation
is able to recover the 3D scenes with high fidelity while maintaining a strong
geometric style defined by the strokes.

###  3.5 Training Strategy

Now we have completed our rendering framework for 3D strokes. We can render 2D
images with a given camera pose from parameters of $N$ strokes
$\\{\theta_{s},\theta_{c},\theta_{\sigma}\\}_{i},\;i\in\\{1,2,\cdots,N\\}$,
and use gradient descent to optimize these 3D stroke parameters. However, in
the previous analysis, we mentioned that the magnitude of gradients vanishes
for faraway points from the existing boundary of strokes, thus leaving shape
parameters unable to gain enough optimization. To this end, we propose a
training scheme for better stroke initialization and update.

##### Error Field.

When adding 3D strokes to an existing scene, we wish to initialize the new
stroke at the location where it is most needed. Therefore, during the training
of the stroke field, we additionally train an error field
$e(\mathbf{x})\in\mathbb{R}^{+}$ to describe the “error rate” at each position
of the reconstructed scene, with the region of highest error considered to be
the locations most in need of additional 3D strokes. The error value
$E(\mathbf{r})$ along each ray $\mathbf{r}$ is obtained by the modified volume
rendering formula:

$E(\mathbf{r})=\int_{t_{n}}^{t_{f}}e(\mathbf{r}(t))\exp(\int_{t_{n}}^{t}e(\mathbf{r}(s))ds)dt.$
(10)

We define the error as the squared difference in color between the stroke
field rendering result and the ground truth. Moreover, we prefer the error
field to conservatively estimate errors, therefore we use a coefficient $k>1$
to amplify the loss where the error is underestimated. Assuming the rendering
color of stroke field is $C(\mathbf{r})$, and the GT color is
$C_{\text{gt}}(\mathbf{r})$, then the training loss for the error field is

$\mathcal{L}_{err}=|d|k^{\max(-sgn(d),0)},\;d=E(\mathbf{r})-\|C(\mathbf{r})-C_{\text{gt}}(\mathbf{r})\|_{2},$
(11)

where $sgn(\cdot)$ is the sign function. At the same time, for views that are
less trained, we apply a regularization loss to the error values. For $N$
total sampling points $\\{\mathbf{x}_{i}\\}$ in a batch, the regularization
loss is $\mathcal{L}_{err\,reg}=\sum_{i}^{N}e(\mathbf{x}_{i})$.

##### Stroke Initialization and Update.

When training a scene with $N$ strokes, we randomly initialize
$N_{\text{start}}$ strokes and then incrementally add the remaining strokes
step by step to the scene. We also decay the initial size of strokes based on
the current number of strokes, i.e., we progressively use smaller strokes as
the number of existing strokes increases. When initializing the $i$-th stroke,
we uniformly sample $M$ coordinates $\\{\mathbf{p}_{i}\\}$ within the bounding
box of the scene, calculate the error values $\\{e_{i}\\}$ in the error field,
and select the sampling point with the highest error as the initial position
for the 3D stroke. For primitives, we set their translation vector to this
initial position. For spline curves, we set their control points to this
initial position plus a random offset vector sampled within the current stroke
size volume. This initialization scheme, guided by the error field,
effectively mitigates optimization issues caused by flat gradients.
Additionally, we reinitialize strokes with near-zero density by resampling the
current error field. This simple technique recycles some of the brush strokes
that might have been trapped in local minima due to poor initialization.

##### Training Objective.

Similar to NeRF, we use multi-view images as the primary supervision. We
randomly sample rays from images in the training set and optimize the error
between the rendering result of each ray and the ground truth. We use the
Charbonnier distance for color loss:

$\mathcal{L}_{color}=\sqrt{\|C(\mathbf{r})-C_{gt}(\mathbf{r})\|_{2}^{2}+\epsilon}$
(12)

For scenes with a ground truth mask, we also impose mask supervision based on
the ray’s opacity
$O(\mathbf{r})=\exp(\int_{t_{n}}^{t_{f}}\sigma(\mathbf{r}(t))dt)$ and the
ground truth $M(\mathbf{r})\in\\{0,1\\}$,

$\mathcal{L}_{mask}=\sqrt{\|O(\mathbf{r})-M_{gt}(\mathbf{r})\|_{2}^{2}+\epsilon}$
(13)

Additionally, we regularize the density parameter of $N$ strokes as
$\mathcal{L}_{den\,reg}=\sum_{i}^{N}|\theta_{\sigma}^{i}|$. Combining the
above losses, our total loss is given by:

$\begin{split}\mathcal{L}&=\lambda_{color}\mathcal{L}_{color}+\lambda_{mask}\mathcal{L}_{mask}\\\
&+\lambda_{den\,reg}\mathcal{L}_{den\,reg}+\lambda_{err}\mathcal{L}_{err}+\lambda_{err\,reg}\mathcal{L}_{err\,reg}\end{split}$
(14)

Figure 5: Results of different types of 3D strokes on the face-forwarding
scenes. Each scene contains 1000 strokes. Figure 6: Results with different
stroke numbers. By changing the total number of strokes, our stroke-based
representation can capture both the abstraction and fine-grained details of
the original scenes.

##  4 Experiments

We test our method on both synthetic objects from Blender [18] and face-
forwarding scenes from LLFF [17]. Each dataset comprises dozens to hundreds of
multi-view images with given camera poses. We show our reconstruction results
based on various types of vectorized 3D strokes, compare the stylization
results of our method with other image-to-painting methods, and conduct an
ablation study.

###  4.1 Stroke-based scene reconstruction

Fig. 4 shows a sample of the stylized 3D synthetic scenes reconstructed by our
method with several types of 3D strokes. Our method is able to synthesize 3D
scenes based on vectorized 3D strokes, providing strong geometry and
appearance stylization while maintaining the abstracted shapes and colors of
the original scenes. Fig. 5 compares reconstruction results of different types
of 3D strokes on the same scene. Cubic Bézier and Ellipsoid generally have the
best reconstruction fidelity compared to other types of strokes. We list the
quantitative metrics of different strokes in the supplementary material. Fig.
6 demonstrates the abstract-to-detail painting results by varying the number
of total 3D strokes that are used to reconstruct the scene. We can observe
that the lower number of strokes approximates the original scene with a strong
shape-related geometric style, while the higher number of strokes leads to
reconstruction results with higher fidelity.

Figure 7: Comparison with 2D image-to-painting methods.

###  4.2 Comparison with other methods

Our method lies in a different category of stylization compared to the common
NeRF-based stylization works, which usually take a reference image as input,
and transfer the textural style onto the appearance of a trained radiance
field. Our method does not require extra reference images, instead, the style
of our method is intrinsic and contained in the specification of various 3D
strokes.

We compare our method with 2D image-to-painting methods. Fig. 7 shows the
vectorized reconstruction results of our method and other 2D vectorized image
representation methods, including diffvg [13] and Stylized Neural Painting
[32]. We use multi-view images to train our 3D stroke representation, while
for 2D methods we select images of the specific viewpoints and apply the
vectorization tracing. It can be observed that all these vectorized methods
are able to reconstruct original images with obvious stroke styles.
Nevertheless, our method can recover more details and strictly maintain the
multi-view consistency across different views, which is guaranteed by its 3D
representation. The results are best viewed in our supplementary video.

###  4.3 Ablation studies

##### Use of adaptive $\delta$ in region function.

We investigate the effects of employing an adaptive $\delta$ in differentiable
rendering. The results in Tab. 2 show that omitting adaptive $\delta$ leads to
reduced performance in novel view synthesis. Additionally, Fig. 8 demonstrates
that constant $\delta$ values, whether low or high, result in sub-optimal or
overly blurry scenes.

##### Use of error field.

In Sec. 3.3, we note that loss gradients are concentrated near the 3D stroke
boundary, creating many local minima that heavily influence optimization based
on shape initialization. Tab. 2 demonstrates that using the error field for
stroke shape initialization significantly improves scene reconstruction
fidelity.

##### Choice of composition function.

We compare different composition approaches in Sec. 3.4. We can observe in
Fig. 9 that the ‘overlay’ composition leads to the best reconstruction
quality, while the ‘softmax’ composition also reconstructs the scene well. The
latter can be chosen if order invariant painting is required.

Method | PSNR$\uparrow$ | SSIM$\uparrow$ | LPIPS$\downarrow$
---|---|---|---
Ours (constant low $\delta$) | $21.30$ | $0.793$ | $0.177$
Ours (constant high $\delta$) | $21.86$ | $0.796$ | $0.235$
Ours (w/o adaptive $\delta$) | $22.50$ | $0.815$ | $0.179$
Ours (w/o error field) | $21.79$ | $0.811$ | $0.192$
Ours (full) | $\mathbf{23.13}$ | $\mathbf{0.825}$ | $\mathbf{0.171}$

Table 2: Novel view synthesis metrics of ablation studies. Figure 8:
Comparison of ablation study results. Figure 9: Comparison of stroke
composition methods.

###  4.4 Applications

We explore various applications based on our vectorized 3D scene
representations, including color transfer and text-driven zero-shot scene
drawing.

##### Color Transfer.

As the 3D stroke representation has separate shape and appearance parameters,
we can fix the geometry and only fine-tune the color parameter to achieve
color transfer effects. We adopt the style loss in common style transfer
works, which matches the gram matrix of the feature map outputs from the 4-th
and 9-th layers of a pre-trained VGG network. The reference style images and
color transfer results are demonstrated in Fig. 10.

Figure 10: Color transfer results. The left column is the input and style
target images respectively.

##### Text-driven Scene Drawing.

We also explore using the vectorized 3D stroke representation to achieve scene
creation tasks under a text-guided zero-shot generation framework. We use CLIP
[23], a vision-language model that embeds the 2D images and text prompts into
the same embedding space. Following the generative setup in DreamField [11],
we sample camera poses following a circular path around the scene origin and
render a large image patch, and then we optimize the distance between CLIP
embeddings of the synthesized image and the text prompt. The generated 3D
drawings of different objects are shown in Fig. 11. We leave more details of
the training setup in the supplementary material.

Figure 11: Text-driven 3D drawings of common objects with 300 strokes of cubic
Bézier curve.

##  5 Discussion

We present a novel method to stylize a 3D scene from multi-view 2D images.
Different from NeRF-based representations, our method represents the scene as
vectorized 3D strokes, mimicking human painting during scene reconstruction
process. We demonstrate that this stroke-based representation can successfully
stylize 3D scenes with large geometry and appearance transformations, which
was not achieved with previous NeRF stylization approaches.

##### Limitations and future works.

Our method uses a stroke setting that demands manual effort to design stroke
shapes and appearances. The 3D strokes can be further learned with a
generative framework to create a variety of stroke types, like ink and oil
brushes, with generated and more detailed SDFs. Moreover, our method may
require numerous strokes to represent very complex scenes, partly due to the
existence of many local minima during the optimization. Incorporating a
globally aware loss, like the optimal transport loss in [32] into 3D space,
may enhance the convergence efficiency of our method, which we leave as future
work.

## References

  * Barron et al. [2021] Jonathan T Barron, Ben Mildenhall, Matthew Tancik, Peter Hedman, Ricardo Martin-Brualla, and Pratul P Srinivasan.  Mip-nerf: A multiscale representation for anti-aliasing neural radiance fields.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 5855–5864, 2021.
  * Barron et al. [2023] Jonathan T Barron, Ben Mildenhall, Dor Verbin, Pratul P Srinivasan, and Peter Hedman.  Zip-nerf: Anti-aliased grid-based neural radiance fields.  _arXiv preprint arXiv:2304.06706_ , 2023.
  * Chiang et al. [2022] Pei-Ze Chiang, Meng-Shiun Tsai, Hung-Yu Tseng, Wei-Sheng Lai, and Wei-Chen Chiu.  Stylizing 3d scene via implicit representation and hypernetwork.  In _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , pages 1475–1484, 2022.
  * Chittenden [2018] Tara Chittenden.  Tilt brush painting: Chronotopic adventures in a physical-virtual threshold.  _Journal of contemporary painting_ , 4(2):381–403, 2018.
  * Fan et al. [2022] Zhiwen Fan, Yifan Jiang, Peihao Wang, Xinyu Gong, Dejia Xu, and Zhangyang Wang.  Unified implicit neural stylization.  In _European Conference on Computer Vision_ , pages 636–654. Springer, 2022.
  * Gatys et al. [2015] Leon A Gatys, Alexander S Ecker, and Matthias Bethge.  A neural algorithm of artistic style.  _arXiv preprint arXiv:1508.06576_ , 2015.
  * Gatys et al. [2016] Leon A Gatys, Alexander S Ecker, and Matthias Bethge.  Image style transfer using convolutional neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2414–2423, 2016.
  * Haeberli [1990] Paul Haeberli.  Paint by numbers: Abstract image representations.  In _Proceedings of the 17th annual conference on Computer graphics and interactive techniques_ , pages 207–214, 1990.
  * Hertzmann [1998] Aaron Hertzmann.  Painterly rendering with curved brush strokes of multiple sizes.  In _Proceedings of the 25th annual conference on Computer graphics and interactive techniques_ , pages 453–460, 1998.
  * Huang et al. [2022] Yi-Hua Huang, Yue He, Yu-Jie Yuan, Yu-Kun Lai, and Lin Gao.  Stylizednerf: Consistent 3d scene stylization as stylized nerf via 2d-3d mutual learning.  In _Computer Vision and Pattern Recognition (CVPR)_ , 2022.
  * Jain et al. [2022] Ajay Jain, Ben Mildenhall, Jonathan T. Barron, Pieter Abbeel, and Ben Poole.  Zero-shot text-guided object generation with dream fields.  2022\.
  * Johnson et al. [2016] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  In _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part II 14_ , pages 694–711. Springer, 2016.
  * Li et al. [2020] Tzu-Mao Li, Michal Lukáč, Michaël Gharbi, and Jonathan Ragan-Kelley.  Differentiable vector graphics rasterization for editing and learning.  _ACM Transactions on Graphics (TOG)_ , 39(6):1–15, 2020.
  * Liu et al. [2023] Kunhao Liu, Fangneng Zhan, Yiwen Chen, Jiahui Zhang, Yingchen Yu, Abdulmotaleb El Saddik, Shijian Lu, and Eric P Xing.  Stylerf: Zero-shot 3d style transfer of neural radiance fields.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8338–8348, 2023.
  * Liu et al. [2021] Songhua Liu, Tianwei Lin, Dongliang He, Fu Li, Ruifeng Deng, Xin Li, Errui Ding, and Hao Wang.  Paint transformer: Feed forward neural painting with stroke prediction.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 2021.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Mildenhall et al. [2019] Ben Mildenhall, Pratul P Srinivasan, Rodrigo Ortiz-Cayon, Nima Khademi Kalantari, Ravi Ramamoorthi, Ren Ng, and Abhishek Kar.  Local light field fusion: Practical view synthesis with prescriptive sampling guidelines.  _ACM Transactions on Graphics (TOG)_ , 38(4):1–14, 2019.
  * Mildenhall et al. [2021] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Communications of the ACM_ , 65(1):99–106, 2021.
  * Müller et al. [2022] Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller.  Instant neural graphics primitives with a multiresolution hash encoding.  _ACM Transactions on Graphics (ToG)_ , 41(4):1–15, 2022.
  * Nakano [2019] Reiichiro Nakano.  Neural painters: A learned differentiable constraint for generating brushstroke paintings.  _arXiv preprint arXiv:1904.08410_ , 2019.
  * Nguyen-Phuoc et al. [2022] Thu Nguyen-Phuoc, Feng Liu, and Lei Xiao.  Snerf: Stylized neural implicit representations for 3d scenes.  _ACM Trans. Graph._ , 41(4), 2022.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Radford et al. [2021] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al.  Learning transferable visual models from natural language supervision.  In _International conference on machine learning_ , pages 8748–8763. PMLR, 2021.
  * Risser et al. [2017] Eric Risser, Pierre Wilmot, and Connelly Barnes.  Stable and controllable neural texture synthesis and style transfer using histogram losses.  _arXiv preprint arXiv:1701.08893_ , 2017.
  * Seo et al. [2009] SangHyun Seo, JinWan Park, and KyungHyun Yoon.  A painterly rendering based on stroke profile and database.  In _CAe_ , pages 9–16, 2009.
  * Sheffield [1985] Charles Sheffield.  Selecting band combinations from multi spectral data.  _Photogrammetric Engineering and Remote Sensing_ , 58(6):681–687, 1985.
  * Shen et al. [2018] Falong Shen, Shuicheng Yan, and Gang Zeng.  Neural style transfer via meta networks.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 8061–8069, 2018.
  * Shiraishi and Yamaguchi [2000] Michio Shiraishi and Yasushi Yamaguchi.  An algorithm for automatic painterly rendering based on local source image approximation.  In _Proceedings of the 1st international symposium on Non-photorealistic animation and rendering_ , pages 53–58, 2000.
  * Yariv et al. [2021] Lior Yariv, Jiatao Gu, Yoni Kasten, and Yaron Lipman.  Volume rendering of neural implicit surfaces.  _Advances in Neural Information Processing Systems_ , 34:4805–4815, 2021.
  * Zhang et al. [2022] Kai Zhang, Nick Kolkin, Sai Bi, Fujun Luan, Zexiang Xu, Eli Shechtman, and Noah Snavely.  Arf: Artistic radiance fields.  In _European Conference on Computer Vision_ , pages 717–733. Springer, 2022.
  * Zhao and Zhu [2011] Mingtian Zhao and Song-Chun Zhu.  Customizing painterly rendering styles using stroke processes.  In _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on non-photorealistic animation and rendering_ , pages 137–146, 2011.
  * Zou et al. [2021] Zhengxia Zou, Tianyang Shi, Shuang Qiu, Yi Yuan, and Zhenwei Shi.  Stylized neural painting.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 15689–15698, 2021.

Supplementary Material

This supplementary document provides additional details of the 3D strokes in
Sec. 6, implementations in Sec.7, additional comparisons of various 3D strokes
in Sec. 8, and application training setups in Sec. 9. Please also watch our
accompanying video for an animated visualization of stylization results.

##  6 Details of 3D Strokes

###  6.1 Transformation of Basic Primitives

In Sec. 3.2.1 of the main paper, we use a transformation matrix to map the
coordinates in the shared scene space into the canonical space of each unit
signed distance field. Here we provide the construction details of the
transformation matrix. Given a translation vector
$\mathbf{t}=(t_{x},t_{y},t_{z})$, an Euler angle rotation vector
$\mathbf{r}=(r_{x},r_{y},r_{z})$, and a scale vector
$\mathbf{s}=(s_{x},s_{y},s_{z})$, we first construct the matrices for each
term respectively, then combine them in the order of scale, rotation, and
translation to get the final transformation matrix $M$:

$\begin{split}T&=\begin{bmatrix}1&0&0&t_{x}\\\ 0&1&0&t_{y}\\\ 0&0&1&t_{z}\\\
0&0&0&1\end{bmatrix},\\\ R_{x}&=\begin{bmatrix}1&0&0&0\\\ 0&\cos r_{x}&-\sin
r_{x}&0\\\ 0&\sin r_{x}&\cos r_{x}&0\\\ 0&0&0&1\end{bmatrix},\\\
R_{y}&=\begin{bmatrix}\cos r_{y}&0&\sin r_{y}&0\\\ 0&1&0&0\\\ -\sin
r_{y}&0&\cos r_{y}&0\\\ 0&0&0&1\end{bmatrix},\\\ R_{z}&=\begin{bmatrix}\cos
r_{z}&-\sin r_{z}&0&0\\\ \sin r_{z}&\cos r_{z}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{bmatrix},\\\ S&=\begin{bmatrix}s_{x}&0&0&0\\\ 0&s_{y}&0&0\\\
0&0&s_{z}&0\\\ 0&0&0&1\end{bmatrix},\\\ M&=TR_{z}R_{y}R_{x}S\end{split}$ (15)

In cases where composited primitives do not utilize the full set of
transformation components—translation, rotation, and scale—we omit the
respective term in the $M$ formula. For uniform scaling, represented by a
scalar scaling factor $s$, we set $s_{x}$, $s_{y}$, and $s_{z}$ all equal to
$s$.

###  6.2 Complete List of 3D Strokes

Stroke Name | Base SDF | Translation | Rotation | Uniform Scale | Anisotropic Scale
---|---|---|---|---|---
Sphere | Unit Sphere | ✓ | ✗ | ✓ | ✗
Ellipsoid | Unit Sphere | ✓ | ✓ | ✗ | ✓
Axis-aligned Cube | Unit Cube | ✓ | ✗ | ✓ | ✗
Oriented Cube | Unit Cube | ✓ | ✓ | ✓ | ✗
Axis-aligned Box | Unit Cube | ✓ | ✗ | ✗ | ✓
Oriented Box | Unit Cube | ✓ | ✓ | ✗ | ✓
Round Cube | Unit RoundCube | ✓ | ✓ | ✓ | ✗
Round Box | Unit RoundCube | ✓ | ✓ | ✗ | ✓
Line | Unit Line | ✓ | ✓ | ✓ | ✗
Triprism | Unit Triprism | ✓ | ✓ | ✓ | ✗
Octahedron | Unit Octahedron | ✓ | ✓ | ✓ | ✗
Tetrahedron | Unit Tetrahedron | ✓ | ✓ | ✓ | ✗
Table 3: The basic shapes and transformations used in all 3D strokes.

Utilizing various combinations of basic geometric shapes in unit space, along
with transformations including translation, rotation, and scale, enables the
creation of a diverse palette of 3D strokes. These strokes exhibit distinct
geometric and aesthetic stylization. The full assortment of these 3D strokes
is detailed in Tab. 3.

###  6.3 Spline Curves

####  6.3.1 Polynomial splines

In Sec. 3.2.2 of the main paper, we use three types of different polynomial
curves that are commonly used in computer graphics. Specifically, we use the
quadratic Bézier, cubic Bézier, and Catmull Rom spline, respectively. We
provide the concrete definition of these curves below. All points defined here
are vectors in the 3D scene space.

##### Quadratic Bézier Spline.

A quadratic Bézier spline is defined by three control points, the start point
$\mathbf{P}_{0}$, the end point $\mathbf{P}_{2}$, and middle point
$\mathbf{P}_{1}$ that is also tangent to the $\mathbf{P}_{0}$ and
$\mathbf{P}_{2}$. Note that the curve only blends toward but does not pass the
middle point $\mathbf{P}_{1}$. The parametric form is given by:

$\mathbf{C}(t;\mathbf{P}_{0},\mathbf{P}_{1},\mathbf{P}_{2})=(1-t)^{2}\mathbf{P}_{0}+2(1-t)t\mathbf{P}_{1}+t^{2}\mathbf{P}_{2}$
(16)

##### Cubic Bézier Spline.

A cubic Bézier spline introduces an additional control point compared with the
quadratic spline. This spline is defined by four points: the start point
$\mathbf{P}_{0}$, two control points $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$,
and the end point $\mathbf{P}_{3}$. The curve starts at $\mathbf{P}_{0}$ and
ends at $\mathbf{P}_{3}$, with $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$
influencing its shape. The parametric equation of a cubic Bézier spline is:

$\begin{split}\mathbf{C}(t;\mathbf{P}_{0},\mathbf{P}_{1},\mathbf{P}_{2},\mathbf{P}_{3})&=(1-t)^{3}\mathbf{P}_{0}+3(1-t)^{2}t\mathbf{P}_{1}\\\
&+3(1-t)t^{2}\mathbf{P}_{2}+t^{3}\mathbf{P}_{3}\end{split}$ (17)

##### Catmull Rom Spline.

The Catmull Rom spline is another form of cubic spline, notable for its
ability to interpolate its control points. This spline is defined by a series
of points, with the curve passing through each of these points except the
first and last. One feature of the Catmull Rom spline is that the tangent at
each point is determined by the line connecting the previous and next points,
ensuring a smooth transition. In our implementation, we specifically use the
centripetal Catmull-Rom spline, a variant of the standard Catmull-Rom spline.
This type of spline is particularly advantageous for avoiding the issue of
self-intersecting loops in the curve, which are common in the uniform and
chordal Catmull-Rom splines. The parametric form of the Catmull Rom spline,
for a segment between $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$, with
$\mathbf{P}_{0}$ and $\mathbf{P}_{3}$ influencing the shape, is defined as:

$\begin{split}\mathbf{C}(t;\mathbf{P}_{0},\mathbf{P}_{1},\mathbf{P}_{2},\mathbf{P}_{3})&=\frac{1}{2}[(2\mathbf{P}_{1})\\\
&+(-\mathbf{P}_{0}+\mathbf{P}_{2})t\\\
&+(2\mathbf{P}_{0}-5\mathbf{P}_{1}+4\mathbf{P}_{2}-\mathbf{P}_{3})t^{2}\\\
&+(-\mathbf{P}_{0}+3\mathbf{P}_{1}-3\mathbf{P}_{2}+\mathbf{P}_{3})t^{3}]\end{split}$
(18)

####  6.3.2 Nearest Point Finding

As mentioned in Sec. 3.2.2 of the main paper, we need to locate the nearest
point on the spline in order to compute the SDF of the spline curve. This
involves solving the following equation that finds the $t$ value with the
minimized distance to the query point $\mathbf{p}\in\mathbb{R}^{3}$ in a
differential way:

$t^{*}=\underset{t}{\arg\min}\|\mathbf{C}(t,\theta_{s}^{\text{curve}})-\mathbf{p}\|_{2},\;s.t.\;0\leq
t\leq 1\;,$ (19)

where $\theta_{s}^{\text{curve}}$ denotes the parameters of the spline curve.
An analytical solution might exist for some specific formulations of the
splines. However, for more versatility, we use a general approximation
solution that can be adapted to any parametric spline curve in the main paper.
$K+1$ samples are uniformly selected on the curve to form $K$ line segments,
and the distance of each line segment to the query point is calculated to find
$t^{*}$. Assuming the line segment is given as
$L(t;\mathbf{A},\mathbf{B})=(1-t)\mathbf{A}+t\mathbf{B},\;t\in[0,1]$, where
$\mathbf{A}$ and $\mathbf{B}$ are the two endpoints of the line segment, the
distance from point $\mathbf{p}$ to this line segment $L$ can be calculated
as:

$d_{L}(\mathbf{p},\mathbf{A},\mathbf{B})=\|\mathbf{p}-L(t^{\prime};\mathbf{A},\mathbf{B})\|_{2}$
(20)

where $t^{\prime}$ is the $t$ value that gives the nearest point on the line
segment:

$t^{\prime}(\mathbf{p},\mathbf{A},\mathbf{B})=\min(\max(\frac{(\mathbf{p}-\mathbf{A})\cdot(\mathbf{B}-\mathbf{A})}{(\mathbf{B}-\mathbf{A})\cdot(\mathbf{B}-\mathbf{A})},0),1)$
(21)

We then compute $t^{*}$ using Algorithm 1. The computation complexity can be
easily controlled by adjusting the number of $K$, where a higher $K$ leads to
a more accurate approximation but at the cost of higher computation.

Input: $\mathbf{C}$: parametric spline function defining the coordinates at
distance $t\in[0,1]$.

$\theta_{s}^{curve}$: parameters of the spline curve.

$K$: number of line segments used.

$\mathbf{p}$: coordinate of the query point.

Result: $t^{*}$: the $t$ value of the nearest point on the spline curve.

$t^{*}\leftarrow 0$;

$t_{\text{start}}\leftarrow 0$;

$d_{\text{min}}\leftarrow\infty$;

$\mathbf{A}\leftarrow\mathbf{C}(t_{\text{start}},\theta_{s}^{curve})$;

for  _$i\leftarrow 1$ to $K$_ do

$t_{\text{end}}\leftarrow i/K$;

$\mathbf{B}\leftarrow\mathbf{C}(t_{\text{end}},\theta_{s}^{curve})$;

if  _$d_{L}(\mathbf{p},\mathbf{A},\mathbf{B}) <d_{\text{min}}$_ then

$d_{\text{min}}=d_{L}(\mathbf{p},\mathbf{A},\mathbf{B})$;

$t^{*}=t_{\text{start}}+(t_{\text{end}}-t_{\text{start}})\cdot
t^{\prime}(\mathbf{p},\mathbf{A},\mathbf{B})$;

end if

$t_{\text{start}}\leftarrow t_{\text{end}}$;

$\mathbf{A}\leftarrow\mathbf{B}$;

end for

Algorithm 1 Find the approximated nearest point on a given spline curve by
uniformly sample $K$ line segments on the curve.

##  7 Implementation Details

We implement the stroke field using Pytorch [22] framework and implement the
strokes using fused CUDA kernels to accelerate training and reduce GPU memory
usage. We transform sampled coordinates to the canonical volume according to
the scene bounding box and use the normalized scene coordinates as inputs to
the stroke field. Like ZipNeRF [2], we use a proposal network based on hash
grid representation [19] to facilitate ray sampling. Specifically, for each
ray of a pixel, we first sample 32 points using the proposal MLP to obtain
sampling weights. We resample 32 points and compute the stroke field’s density
and color on each point, and use the same volumetric rendering formula in NeRF
to acquire the final pixel color.

We train 15k steps using 500 strokes for scenes with a single object, and
train 25k steps using 1000 strokes for face-forwarding scenes. We employ the
AdamW [16] optimizer with betas $(0.9,0.99)$, setting the learning rate to
$0.01$ and exponentially decays to $0.0003$. We start with $k_{\delta}=7$ and
gradually decay it to $k_{\delta}=1$ during training. Additionally, as brushes
are progressively added to the scene, we start with 25% of the sampling points
and gradually use all the sampling points when training reaches 80%. We set
$\lambda_{color}=1$, $\lambda_{mask}=0.02$, $\lambda_{den\,reg}=0.0001$,
$\lambda_{err}=0.1$, and $\lambda_{err\,reg}=0.001$ in our training setup.

##  8 Quantitative comparison of strokes

In Sec. 4.1 of the main paper, we conduct a qualitative comparison between the
visual effects of several selected 3D strokes. Additionally, this section
includes a comprehensive quantitative analysis of all 3D strokes, as detailed
in Tab. 4. The metrics are measured using 500 strokes on object scenes and
1000 strokes on face-forwarding scenes. As shown in the table, the ellipsoid
stroke typically demonstrates superior fidelity in scene reconstruction,
followed closely by the cubic Bézier curve.

Table 4: Averaged quantitative metrics of reconstruction results of different 3D strokes. |  PSNR$\uparrow$ |  SSIM$\uparrow$ |  LPIPS$\downarrow$
---|---|---|---
Sphere | $19.72$ | $0.591$ | $0.416$
Ellipsoid | $\mathbf{21.78}$ | $\mathbf{0.687}$ | $\mathbf{0.283}$
Cube | $19.65$ | $0.593$ | $0.389$
Axis-aligned Box | $19.99$ | $0.588$ | $0.408$
Oriented Box | $20.66$ | $0.637$ | $0.321$
Round Cube | $20.17$ | $0.623$ | $0.380$
Tetrahedron | $20.30$ | $0.625$ | $0.372$
Octahedron | $20.09$ | $0.614$ | $0.380$
Triprism | $20.68$ | $0.640$ | $0.351$
Line | $20.74$ | $0.641$ | $0.347$
Quadratic Bézier | $21.32$ | $0.676$ | $0.311$
Cubic Bézier | $21.64$ | $\mathbf{0.687}$ | $0.308$
Catmull-Rom | $21.47$ | $0.675$ | $0.324$

##  9 Training setup of applications

###  9.1 Color Transfer

In Sec. 4.4 of the main paper, we transfer the color distribution from a
reference style image to a trained stroke-based 3D scene. We adopt perceptual
style loss, which extracts the gram matrix [7] at specific layers of a pre-
trained VGG16 network [12], and computer the difference between the rendered
RGB image and the target style image. Since computing style loss requires an
image rather than individual pixels as input, we randomly render 32x32 chunks
of original images under the given camera poses. We add this style loss with a
weight $\lambda_{style}=0.25$ to the total loss and fine-tune the color
parameters of trained strokes for 5k iterations.

###  9.2 Text-driven scene drawing

In Sec. 4.4 of the main paper, we use the vision-language model CLIP [23]
(ViT-B/32) to achieve scene drawing based on a given text prompt.
Specifically, we minimize the loss between the text embedding of the given
prompt and the image embedding of the rendered RGB image. When rendering the
images, we sample camera poses in a circular path looking at the scene’s
origin with azimuth angle in $[0,360]$ and elevation angle in $[40,105]$, and
render an image chunk of size $128\times 128$. We consider the azimuth angle
starting at zero as the front view and adjust the CLIP guidance scale for
larger azimuth angles accordingly. This adjustment results in generation
outcomes that are more coherent with the viewpoint.

For text guidance, we use the template “realistic 3D rendering painting of
[OBJECT]”, where [OBJECT] is substituted with the specific object description
we aim to generate. Additionally, we discovered that incorporating a
silhouette loss alongside the RGB loss significantly enhances the quality of
the generated geometric shapes. This is achieved by generating another image
embedding from the rendered opacity and encouraging lower loss between this
silhouette image embedding with the text embedding.


# The Constraining Capability of BNS Dark Sirens Observed by the LIGO
Gravitational Wave Detector on the Hubble Constant

Xu Chao  Astronomical Society of WIT, Wuhan Institute of Technology, Wuhan
430205,China School of Chemical Engineering and Pharmacy, Wuhan Institute of
Technology, Wuhan 430205,China

Abstract

The Hubble Constant observed at high redshift and low redshift are
inconsistent, representing one of the urgent issues to be resolved in the
field of cosmology. The discovery of gravitational waves opens a new window
for addressing this problem. For instance, the GW170817 event, through the
coordinated observation of electromagnetic and gravitational wave signals,
allows for constraints to be imposed from a completely new perspective.
However, the number of gravitational wave events where both electromagnetic
and gravitational wave signals are observed simultaneously is too small,
making it difficult to enhance the precision through statistical methods. In
this paper, we use dark sirens as the subjects of study. Through the standard
gravitational wave data simulation and the analysis process, we analyze the
constraints a typical binary neutron star merger event can place on the Hubble
Constant. We simulated a random event and found that it an provide an error of
+0.04-0.05 for the Hubble Constant. By combining multiple events, this
constraint can be improved.

Keyword: Gravitational waves; Dark sirens; Hubble Constant; Binary neutron
star

##  1 Introduction

In 1929, Edwin Hubble unveiled a groundbreaking discovery, establishing
through the study of the distance-redshift relationship of galaxies that the
universe is in a state of continuous expansion[1]. He formulated the
relationship between the recession velocity of galaxies and their distance as
$v=H_{0}r$, where $H_{0}$ is known as the Hubble constant, a pivotal parameter
in modern cosmology denoting the rate of cosmic expansion.

The precise constraint of the Hubble constant can be ascertained through
various methodologies, including employing the distance ladder involving
Cepheid variables and Type Ia supernovae, among other observational
techniques[2, 3, 4, 5]. However, a notable discrepancy has surfaced concerning
the values obtained from different approaches. Current measurements of the
Hubble constant fluctuate significantly between 67.4 and 76.8
$\rm{km}\cdot\rm{s^{-1}}\cdot\rm{Mpc^{-1}}$[6], a variation that vastly
surpasses the standard uncertainty. This discrepancy, termed the “Hubble
Tension”, poses a substantial challenge to the existing cosmological paradigm,
hinting at the potential necessity for new physics beyond the standard
cosmological model or perhaps attributing to measurement errors[7]. Numerous
studies are fervently working to enhance the precision of these measurements,
aspiring to delineate a more definitive conclusion in the foreseeable
future[8, 9].

The advent of gravitational wave (GW) astronomy offers a fresh and promising
avenue for study astronomy, cosmology, and physics[10, 11]. Unlike
electromagnetic signals, gravitational waves are not obstructed by matter,
potentially relaying rich information from distant astronomical entities.
Moreover, gravitational wave observation raised a new method in measuring the
luminosity distance of cosmic sources. It will becomes a powerful tool if the
corresponding electromagnetic (EM) signal is observed. On 17 August 2017 the
Advanced LIGO and Virgo detectors observed GW170817[12] and its
electromagnetic counter part GRB 170817A[13], they determine the Hubble
constant to be
$70.0^{+12.0}_{-8.0}\rm{km}\cdot\rm{s^{-1}}\cdot\rm{Mpc^{-1}}$[14]. However,
only a fraction of these compact binary sources have EM counter part, making
it hard to enhance the precision through statistical methods. Fortunately, we
can still extract the redshift information from binary neutron stars(BNS)[15].
Such BNS gravitational wave sources without electromagnetic counterparts as
“BNS dark sirens”. To fully harness the potential of dark sirens, innovative
approaches are needed to measure redshift information. In 2012, Messenger and
Read proposed that the tidal deformations of neutron stars could contribute
additional information to the phase evolution of gravitational waveforms,
breaking the degeneracy between mass and redshift[15]. This allows for the
acquisition of redshift information from binary neutron star merger events
exclusively through gravitational wave observations, facilitating the
establishment of a relationship between luminosity distance ($D_{L}$) and
redshift ($z$), and subsequently constraining the Hubble constant.

This study predominantly focuses on the gravitational wave dark sirens
resulting from BNS mergers, endeavoring to foresee the constraints this method
can impose on the Hubble constant through predictions based on the LIGO
gravitational wave detector. The paper is structured as follows: Section 2
delineates the fundamental theories and computational methods concerning the
binary neutron star gravitational wave dark sirens; Section 3 Section 3
introduces the pertinent methods and tools for gravitational wave parameter
analysis and discusses the error analysis of the computed results by Bayesian
methods in this section; finally, Section 4 encapsulates the entire study.
Given the scarcity of gravitational wave events, this study leverages
simulated data based on the LIGO detector for analysis. To facilitate
computation, we adopt a unit system where $c=G=1$.

##  2 Theory

According to the theory of general relativity, the gravitational waves
generated by the merger of binary compact object can be described by two
independent components, $h_{+}$ and $h_{\times}$ in transverse-traceless
gauge. In the frequency domain, the response of gravitational waves can be
described by these two independent components[16],

$h(f)=F_{+}(\theta;\phi;\psi)h_{+}(f)+F_{\times}(\theta;\phi;\psi)h_{\times}(f)$
(1)

where $F_{+}$ and $F_{\times}$ are the antenna pattern functions of the
detector, ($\theta$, $\phi$) represent the angles indicating the position of
the gravitational wave source in the sky relative to the detector, and $\psi$
is the polarization angle of the gravitational wave. The components $h_{+}(f)$
and $h_{\times}(f)$ can be obtained through post-Newtonian approximations or
numerical relativity methods. A significant amount of work has been dedicated
to the research of waveform construction[17, 18, 19]. Generally, for binary
black hole mergers, $h_{+/\times}(f)$ are functions of $\vec{\theta}$ =
$\vec{\theta}(m_{i},a_{i},\theta_{i},\phi_{jl},\phi_{12},\theta_{jn},D_{L},\psi,\phi,t_{0},\alpha,\delta)$,
where $m_{i}$ is the detector-frame mass of the $i$-th object, $a_{i}$ is the
dimensionless spin magnitude of the $i$-th object, $\theta_{i}$ is the zenith
angle between the spin and orbital angular momenta for the ith object ( where
i = 1, 2), $\phi_{jl}$ is the difference between total and orbital angular
momentum azimuthal angles, $\phi_{12}$ is the difference between the azimuthal
angles of the individual spin vector projections onto the orbital plane,
$\theta_{jn}$ is the zenith angle between the total angular momentum and the
line of sight, $D_{L}$ is the luminosity distance to the source, $\psi$ is the
polarization angle of the source, $\phi$ is the binary phase at a reference
frequency, $t_{0}$ is the GPS reference time at the geocenter (typically
merger time), $\alpha$ and $\delta$ are right ascension and declination of the
source.

For the merger of binary neutron stars, an additional parameter $\lambda$ is
introduced to quantify the equation of state properties of neutron stars.
During their orbital motion, the neutron stars experience tidal deformations
due to tidal forces, which contribute additional phases to the gravitational
waves[15] (ignoring the spin of the neutron stars),

$\Psi(f)=\Psi_{\text{BBH}}(f)+\Psi_{\text{tidal}}(f)$ (2)

where $\Psi_{\text{BBH}}(f)$ is the phase term in the waveform for binary
black holes,

$\Psi_{\text{tidal}}(f)=\sum_{j=1}^{2}\frac{3\lambda_{j}(1+z)^{5}}{128\eta
M^{5}}\left[-\frac{24}{\chi_{j}}(1+\frac{11\eta}{\chi_{j}}v^{5})-\frac{5}{28\chi_{j}}(3179-919\chi_{j}-2286\chi_{j}^{2}+260\chi_{j}^{3})v^{7}\right]$
(3)

In this work, we choose the SLy model[20], which aligns well with current
observations, as the reference equation of state for neutron stars[21], and
use a linear fitting function to relate the tidal deformability $\lambda$ to
the neutron star mass $m$,

$\lambda_{j}=Bm_{j}+C$ (4)

Here, $B$ and $C$ are two tidal effect parameters; in this paper, we take
$B=-1.99$ and $C=4.46$[22]. This equation encodes redshift information into
the gravitational wave waveform. In subsequent parameter inference, we can
infer $D_{L}$ and $z$ simultaneously. Knowing both $D_{L}$ and $z$, we can
constrain the Hubble Constant using cosmological models.

##  3 Method

The gravitational wave signals detected by the LIGO (Laser Interferometer
Gravitational-Wave Observatory) detectors are typically buried within
noise[23]. To extract the signals, we often employ a method called Match-
Filter[24].In this work, we utilize the IMRPhenomD waveform to generate
signals[25, 26] and incorporate the additional phase introduced in the
previous section into the waveform as our gravitational wave signal $s(t)$.We
generate simulated noise $n(t)$ using the power spectral density (PSD) of
LIGO.The noise is then added to the simulated gravitational wave signal,
resulting in the simulated data$d(t)=n(t)+s(t)$. Standard gravitational wave
signal processing techniques are applied to process the simulated data,
followed by Bayesian analysis to extract the parameters from the gravitational
wave signal.

For computational convenience, we have chosen a binary system of
$m_{1}=m_{2}=1.4M_{\odot}$ for the gravitational wave event, which is also a
typical neutron star mass[27].Due to the fact that neutron star mergers
typically occur at the end of neutron star evolution, we assume that their
spin is zero. In this work, we consider gravitational wave sources without
electromagnetic counterparts, so we assume $\theta_{\rm jn}=0$. The luminosity
distance $D_{L}$ is choosed as $40$Mpc, which is about the distance of
GW170817[12]. We randomly select parameters such as position, polarization,
and phase to generate simulated gravitational wave signals and the redshift
used for the injection is calculated using the results of Planck18 as the
fiducial cosmological model.

In the parameter inference of gravitational waves, we employ Bayesian analysis
to extract the parameters from the gravitational wave signals[28, 29, 30]:

$p(\vec{\theta}|d)=\frac{\mathcal{L}(\vec{\theta})\pi(\vec{\theta})}{\mathcal{Z}}$
(5)

where $\pi(\theta)$ represents the prior information,
$\mathcal{L}(\theta)$represents the likelihood function, $p(\theta|d)$ is the
posterior probability density, and $\mathcal{Z}$ is a normalization constant.

In our work, for computational convenience, we fix other parameters except
luminosity distance, redshift, and merger time. We will complete the full
Bayesian analysis in future work.Our results are shown in Figure 1. From this
figure, it can be observed that the data processing and parameter inference
methods can provide good constraints on the luminosity distance
$D_{L}=38.67^{+2.85}_{-2.41}$, but they are unable to provide accurate
estimates for the redshift. We obtained a relatively large error in the
inferred redshift, which is $z=0.02^{+0.03}_{-0.02}$。This is consistent with
our expectations, as the sensitivity of LIGO detectors is currently not
sufficient to accurately obtain redshift information using this method.

Figure 1: The results of Bayesian analysis

Once the redshift and luminosity distance are known, the Hubble constant can
be calculated directly by $H_{0}=v/r=zc/D_{L}$, where $c$ represents the speed
of light, $z$ refers to the redshift of the source, and $D_{L}$ represents the
luminosity distance of the source. The result is shown in Figure 2. In this
figure, we can see the uncertainty is very small, i.e.,
$67.24^{+0.04}_{-0.05}\rm{km}\cdot\rm{s^{-1}}\cdot\rm{Mpc^{-1}}$.

Figure 2: The uncertainty of Hubble constant

##  4 Conclusion

Gravitational waves from neutron star mergers events encode absolute
luminosity distance and redshift chirp mass in their amplitude, while
additional contributions in the phase can decouple the chirp mass from the
redshift. This allows establishing a relationship between luminosity distance
( $D_{L}$ ) and redshift ($z$), thereby constraining the Hubble constant
($H_{0}$). In this paper, we studied the observational method based on
gravitational waves from binary neutron star mergers and provided predictions
for the constraining power of this method on the Hubble constant.We considered
the LIGO-Hanford (H1) and LIGO-Livingston (L1) detector network, simulated
gravitational wave signals from binary neutron star mergers using their
instrument parameters, and used the Injection-Inference method to study the
constraining power of LIGO on the Hubble Constant.

Our results indicate that, assuming the knowledge of parameters other than
luminosity distance, redshift, and merger time, the two detectors of LIGO
possess strong observational properties for the Hubble Constant. However, the
impact of other parameters on the Hubble Constant warrants further exploration
in subsequent articles.

In the future, with the construction of third-generation gravitational wave
detectors, we will be able to detect a larger number of gravitational wave
signals from binary neutron star mergers.The BNS dark sirens will still play
an important role. We expect that in the future, BNS dark sirens will find
widespread applications and provide valuable constraints on the Hubble
Constant. These independent observations, separate from electromagnetic
radiation, will complement existing observations and contribute to a more
comprehensive understanding of the universe.

## References

  * [1] E. A. Hubble.  Hubble, E.: A relation between distance and radial velocity among extra-galactic nebulae.  Proceedings of the National Academy of Sciences, 15(3):168–173, 1929.
  * [2] Jaeger T De, B E Stahl, W Zheng, A V Filippenko, A G Riess, and L Galbany.  A measurement of the Hubble constant from Type II supernovae.  Monthly Notices of the Royal Astronomical Society, 496(3):3402–3411, 2020.
  * [3] Adam G. Riess, Lucas M. Macri, Samantha L. Hoffmann, Dan Scolnic, Stefano Casertano, Alexei V. Filippenko, Brad E. Tucker, Mark J. Reid, David O. Jones, and Jeffrey M. Silverman.  A 2.4% Determination of the Local Value of the Hubble Constant.  The Astrophysical Journal, 826(1):56, 2016.
  * [4] J. W. Henning, J. T. Sayre, C. L. Reichardt, P. A. R. Ade, A. J. Anderson, J. E. Austermann, J. A. Beall, A. N. Bender, B. A. Benson, and L. E. Bleem.  Measurements of the Temperature and E-Mode Polarization of the CMB from 500 Square Degrees of SPTpol Data.  The Astrophysical Journal, 852(2), 2017.
  * [5] Beutler Florian, Blake Chris, Colless Matthew, Jones D. Heath, Staveley Smith Lister, Campbell Lachlan, Parker Quentin, Saunders Will, and Watson Fred.  The 6dF Galaxy Survey: baryon acoustic oscillations and the local Hubble constant.  Monthly Notices of the Royal Astronomical Society, 416(4):3017–3032, 2011.
  * [6] Vivian Poulin, Tristan L Smith, Tanvi Karwal, and Marc Kamionkowski.  Early Dark Energy Can Resolve The Hubble Tension, 2019.
  * [7] Edvard Mörtsell and Suhail Dhawan.  Does the Hubble constant tension call for new physics?  Journal of Cosmology and Astroparticle Physics, 2018(09):025–025, 2018.
  * [8] Stephen M Feeney, Hiranya V Peiris, Andrew R Williamson, Samaya M Nissanke, Daniel J Mortlock, Justin Alsing, and Dan Scolnic.  Prospects for resolving the hubble constant tension with standard sirens.  Phys. Rev. Lett., 122:061105, 2019.
  * [9] Riess, G. Adam, Casertano, Stefano, Yuan, Wenlong, Macri, M. Lucas, and Scolnic.  Large magellaniccloud cepheid standards provide a 1 The Astrophysical Journal, 876(1):85, 2019.
  * [10] Christian Corda.  The future of gravitational theories in the era of the gravitational wave astronomy.  International Journal of Modern Physics D, 27(05):1850060, 2018.
  * [11] Rong Gen Cai, Zhoujian Cao, Zong Kuan Guo, Shao Jiang Wang, and Tao Yang.  The gravitational wave physics.  National Science Review, 4(5), 2017.
  * [12] The Ligo Scientific Collaboration and The Virgo Collaboration.  Gw170817: Observation of gravitational waves from a binary neutron star inspiral.  Physical Review Letters, 119:161101, 2017.
  * [13] A. Goldstein, P. Veres, E. Burns, M. S. Briggs, and M. Stanbro.  An Ordinary Short Gamma-Ray Burst with Extraordinary Implications: Fermi -GBM Detection of GRB 170817A.  Astrophysical Journal, 848(2):L14, 2017.
  * [14] B. P. Abbott et al.  A gravitational-wave standard siren measurement of the Hubble constant.  Nature, 551(7678):85, 2017.
  * [15] Chris Messenger and Jocelyn Read.  Measuring a cosmological distance-redshift relationship using only gravitational wave observations of binary neutron star coalescences.  Phys. Rev. Lett., 108:091101, 2012.
  * [16] Sathyaprakash, B., S., Schutz, Bernard, and F.  Physics, astrophysics and cosmology with gravitational waves.  Living Reviews in Relativity, 12, 2009.
  * [17] L. Blanchet, T. Damour, and G. Schaefer.  Post-newtonian hydrodynamics and post-newtonian gravitational wave generation for numerical relativity.  Monthly Notices of the Royal Astronomical Society, 242(3):289–305, 1990.
  * [18] Lee Samuel Finn and David F. Chernoff.  Observing binary inspiral in gravitational radiation: One interferometer, 1993.
  * [19] K. G. Arun, Bala R Iyer, B. S. Sathyaprakash, and Pranesh A Sundararajan.  Parameter estimation of inspiralling compact binaries using 3.5 post-newtonian gravitational wave phasing: The nonspinning case.  Physical Review D, 71:084008, 2005.
  * [20] F. Douchin and P. Haensel.  A unified equation of state of dense matter and neutron star structure, 2001.
  * [21] The Ligo Scientific Collaboration and The Virgo Collaboration.  Gw170817: Observation of gravitational waves from a binary neutron star inspiral.  Physical Review Letters, 119:161101, 2017.
  * [22] Bo Wang, Zhenyu Zhu, Ang Li, and Wen Zhao.  Comprehensive analysis of the tidal effect in gravitational waves and implication for cosmology.  The Astrophysical Journal Supplement Series, 250(1):6, 2020.
  * [23] A. Abramovici, W. E. Althouse, R. W. P. Drever, Y. Gursel, S. Kawamura, F. J. Raab, D. Shoemaker, L. Sievers, R. E. Spero, and K. S. Thorne.  Ligo: The laser interferometer gravitational-wave observatory.  Science, 256(5055):325–333, 1992.
  * [24] Helstrom, Carl, W., Wilson, Fred, and L.  Statistical theory of signal detection.  Phys. Today, 23(6):73–73, 1970.
  * [25] Sebastian Khan, Sascha Husa, Mark Hannam, Frank Ohme, Michael Pürrer, Xisco Jiménez Forteza, and Alejandro Bohé.  Frequency-domain gravitational waves from non-precessing black-hole binaries. ii. a phenomenological model for the advanced detector era.  Physical Review D, 93:044007, 2016.
  * [26] Sascha, Husa, Sebastian, Khan, Mark, Hannam, Michael, Pürrer, Frank, and Ohme.  Frequency-domain gravitational waves from nonprecessing black-hole binaries. i. new numerical waveforms and anatomy of the signal.  Physical Review D, 93(4):44006–44006, 2016.
  * [27] S. and Chandrasekhar.  The maximum mass of ideal white dwarfs.  Journal of Astrophysics and Astronomy, 74:81, 1994.
  * [28] Thomas Bayes.  An essay towards solving a problem in the doctrine of chances.  Philosophical Transactions of the Royal Society of London, 53:370–418, 1763.
  * [29] Eric Thrane and Colm Talbot.  An introduction to bayesian inference in gravitational-wave astronomy: parameter estimation, model selection, and hierarchical models—corrigendum.  Publications of the Astronomical Society of Australia, 37, 2020.
  * [30] Rafael Rebolo.  Bayesian Astrophysics.  Canary Islands Winter School of Astrophysics. Cambridge University Press, 2018.