11institutetext: Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan
281, 9000 Gent, Belgium
11email<EMAIL_ADDRESS>
# Tetrahedral grids in Monte Carlo radiative transfer

Arno Lauwers  Maarten Baes  Peter Camps  Bert Vander Meulen

###### Abstract

Context. To understand the structures of complex astrophysical objects, 3D
numerical simulations of radiative transfer processes are invaluable. For
Monte Carlo radiative transfer, the most common radiative transfer method in
3D, the design of a spatial grid is important and non-trivial. Common choices
include hierarchical octree and unstructured Voronoi grids, each of which has
advantages and limitations. Tetrahedral grids, commonly used in ray-tracing
computer graphics, can be an interesting alternative option.

Aims. We aim to investigate the possibilities, advantages, and limitations of
tetrahedral grids in the context of Monte Carlo radiative transfer. In
particular, we want to compare the performance of tetrahedral grids to other
commonly used grid structures.

Methods. We implemented a tetrahedral grid structure, based on the open-source
library TetGen, in the generic Monte Carlo radiative transfer code SKIRT.
Tetrahedral grids can be imported from external applications or they can be
constructed and adaptively refined within SKIRT. We implemented an efficient
grid traversal method based on Plücker coordinates and Plücker products.

Results. The correct implementation of the tetrahedral grid construction and
the grid traversal algorithm in SKIRT were validated using 2D radiative
transfer benchmark problems. Using a simple 3D model, we compared the
performance of tetrahedral, octree, and Voronoi grids. With a constant cell
count, the octree grid outperforms the tetrahedral and Voronoi grids in terms
of traversal speed, whereas the tetrahedral grid is poorer than the other
grids in terms of grid quality. All told, we find that the performance of
tetrahedral grids is relatively poor compared to octree and Voronoi grids.

Conclusions. Although the adaptively constructed tetrahedral grids might not
be favourable in most media representative of astrophysical simulation models,
they still form an interesting unstructured alternative to Voronoi grids for
specific applications. In particular, they might prove useful for radiative
transfer post-processing of hydrodynamical simulations run on tetrahedral or
unstructured grids.

###### Key Words.:

radiative transfer

##  1 Introduction

Nearly all our astronomical knowledge is based on the observation and detailed
analysis of radiation we obtain from astronomical sources over the
electromagnetic spectrum. Through processes like absorption, scattering, and
re-emission, this radiation is affected and altered by the medium it passes
through. As a result, detailed modelling of these interactions, that is,
radiative transfer simulations, is required. This is particularly important if
we want to generate synthetic observations for snapshots of
(magneto)hydrodynamical simulations that can directly be compared to
observations.

The 3D radiative transfer problem is a notoriously hard nut to crack, and
various approaches have been proposed. The Monte Carlo method (Whitney 2011;
Steinacker et al. 2013; Noebauer & Sim 2019) is the most widely used method
for various types of radiative transfer problems. The method is intuitive,
conceptually straightforward, and, contrary to other methods such as finite
difference techniques, easily applicable to 3D problems. In the past few
years, advanced and powerful Monte Carlo radiative transfer codes have been
developed for a variety of radiative transfer applications, such as dust,
resonant line, atomic and molecular line, X-ray, UV ionising radiation, and
neutrino radiative transfer (e.g. Pinte et al. 2006; Dullemond et al. 2012;
Reissl et al. 2016; Harries et al. 2019; Li et al. 2020; Camps & Baes 2020;
Smith et al. 2023, to name just a few generic examples).

Three-dimensional Monte Carlo radiative transfer codes generally need a grid
onto which the medium is discretised. This grid contains a large number of
cells within which all physical properties are considered uniform. Early Monte
Carlo radiative transfer simulations often employed simple grids such as
spherical shells or regular Cartesian grids. For modern applications, and in
particular for the radiative transfer post-processing of hydrodynamical
simulations, this approach is not ideal, as these simulations typically cover
a large dynamic range in densities. The standard approach to dealing with this
complexity is to use hierarchical Cartesian grids, such as octree or binary
trees. Many Monte Carlo radiative transfer codes are equipped with such grids
(Kurosawa & Hillier 2001; Saftly et al. 2013, 2014; Juvela 2019).

More recently, unstructured Voronoi grids have gained popularity. Several
Monte Carlo radiative transfer codes have now been equipped with Voronoi grids
(Camps et al. 2013; Hubber et al. 2016; Koepferl et al. 2017; Smith et al.
2020; Byrohl et al. 2021; Tessore et al. 2021). The use of Voronoi grids in
radiative transfer, and in computational methods in general, has several
advantages. Voronoi meshes have the ability to deal with arbitrary complex
geometries with a large dynamical range, without any preferential directions.
Since a Voronoi grid is completely determined by the position of the
generating nodes, it is easy to adapt the cell sizes and grid resolution
according to user-defined criteria. In general, Voronoi grids allow the
dynamical range in densities to be closely reflected with as few cells as
needed.

Voronoi grids are not the only type of unstructured meshes used in
computational geometry. Tetrahedral meshes — tessellations in which all cells
are tetrahedra — are an interesting alternative. Tetrahedral grids are very
popular in computer graphics and 3D reconstruction applications (e.g. Georgii
& Westermann 2006; Marmitt & Slusallek 2006; Maria et al. 2017; Morrical et
al. 2023). A tetrahedron is the most elementary polyhedron: each tetrahedron
has only four vertices and six edges, and is bounded by four triangular faces.
In comparison, an average Voronoi cell in a 3D Voronoi tessellation has
approximately 27 vertices and 41 edges, and is bounded by approximately 16
convex polygonal faces (van de Weygaert 1994). Tetrahedral grids are therefore
‘simpler’ than Voronoi grids and more economical in terms of storage.

In addition to their ability to efficiently represent media with a large
dynamical range, the most important characteristic of grids in the context of
Monte Carlo radiative transfer is whether grid traversal can be implemented in
a robust and effective manner. Grid traversal essentially comes down to
determining, for a path defined by a starting position and direction, the
ordered list of cells crossed along this path and the individual path length
contributions covered within each individual cell. In many Monte Carlo
radiative transfer simulations, grid traversal is responsible for a large
fraction of the computational time.

For Voronoi grids, grid traversal is non-trivial and requires the calculation
of the intersection points of the path and each of the faces of the cell
(Camps et al. 2013; Hubber et al. 2016). One might expect that grid traversal
through a tetrahedral grid would be more efficient than through a Voronoi
grid. A tetrahedron has only four faces, so even using a naive algorithm in
which we calculate the intersection points of the path and each of the faces,
it can probably be quite efficient (remember that Voronoi cells have, on
average, almost four times as many faces). Moreover, efficient grid traversal
algorithms have been developed for tetrahedral grids in the field of computer
graphics (Platis & Theoharis 2003; Marmitt & Slusallek 2006; Lagae & Dutré
2008; Maria et al. 2017).

In this work, we explored the use of tetrahedral grids in Monte Carlo
radiative transfer. More specifically, we implemented tetrahedral grids and an
advanced grid traversal algorithm in the Monte Carlo radiative transfer code
SKIRT (Camps & Baes 2015, 2020). In Sects. 2 and 3 we discuss the construction
and traversal of tetrahedral grids in SKIRT, and in Sect 4 we validate our
implementation using established benchmark models. We then compare the
efficiency of tetrahedral and other grids in Sect. 5 and discuss these results
in a broader context in Sect. 6. Finally, we present our conclusions in Sect.
7.

##  2 Tetrahedral grid construction in SKIRT

###  2.1 The SKIRT radiative transfer code

SKIRT (Camps & Baes 2015, 2020) is a generic 3D Monte Carlo radiative transfer
code. Originally conceived as a dust radiative transfer code (Baes et al.
2003, 2011), it has transformed into a more generic tool that can be applied
to various radiative transfer applications. Apart from dust radiative
transfer, it can now perform Ly$\alpha$ resonant line radiative transfer
(Camps et al. 2021), non-local thermodynamic equilibrium line radiative
transfer for selected ions, atoms, and molecules (Matsumoto et al. 2023),
X-ray radiative transfer (Vander Meulen et al. 2023), and Hi radiative
transfer (Gebek et al. 2023). The efficiency of the code is maximised through
many optimisation mechanisms and a hybrid parallelisation strategy (Baes 2008;
Steinacker et al. 2013; Baes et al. 2016, 2022; Verstocken et al. 2017).

The SKIRT code is equipped with a diverse suite of options to discretise the
medium. Apart from the standard spherical, cylindrical and Cartesian grids,
the user can select octree grids (Saftly et al. 2013), $k$-d tree grids
(Saftly et al. 2014), and Voronoi grids (Camps et al. 2013). The availability
of these different options makes SKIRT a suitable environment to investigate
the accuracy and efficiency of tetrahedral grids.

###  2.2 Delaunay triangulation

A triangulation of a set of points is a subdivision of the convex hull of
those points into non-overlapping simplices. In 3D this is sometimes called a
tetrahedralisation as the simplices are tetrahedra. The Delaunay triangulation
(DT) of a set of 3D points $\\{\boldsymbol{p}_{i}\\}$ is the set of non-
overlapping tetrahedra whose vertices are points, $\boldsymbol{p}_{i}$, that
do not lie inside the circumsphere of any tetrahedron in the triangulation.
This last condition is also called the Delaunay condition. This ensures a
unique triangulation if and only if no more than four points lie on the same
sphere. However, even if five or more points lie on the same sphere, we can
usually still find a non-unique triangulation. An example of a DT can be seen
in Fig. 1.

Figure 1:  DT of a set of 30 random points drawn from a uniform distribution.

Generating the DT from a set of given input of points is a well-known problem
in computational geometry (Borouchaki & Lo 1995; Cignoni et al. 1998). The 2D
DT has the useful property that it maximises the smallest angle in each
triangle, resulting in regular triangles with nearly square bounding boxes. On
the other hand, the 3D DT does not have this property; cells tend to be
elongated and very irregular (Musin 1997; Chew 1997; Maur 2002). For many
applications an isotropic grid is preferred, making regular tetrahedra the
more desirable choice. In order to improve the quality of a tetrahedral grid
one can move the vertices, change the Delaunay condition, or add so-called
Steiner points in appropriate locations. Several techniques have already been
developed for this purpose, such as Delaunay refinement (Shewchuk 1970).

###  2.3 TetGen

Given the intricacy of the problem, we chose to use the third-party library
TetGen (Si 2015) for the construction of tetrahedral grids in SKIRT. TetGen is
a C++ quality tetrahedral mesh generator that is widely used in computational
fluid dynamics simulations (e.g. Shirokov 2021; Mayer et al. 2009; Bonfiglioli
et al. 2013).

TetGen uses a constrained Delaunay refinement algorithm (Shewchuk & Si 2014)
and offers various features for creating a robust computational grid to
perform radiative transfer on. First of all is the ability to generate
constrained triangulations, forcing grid cells to adhere to a specified input
surface. In our case this would typically be the simulation’s bounding box.
This means that the simulation’s domain is fully covered with grid cells,
rather than only partially by the convex hull of the triangulation.
Additionally, and most importantly, TetGen supports a mesh sizing function,
allowing users to define a function in the grid space that determines the cell
refinement. An example of this could be a scalar density field to determine
the desired resolution. In practice, TetGen achieves this by letting the user
set a function tetunsuitable, which should return true if the tetrahedron must
be optimised. This is usually done by applying various local mesh operations,
such as adding Steiner points. The tetunsuitable function can represent any
condition chosen by the user, such as a maximum ratio of edge lengths, volume,
or dihedral angle. If desired, it could also be made anisotropic. It should
also be noted that TetGen performs additional relaxing and smoothing
operations after applying the mesh sizing function to the entire domain. As a
result, the user-defined conditions are not strictly met in all cells. While
these operations can be disabled, doing so typically results in a lower-
quality grid.

###  2.4 TetGen in SKIRT

For our specific purpose of discretising a given density field in the SKIRT
radiative transfer code, we needed a recipe on which the refinement algorithm
is based. Mimicking the strategy applied for hierarchical grids in SKIRT
(Saftly et al. 2013, 2014), we chose limiting the total mass fraction inside a
cell as our default option. This naturally results in a higher resolution in
regions with a larger density. An example of this adaptive grid construction
is displayed in Fig. 2. We employed an additional subdivision criterion, which
is also used in the hierarchical tree grids in SKIRT as described by Saftly et
al. (2014). We subdivided the cell if the quantity

$q=\begin{cases}\;\dfrac{\rho_{\text{max}}-\rho_{\text{min}}}{\rho_{\text{max}}}&\text{if
}\rho_{\text{max}}>0,\\\ \;0&\text{if }\rho_{\text{max}}=0.\end{cases}$ (1)

exceeds some user-defined threshold $q_{\text{max}}\leq 1$. In this
expression, $\rho_{\text{min}}$ and $\rho_{\text{max}}$ represent the minimum
and maximum values of the density in the cell, which are estimated by
generating a large number of random positions in the cell (see Sect. 2.5).
This additional criterion avoids the use of cells with large density
gradients.

Figure 2: TetGen’s mesh sizing function applied in SKIRT. Right: Cut through
the density field for a dusty spiral model. Left: Cut through the
triangulation built using TetGen by constraining the minimum mass fraction per
tetrahedron to $10^{-5}$.

In addition to adaptively constructing the grid, tetrahedral grids can be
directly imported from hydrodynamical simulations or generated using other
tools. Furthermore, when supplied with particle data, such as from a smoothed-
particle hydrodynamic simulation, TetGen can generate the corresponding DT of
the point set. The uniform density in each cell can then be derived by
interpolating the input densities at the vertices.

###  2.5 Generating random points

After the construction of a grid, the properties of each cell, such as the
mass density, must be determined. In the SKIRT radiative transfer code, this
is achieved by sampling a user-defined number of random points in each cell
and averaging the underlying density at these locations. This requires an
algorithm to efficiently sample random points in an arbitrary tetrahedron. A
simple option could consist of generating random points inside the bounding
box of the tetrahedron and subsequently checking whether they fall inside the
tetrahedron. However, more efficient techniques for generating uniformly
distributed random points within a tetrahedron exist, and we used the
algorithm described by Visual Computing Lab (2011).

##  3 Grid traversal

###  3.1 Plücker coordinates

Plücker coordinates offer a method for representing directed lines in 3D space
using 6D vectors (Shoemake 1998; Platis & Theoharis 2003). Given a directed
line $R$ defined by a point $\boldsymbol{r}$ and a direction $\boldsymbol{k}$,
its Plücker coordinates can be expressed as

$\displaystyle\boldsymbol{\pi}_{R}=\\{\boldsymbol{k}:\boldsymbol{k}\times\boldsymbol{r}\\}=\\{\boldsymbol{U}_{R}:\boldsymbol{V}_{R}\\}.$
(2)

We can now define the permuted inner product or Plücker product of two
directed lines $R$ and $S$ as

$\displaystyle\boldsymbol{\pi}_{R}\odot\boldsymbol{\pi}_{S}=\boldsymbol{U}_{R}\cdot\boldsymbol{V}_{S}+\boldsymbol{U}_{S}\cdot\boldsymbol{V}_{R}.$
(3)

The sign of this product determines the relative orientation between two
directed lines. We have

$\begin{cases}\;\boldsymbol{\pi}_{R}\odot\boldsymbol{\pi}_{S}>0&\iff\;\textrm{$S$
goes anti-clockwise around $R$,}\\\
\;\boldsymbol{\pi}_{R}\odot\boldsymbol{\pi}_{S}<0&\iff\;\textrm{$S$ goes
clockwise around $R$,}\\\
\;\boldsymbol{\pi}_{R}\odot\boldsymbol{\pi}_{S}=0&\iff\;\textrm{$S$ intersects
or is parallel to $R$.}\end{cases}$ (4)

This allows us to efficiently find the exit face when tracing a photon path
through a tetrahedron, as will be discussed in the next section.

###  3.2 Traversal algorithm

We use the term ‘ray’ to indicate a given photon path represented by its
Plücker coordinates. Considering the edges of a tetrahedron as directed lines
from one vertex to another, a ray exits a face if the clockwise-
ordered111Viewed from the fourth vertex or also from inside. edges yield a
negative Plücker product with the ray. In other words, if the Plücker products
computed around clockwise-ordered edges move clockwise around the ray, we have
an exiting ray. Conversely, if they are all anti-clockwise, indicated by
positive Plücker products, we have an entering ray.

A naive approach would be to calculate these three Plücker products for each
of the four faces. However, despite the edges having a different orientation
when associated with each face, there is no need to recalculate the products,
as flipping the orientation of the edge simply flips the sign of the resulting
Plücker product. We can thus reduce the number of required Plücker products to
six, one for each edge. This approach was applied by Platis & Theoharis
(2003). However, we usually know the entry face of the ray and thus only need
to check the other three faces. This, with some additional optimisations, was
used by Marmitt & Slusallek (2006) to reduce the number of Plücker products to
$2.67$ on average. Other approaches exist, such as using three to six scalar
triple products, as described in Lagae & Dutré (2008).

On the other hand, Maria et al. (2017) described a very efficient method that
only requires two Plücker products. This is illustrated in Fig. 3. We label
the faces using the opposite vertices and use face $0$ as the entry face. We
started by computing the Plücker product of the edge $1\rightarrow 0$ with the
ray. If this results in a (anti-)clockwise orientation, we next used the edge
that is (anti-)clockwise viewed from vertex $0$. We computed the Plücker
product of this new edge and the ray. If these two calculated Plücker products
have an opposite sign, then the exit face is the face bounded by these two
edges. This can be understood if we remember that the edges of the entry face
have a known orientation with the ray222Anti-clockwise orientation with
clockwise ordered edges (opposite from the leaving ray)., as it is the
entering face. Using this additional information for the last edge, means we
have found three edges with a clockwise orientation in a clockwise order. If
instead they have an equal sign, then the remaining edge, which we do not
check, will have an opposite-signed Plücker product. This will make face $1$
the exit face, which can be understood with an argument analogous to the
above.

Figure 3: Algorithm, as described by Maria et al. (2017), that finds the exit
face (coloured) using only two Plücker products (left) and the corresponding
decision tree (right). The leaf nodes in the tree are the possible exit faces,
numbered according to the opposite vertices.

Since all the possible edges we have to check intersect at the same point,
vertex $0$, we can shift our coordinate frame such that the origin coincides
with this vertex. This will speed up the computation by ensuring the Plücker
coordinates, (2), of the edges have a vanishing moment $\boldsymbol{V}_{E}$.
We will, however, need to shift the ray to this new coordinate system by
updating its moment to

$\displaystyle\boldsymbol{V}_{R}=\boldsymbol{k}\times\left(\boldsymbol{r}-\boldsymbol{v}_{0}\right),$
(5)

where $\boldsymbol{v}_{0}$ are the coordinates of vertex $0$. Once the exit
face is found we can find the distance to this face using

$\displaystyle
s=\dfrac{\boldsymbol{n}\cdot\left(\boldsymbol{r}-\boldsymbol{v}_{0}\right)}{\boldsymbol{n}\cdot\boldsymbol{k}}$
(6)

where $\boldsymbol{n}$ is normal to the intersected face. We can then move the
ray forwards by a distance $s$.

However, when starting the algorithm we did not have the entering face
information. Thus, for the first step, we needed to take a slightly different
approach, which was also suggested by Maria et al. (2017). Once we determined
the tetrahedron containing the ray position, we calculated the Plücker product
with an arbitrary edge. This orientation indicates which of the two adjacent
faces cannot be the exit face. Subsequently, we can apply the previous
algorithm, using this face as the ‘fake’ entry face. The algorithm will still
work even if this face is not the true entering face.

We should also mention that there are a few (rarely occurring) edge cases,
such as the ray going through a vertex or edge, that need to be handled with
care. We used a simple plane intersection test whenever the leaving face might
be ambiguous to resolve these edge cases.

###  3.3 Locating the ray origin

A Monte Carlo radiative transfer simulation differs from a traditional ray-
tracing simulation in that each ray’s origin is generated randomly, so that
the cell containing this position needs to be determined every time a ray is
fired. To accomplish this, we adopted a strategy inspired by Camps et al.
(2013) and constructed a $k$-d search tree using the centroids of the
tetrahedra. When we wanted to find the cell containing a given point, we
located the closest centroid using the $k$-d tree and used our traversal
algorithm to traverse from the found centroid towards the given input point.
It is important to note that the nearest centroid often does not belong to the
tetrahedron containing the input point. If, during traversal, the distance
travelled exceeds the distance between the found centroid and the input point,
we have found the desired tetrahedron. In the very rare event the algorithm
fails, we can resort to a loop to check all tetrahedra.

It is also a viable option to use a breadth-first search over the neighbours
around the found cell. We conducted several tests and found that using the
traversal algorithm is very reliable and the most efficient.

##  4 Validation

To verify our implementation, we performed two dust radiative transfer
benchmarks, comparing our results with those calculated using a regular
cylindrical SKIRT grid as well as with results from other radiative transfer
codes.

Figure 4: Results from the benchmark described by Pascucci et al. (2004). The
simulation was performed on both a regular (cylindrical) and a tetrahedral
grid in SKIRT and using six other radiative transfer codes. The figure shows
the SED for two inclinations, $i=12.5^{\circ}$ and $i=77.5^{\circ}$, relative
to the disc, and for optical depths of $0.1$, $1$, $10$, and $100$. The solid
and dotted lines represent the results using our grid and the regular grid in
SKIRT, respectively. The coloured areas represent the minimum and maximum from
the other radiative transfer codes’ benchmark results.

Pascucci et al. (2004): The Pascucci et al. (2004) benchmark features a
central T Tauri star embedded in a circumstellar disc with an inner cavity
free of dust. The density structure of the disc represents that of a massless
Keplerian disc. The challenge of the radiative transfer problem is that the
density distribution provides a significant density gradient in the inner part
of the disc. The benchmark performs a panchromatic simulation and focuses on
the spectral energy distributions (SEDs); it is thus less sensitive to the
individual paths taken and will validate other aspects, such as calculated
path lengths and absorption. In the original work, five different radiative
transfer codes participated to the benchmark, and all of them reproduce the
SEDs very well, with relative differences generally below 10%, except for the
most challenging cases (high optical depth and edge-on orientation). We ran
the benchmark simulations with SKIRT, using both tetrahedral and regular
cylindrical grids, with $30$ and half a million cells, respectively. As shown
in Fig. 4, our results reproduce the benchmark results very well, with an
absolute relative error of $1.5\%$ and $2\%$ for $i=12.5^{\circ}$ and
$i=77.5^{\circ}$, respectively.

Figure 5: Results from the benchmark described by Pinte et al. (2009). The
simulation was performed on both a regular cylindrical and a tetrahedral grid
in SKIRT and using three other radiative transfer codes. The polarisation
degree and the difference to the means of all benchmark results are shown for
the two inclinations $i=69.5^{\circ}$ and $i=87.1^{\circ}$ relative to the
disc.

Pinte et al. (2009): The Pinte et al. (2009) benchmark builds on the Pascucci
et al. (2004) benchmark but extends it to more challenging values. The optical
depth extends to larger values (up to $10^{6}$ in the midplane), scattering is
anisotropic instead of isotropic, and scattering polarisation is included.
With images using SKIRT polarisation (Peest et al. 2017), this benchmark
provides a very strong test for the accuracy of the individual photon paths.
The results for this benchmark are depicted in Fig. 5. The SKIRT results
corresponding to the cylindrical and tetrahedral grids are in close agreement,
with an average error of only $0.3\%$ and $0.5\%$ for $i=69.5^{\circ}$ and
$i=87.1^{\circ}$, respectively. The tetrahedral and cylindrical grid
respectively used $41$ and $6$ million cells. The artefacts observed in the
outer regions of the polarisation maps stem from using a photon count that was
not scaled with the amount of cells required by our 3D grid. When attempting
to accurately represent the 2D medium with a 3D grid, a significantly greater
number of cells was necessary compared to the larger cylindrical cells.
However, we did not proportionately increase the photon count, resulting in
fewer photons per cell within the 3D grid. Consequently, in the outer regions
where the flux is typically 20 orders of magnitude lower than in the central
region, some cells end up with zero flux, leading to an undefined polarisation
degree in those areas.

##  5 Comparison to other grids

We next compared the performance of the tetrahedral grid to that of the
unstructured Voronoi grid and the hierarchical octree grid already present in
SKIRT. To this end, we constructed a simple but fully 3D model, and we ran
simulations with various grid types and cell counts. We measured the
simulation run time and evaluated the grid quality using a metric based on
synthetic observations of this model.

The simulation model used for our tests is defined in a cubical spatial
domain. The cube contains a dust medium defined by ten superposed, arbitrarily
placed Plummer spheres with varying radii and dust masses. The spheres are cut
off at the boundaries of the cube. There are also three random point sources
placed at arbitrary locations within the cube and thus embedded in the medium.
We performed monochromatic simulations that generate parallel-projected 100 by
100 pixel images of this system for three different sight lines. We
exclusively used the scattered flux to eliminate the error from the bright
point sources in the images. One of these images is shown in Fig. 6. The dust
density in the cube varies by four orders of magnitude, and the optical depth
from the centre of the cube to its boundary ranges from 2.5 to 450. While this
geometry is quite artificial, for our purposes it should be fairly
representative of many astrophysical simulation models.

We ran simulations of the test model with varying cell counts. Constructing a
Voronoi grid with a given cell count is straightforward: we can simply specify
the number of Voronoi sites to be sampled from the dust distribution. However,
tetrahedral and octree grids are created adaptively. For these grid types, we
varied the maximum dust mass fraction allowed in a grid cell to control the
cell count. A stricter mass fraction limit will cause extra cells, especially
in regions with higher density.

The leftmost panel of Fig. 7 shows the runtime of our test model simulations
as a function of the number of grid cells. As expected, the slopes of the
curves are similar for the three grid types. At constant cell count, the
octree grid significantly outperforms the other grids. This is unsurprising
given the simple Cartesian cell geometry and in view of the extensively
optimised implementation as documented by Saftly et al. (2013, 2014). On the
other hand, at constant cell count, the Voronoi grid is substantially slower
than the tetrahedral grid. Again, this is not surprising because, compared to
tetrahedrons, Voronoi cells have a more complex geometry and on average have
more neighbouring cells.

Comparing run time at constant cell count tells only half of the story. Some
grid types may require fewer cells, and thus a correspondingly shorter
runtime, to achieve the same results. We thus needed a metric to evaluate the
‘quality’ of the spatial grid employed in the simulation. This is ultimately
determined by the accuracy of the generated synthetic observables. We
therefore defined a quality measure based on the images generated by our test
setup. We first produced a reference image for each sight line by running the
simulation using a very high-resolution octree grid and a large number of
photon packets, ensuring that the fluxes in each pixel are numerically
converged. To derive a quality metric for a given test run, we calculated the
absolute value of the relative error on the flux in each pixel, as compared to
the corresponding reference pixel. The final quality metric is then defined as
the mean value of these relative errors for all pixels in the image.

The rightmost panel of Fig. 7 shows this quality metric for the three sight
lines of our test model simulations, as a function of the number of grid
cells. At constant cell count, the tetrahedral grid produces a significantly
higher relative error than the other grid types. Or equivalently, for a given
level of accuracy, the tetrahedral grid requires a substantially larger number
of cells. Notably, the sight line along the Cartesian axes ($i{=}0^{\circ}$
$a{=}0^{\circ}$) shows an overall worse quality, especially for the octree
grid. This can be attributed to discretisation effects caused by having both
the sight line and all cells in the grid lined up with the coordinate axes.
This occurs for both the octree grid and high-resolution octree reference,
which explains why this can also be seen in the unstructured grids.

The grid performance is compared more directly in Fig. 8, where we plot the
relative error as a function of simulation runtime. To avoid an unfair
comparison, we excluded the sight line at zero inclination and the azimuth
from this summary plot, instead averaging the relative errors for the two
other sight lines. For a given runtime, the tetrahedral grid consistently
produces a larger relative error than the other grid types, at least in the
relevant range where the relative error is below 10 per cent. Compared to the
octree grid, the tetrahedral grid requires more cells for the same quality and
it runs more slowly for an equal cell count (see Fig. 7), which causes the
substantial overall performance gap shown in Fig. 8. Compared to the Voronoi
grid, the tetrahedral grid runs faster for equal cell count, but this
advantage is insufficient to compensate for the extra number of cells required
for the same quality.

The conclusions drawn using our smooth test model cannot necessarily be
generalised to simulation models with other properties. For example, we expect
that unstructured grids will perform much better for models with steep density
gradients, resembling hard boundaries such as those often encountered in
computational fluid dynamics, where tetrahedral grids are prevalent. This is
especially true if the gradients are not lined up with the coordinate axes,
putting the Cartesian octree grid at a disadvantage. Furthermore, constructing
the tetrahedral grid using other refinement options could also improve its
performance.

Figure 6: Image of the scattered flux from our test model with ten Plummer
spheres of dust and three point sources. The dynamic range of the flux is
close to six orders of magnitude.

Figure 7: Comparison of the tetrahedral grid with the octree Voronoi and
grids. Left: Simulation runtime of our test model as a function of the number
of grid cells. Right: Mean of the relative flux residuals, compared pixel by
pixel to the corresponding reference image for three sight lines, as a
function of the number of grid cells. Figure 8: Comparison of the tetrahedral
grid with the Voronoi and octree grids, now showing the relative error (the
mean of the pixel-by-pixel flux residuals) as a function of simulation
runtime. The colour indicates the number of grid cells used.

##  6 Discussion

###  6.1 Hierarchical versus unstructured grids

Tetrahedral grids form a potentially interesting spatial grid structure for
Monte Carlo radiative transfer, with advantages and disadvantages with respect
to other grid types. One of the interesting characteristics of tetrahedral
grids is that they are unstructured. In this sense, they are similar to
Voronoi grids, which have been implemented in several Monte Carlo radiative
transfer codes (Camps et al. 2013; Hubber et al. 2016; Koepferl et al. 2017;
Smith et al. 2020; Byrohl et al. 2021; Tessore et al. 2021). One of the
advantages of unstructured grids is that they are more flexible and adaptable,
and the individual grids cells are not aligned along any preferential
directions. In other words, unstructured grids provide a tessellation of the
medium without the rigidity of a Cartesian grid.

Another advantage is that they can be made to conform to any desired geometry
in such a way that a minimum number of grid cells is required for a given
spatially varying resolution. In the context of radiative transfer, this is
important for two different reasons. First, many radiative transfer problems
require the calculation of the medium’s secondary emission. This requires the
measurement of the radiation field in every cell with sufficient statistics,
which, in turn, requires sufficient photon packets passing through each cell.
Secondly, the grid traversal is not the only and not necessarily the most
time-consuming aspect of Monte Carlo radiation transfer. Sometimes the
calculation of the emissivity of the grid cells can be a substantial
computational effort too. A typical example is thermal emission by dust, which
in the case of multiple populations of transiently heated grains can be a
computationally costly operation (Désert et al. 1986; Draine & Li 2001; Camps
et al. 2015). Another example is the calculation of the level populations for
atoms or molecules in the frame of Monte Carlo non-local thermodynamic
equilibrium line radiative transfer (Matsumoto et al. 2023). For both of these
reasons, minimising the number of grid cells for a given resolution is thus
beneficial for the simulation run time.

Unstructured meshes also have some drawbacks. They are generally more complex
to build and store than hierarchical grids and require more memory usage. In
the context of Monte Carlo radiative transfer, the grid traversal time is
particularly important. Grid traversal through a hierarchical grid is
generally fast and it can be optimised using dedicated methods (Saftly et al.
2013). For Voronoi grids the traversal is relatively slow since the grid cells
have many faces, and each of these faces needs to be tested to determine the
exit point for a given photon path (Camps et al. 2013). For tetrahedral grids,
grid traversal is significantly faster than in Voronoi grids. In a tetrahedral
grid, each cell has obviously only four faces, compared to 16 faces on average
for a grid cell in a 3D Voronoi mesh (van de Weygaert 1994). Moreover,
efficient algorithms based on Plücker coordinates and Plücker products have
been designed for grid traversal through tetrahedral grids. For the same
number of cells, however, the traversal is still not as efficient as in
hierarchical octree grids.

###  6.2 Post-processing hydrodynamical simulations

The main reason why unstructured grids are interesting in the context of Monte
Carlo radiative transfer is linked to post-processing hydrodynamical
simulations. Full-scale hydrodynamical simulations are a powerful tool for
studying the evolution of nearly all astronomical objects, and radiative
transfer post-processing of such simulation outputs is needed to generate
realistic mock observables that can be compared to real observations. The
increased popularity of radiative transfer on Voronoi grids is primarily
driven by the emergence of Lagrangian hydrodynamical solvers that operate on
such unstructured Voronoi grids. The AREPO code (Springel 2010; Weinberger et
al. 2020) is the most widely used one, but other solvers have been developed
(Duffell & MacFadyen 2012; Yalinewich et al. 2015; Vandenbroucke & De Rijcke
2016; Chang et al. 2017; Vandenbroucke & Wood 2018) or are being developed
(Schaller et al. 2023). It is beneficial to use the same grid for the
hydrodynamical calculations and the radiative transfer post-processing.
Voronoi-based radiative transfer post-processing has been applied to generate
mock observations for galaxies generated from moving-mesh cosmological
hydrodynamical simulations, such as IllustrisTNG and AURIGA (e.g. Rodriguez-
Gomez et al. 2019; Schulz et al. 2020; Camps et al. 2021; Popping et al. 2022;
Guzmán-Ortega et al. 2023; Bottrell et al. 2024).

The hydrodynamical equations can also be solved on tetrahedral grids. There
has been a long history of tetrahedral grids in computational fluid dynamics
simulations, in particular in engineering and material science applications
(Biswas & Strawn 1994, 1996; Löhner & Baum 1992; Ng et al. 2009; Xie et al.
2014, 2022). The advantage of tetrahedral meshes is that they are very
flexible and suited to fit within complex shapes with hard boundaries, such as
buildings, turbines, or airplane wings. In astrophysical hydrodynamics we
usually have no strong fixed boundaries for the computational domain.
Nevertheless, hydrodynamical simulations on tetrahedral meshes have also been
explored for astrophysical applications (Bruls et al. 1999; Lionello & Schnack
2002; Kulikov et al. 2020; Kulikov 2021a, b). If one wants to post-process
such simulation outputs using radiative transfer on the same grid as used for
the hydrodynamical simulations, the methodology discussed in this paper is
required.

Tetrahedral grids might also prove useful for post-processing moving-mesh
hydrodynamical simulations. However, in Monte Carlo radiative transfer it is
commonly assumed that all physical quantities, such as density, are uniform
within each cell (Whitney 2011; Steinacker et al. 2013). In finite-volume
hydrodynamics, of which moving-mesh hydrodynamics is a special case, the
hydrodynamical simulations are solved using Gudonov’s method, that is, by
considering Riemann problems at all cell boundaries (Toro 1997; Springel
2010). Rather than the total masses in each Voronoi cell, the densities at the
mesh-generating points are fundamental properties, and the assumption of a
uniform density within each cell (the standard assumption in radiative
transfer post-processing) is not consistent with the hydrodynamics method.

In principle, we can avoid this inconsistency by performing the radiative
transfer on the dual Delaunay tetrahedral grid rather than on the original
Voronoi grid. Instead of assuming a uniform density in each of the tetrahedra
cells, we can interpolate the density at each point based on the density
values at the vertices of each tetrahedron, which are, by definition, the
generating points of the moving-mesh simulation. An approach along this line
has been implemented for Cartesian grids in the original version of SKIRT by
Baes et al. (2003). Rather than assigning a uniform density in each cuboidal
grid cell, we used trilinear interpolation based on the values of the density
at the eight vertices. As a result, the optical depth along a given path
increases as a cubic function rather than as a linear function of the distance
travelled along the path within each cell. This increased the run time for
grid traversal by a factor of three. On the other hand, for a fixed number of
cells the accuracy of the density of the interpolated grid is much larger,
which compensates the longer run times.

The main drawback of applying a similar approach to tetrahedral grids is its
complexity. A planar triangulation of $n$ points results in ${\cal{O}}(n)$
triangles, but this does not apply to higher dimensions, where it is super-
linear with an ${\cal{O}}(n^{2})$ complexity.

##  7 Conclusions

The proper discretisation of Monte Carlo radiative transfer simulations is
crucial for obtaining accurate results. However, it is an intricate and
complicated problem, as the ‘ideal’ spatial grid should satisfy different,
sometimes conflicting, criteria. Research into the advantages and
disadvantages of different types of grids has been an important topic in the
framework of the SKIRT Monte Carlo radiative transfer code (e.g., Baes et al.
2003; Camps et al. 2013; Saftly et al. 2013, 2014).

Tetrahedral grids have long been used in engineering applications, and more
recently have been used for ray tracing in computer graphics. Coupled with its
use in fluid and hydrodynamical simulations, these grid are promising tools in
the context of Monte Carlo radiative transfer.

We have successfully implemented a working tetrahedral grid structure within
SKIRT. Our implementation uses the popular TetGen quality tetrahedral mesh
generator (Si 2015), while for the grid traversal we implemented a very
efficient method based on Plücker coordinates and Plücker products (Shoemake
1998; Maria et al. 2017).

The main conclusions of this work are the following:

  * •

Tetrahedral grids are a very flexible grid structure that can be applied in a
variety of simulations. The grid can be adaptively constructed using an
underlying input medium, directly imported from an external source, or built
from particle data.

  * •

The correct implementation of the grid construction and the grid traversal
algorithm were validated using 2D radiative transfer benchmark problems
(Pascucci et al. 2004; Pinte et al. 2009). The SKIRT results based on
tetrahedral grids and on regular 2D cylindrical grids reproduce the benchmark
results very well.

  * •

We compared the performance of tetrahedral, octree, and Voronoi grids using a
simple 3D model. We find that, with a constant cell count, grid traversal
through an octree grid is significantly more efficient than through a
tetrahedral grid, which in turn is more efficient than a Voronoi grid. On the
other hand, the tetrahedral grid has a lower grid quality compared to the
other grid types. Combining these results, we find that, for a given runtime,
the tetrahedral grid is consistently a lower-quality option than the other
grid types.

  * •

Despite the relatively poor performance of tetrahedral grids compared to
Voronoi and octree grids, tetrahedral grids remain an interesting option for
radiative transfer simulations, in particular for post-processing
hydrodynamical simulations run on tetrahedral or unstructured grids.

###### Acknowledgements.

B.V. acknowledges support by the Fund for Scientific Research Flanders (FWO-
Vlaanderen, project 11H2121N).

## References

  * Baes (2008) Baes, M. 2008, MNRAS, 391, 617
  * Baes et al. (2022) Baes, M., Camps, P., & Matsumoto, K. 2022, A&A, 666, A101
  * Baes et al. (2003) Baes, M., Davies, J. I., Dejonghe, H., et al. 2003, MNRAS, 343, 1081
  * Baes et al. (2016) Baes, M., Gordon, K. D., Lunttila, T., et al. 2016, A&A, 590, A55
  * Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22
  * Biswas & Strawn (1994) Biswas, R. & Strawn, R. C. 1994, Applied Numerical Mathematics, 13, 437, 10.1016/0168-9274(94)90007-8
  * Biswas & Strawn (1996) Biswas, R. & Strawn, R. C. 1996, Applied Numerical Mathematics, 20, 337, 10.1016/0168-9274(95)00099-2
  * Bonfiglioli et al. (2013) Bonfiglioli, A., Grottadaurea, M., Paciorri, R., & Sabetta, F. 2013, Computers & Fluids, 73, 162, 10.1016/j.compfluid.2012.12.022
  * Borouchaki & Lo (1995) Borouchaki, H. & Lo, S. 1995, Computer Methods in Applied Mechanics and Engineering, 128, 153, 10.1016/0045-7825(95)00854-1
  * Bottrell et al. (2024) Bottrell, C., Yesuf, H. M., Popping, G., et al. 2024, MNRAS, 527, 6506
  * Bruls et al. (1999) Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233
  * Byrohl et al. (2021) Byrohl, C., Nelson, D., Behrens, C., et al. 2021, MNRAS, 506, 5129
  * Camps & Baes (2015) Camps, P. & Baes, M. 2015, Astronomy and Computing, 9, 20
  * Camps & Baes (2020) Camps, P. & Baes, M. 2020, Astronomy and Computing, 31, 100381
  * Camps et al. (2013) Camps, P., Baes, M., & Saftly, W. 2013, A&A, 560, A35
  * Camps et al. (2021) Camps, P., Behrens, C., Baes, M., Kapoor, A. U., & Grand, R. 2021, ApJ, 916, 39
  * Camps et al. (2015) Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87
  * Chang et al. (2017) Chang, P., Wadsley, J., & Quinn, T. R. 2017, MNRAS, 471, 3577
  * Chew (1997) Chew, L. P. 1997, Guaranteed-quality delaunay meshing in 3d (short version), in Proceedings of the thirteenth annual symposium on Computational geometry, 391–393, 10.1145/262839.263018
  * Cignoni et al. (1998) Cignoni, P., Montani, C., & Scopigno, R. 1998, Computer-Aided Design, 30, 333, 10.1016/S0010-4485(97)00082-1
  * Désert et al. (1986) Désert, F. X., Boulanger, F., & Shore, S. N. 1986, A&A, 160, 295
  * Draine & Li (2001) Draine, B. T. & Li, A. 2001, ApJ, 551, 807, 10.1086/320227
  * Duffell & MacFadyen (2012) Duffell, P. C. & MacFadyen, A. I. 2012, ApJ, 755, 7
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library, 1202.015
  * Gebek et al. (2023) Gebek, A., Baes, M., Diemer, B., et al. 2023, MNRAS, 521, 5645
  * Georgii & Westermann (2006) Georgii, J. & Westermann, R. 2006, IEEE transactions on visualization and computer graphics, 12, 1345–1352, 10.1109/tvcg.2006.110
  * Guzmán-Ortega et al. (2023) Guzmán-Ortega, A., Rodriguez-Gomez, V., Snyder, G. F., Chamberlain, K., & Hernquist, L. 2023, MNRAS, 519, 4920
  * Harries et al. (2019) Harries, T. J., Haworth, T. J., Acreman, D., Ali, A., & Douglas, T. 2019, Astronomy and Computing, 27, 63
  * Hubber et al. (2016) Hubber, D. A., Ercolano, B., & Dale, J. 2016, MNRAS, 456, 756
  * Juvela (2019) Juvela, M. 2019, A&A, 622, A79
  * Koepferl et al. (2017) Koepferl, C. M., Robitaille, T. P., Dale, J. E., & Biscani, F. 2017, ApJS, 233, 1
  * Kulikov (2021a) Kulikov, I. 2021a, Molecular cloud collapse to stellar densities: Models on moving geodesic vs. unstructured tetrahedron vs. nested meshes, in Journal of Physics: Conference Series, Vol. 2028, IOP Publishing, 012001, 10.1088/1742-6596/2028/1/012001
  * Kulikov (2021b) Kulikov, I. 2021b, Journal of Applied and Industrial Mathematics, 15, 277, 10.1134/S1990478921020083
  * Kulikov et al. (2020) Kulikov, I., Vorobyov, E., Chernykh, I., & Elbakyan, V. 2020, Hydrodynamic modeling of self-gravitating astrophysical objects on tetrahedral meshes, in Journal of Physics: Conference Series, Vol. 1640, IOP Publishing, 012003, 10.1088/1742-6596/1640/1/012003
  * Kurosawa & Hillier (2001) Kurosawa, R. & Hillier, D. J. 2001, A&A, 379, 336
  * Lagae & Dutré (2008) Lagae, A. & Dutré, P. 2008, Computer Graphics Forum, 27, 1303, https://doi.org/10.1111/j.1467-8659.2008.01269.x
  * Li et al. (2020) Li, Y., Gu, M. F., Yajima, H., Zhu, Q., & Maji, M. 2020, MNRAS, 494, 1919
  * Lionello & Schnack (2002) Lionello, R. & Schnack, D. D. 2002, MH4D, an MHD Algorithm on an Unstructured Tetrahedral Grid., in AGU Fall Meeting Abstracts, Vol. 2002, SH52A–0492
  * Löhner & Baum (1992) Löhner, R. & Baum, J. D. 1992, International Journal for Numerical Methods in Fluids, 14, 1407, 10.1002/fld.1650141204
  * Maria et al. (2017) Maria, M., Horna, S., & Aveneau, L. 2017, Efficient Ray Traversal of Constrained Delaunay Tetrahedralization, in Proceedings of the 12th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications, Vol. 1, 12th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2017), Porto, Portugal, 236–243, hal-01486575
  * Marmitt & Slusallek (2006) Marmitt, G. & Slusallek, P. 2006, Fast ray traversal of tetrahedral and hexahedral meshes for direct volume rendering, in Proceedings of the Eighth Joint Eurographics / IEEE VGTC Conference on Visualization, EUROVIS’06 (Goslar, DEU: Eurographics Association), 235–242, isbn: 3905673312
  * Matsumoto et al. (2023) Matsumoto, K., Camps, P., Baes, M., et al. 2023, A&A, 678, A175
  * Maur (2002) Maur, P. 2002, Delaunay triangulation in 3d, Tech. Rep. DCSE/TR-2002-02, University of West Bohemia in Pilsen
  * Mayer et al. (2009) Mayer, U. M., Gerstenberger, A., & Wall, W. A. 2009, International journal for numerical methods in engineering, 79, 846, 10.1002/nme.2600
  * Morrical et al. (2023) Morrical, N., Sahistan, A., Güdükbay, U., Wald, I., & Pascucci, V. 2023, IEEE Transactions on Visualization and Computer Graphics, 29, 537, 10.1109/TVCG.2022.3209418
  * Musin (1997) Musin, O. R. 1997, Properties of the Delaunay triangulation, in Proceedings of the thirteenth annual symposium on Computational geometry, 424–426
  * Ng et al. (2009) Ng, Y. L., Yusoff, M. Z., & Shuaib, N. H. 2009, Tetrahedral mesh generator for CFD simulation of complex geometry, in 2009 3rd International Conference on Energy and Environment (ICEE), 330–336, 10.1109/ICEENVIRON.2009.5398626
  * Noebauer & Sim (2019) Noebauer, U. M. & Sim, S. A. 2019, Living Reviews in Computational Astrophysics, 5, 1
  * Pascucci et al. (2004) Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793
  * Peest et al. (2017) Peest, C., Camps, P., Stalevski, M., Baes, M., & Siebenmorgen, R. 2017, A&A, 601, A92
  * Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967
  * Pinte et al. (2006) Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797
  * Platis & Theoharis (2003) Platis, N. & Theoharis, T. 2003, Journal of Graphics Tools, 8, 10.1080/10867651.2003.10487593
  * Popping et al. (2022) Popping, G., Pillepich, A., Calistro Rivera, G., et al. 2022, MNRAS, 510, 3321
  * Reissl et al. (2016) Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2019, MNRAS, 483, 4140
  * Saftly et al. (2014) Saftly, W., Baes, M., & Camps, P. 2014, A&A, 561, A77
  * Saftly et al. (2013) Saftly, W., Camps, P., Baes, M., et al. 2013, A&A, 554, A10
  * Schaller et al. (2023) Schaller, M., Borrow, J., Draper, P. W., et al. 2023, MNRAS, submitted (arXiv:2305.13380)
  * Schulz et al. (2020) Schulz, S., Popping, G., Pillepich, A., et al. 2020, MNRAS, 497, 4773
  * Shewchuk (1970) Shewchuk, J. 1970, Proc. 14th Annu. ACM Sympos. Comput. Geom., 10.1145/276884.276894
  * Shewchuk & Si (2014) Shewchuk, J. & Si, H. 2014, Proceedings of the Annual Symposium on Computational Geometry, 10.1145/2582112.2582138
  * Shirokov (2021) Shirokov, I. A. 2021, Mathematical Models and Computer Simulations, 13, 1148, 10.1134/S2070048221060211
  * Shoemake (1998) Shoemake, K. 1998, Ray Tracing News, 11, http://www.acm.org/tog/resources/RTNews/html/rtnv11n1.html#art3
  * Si (2015) Si, H. 2015, ACM Trans. Math. Softw., 41, 11, 10.1145/2629697
  * Smith et al. (2020) Smith, A., Kannan, R., Tsang, B. T. H., Vogelsberger, M., & Pakmor, R. 2020, ApJ, 905, 27
  * Smith et al. (2023) Smith, A., Safranek-Shrader, C., Bromm, V., et al. 2023, COLT: Monte Carlo radiative transfer and simulation analysis toolkit, Astrophysics Source Code Library, record ascl:2306.034
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791
  * Steinacker et al. (2013) Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63
  * Tessore et al. (2021) Tessore, B., Pinte, C., Bouvier, J., & Ménard, F. 2021, A&A, 647, A27
  * Toro (1997) Toro, E. 1997, Riemann Solvers and Numerical Methods for Fluid Dynamics (Springer, Berlin), 10.1007/978-3-662-03490-3_1
  * van de Weygaert (1994) van de Weygaert, R. 1994, A&A, 283, 361
  * Vandenbroucke & De Rijcke (2016) Vandenbroucke, B. & De Rijcke, S. 2016, Astronomy and Computing, 16, 109
  * Vandenbroucke & Wood (2018) Vandenbroucke, B. & Wood, K. 2018, Astronomy and Computing, 23, 40
  * Vander Meulen et al. (2023) Vander Meulen, B., Camps, P., Stalevski, M., & Baes, M. 2023, A&A, 674, A123
  * Verstocken et al. (2017) Verstocken, S., Van De Putte, D., Camps, P., & Baes, M. 2017, Astronomy and Computing, 20, 16
  * Visual Computing Lab (2011) Visual Computing Lab. 2011, Generating Random Points in a Tetrahedron, accessed: January 23, 2024
  * Weinberger et al. (2020) Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32
  * Whitney (2011) Whitney, B. A. 2011, Bulletin of the Astronomical Society of India, 39, 101
  * Xie et al. (2022) Xie, B., Huang, Y., & Xiao, F. 2022, Journal of Computational Physics, 463, 111299, 10.1016/j.jcp.2022.111299
  * Xie et al. (2014) Xie, B., Ii, S., Ikebata, A., & Xiao, F. 2014, Journal of Computational Physics, 277, 138, https://doi.org/10.1016/j.jcp.2014.08.011
  * Yalinewich et al. (2015) Yalinewich, A., Steinberg, E., & Sari, R. 2015, ApJS, 216, 35


Our vertex-based local editing achieves realistic global context-aware mesh
deformation using a few handle points. The user specifies a few mesh vertices
(red) as handles, their target positions (blue), and a region of influence
(light green).

#  DragD3D: Vertex-based Editing for Realistic Mesh Deformations using 2D
Diffusion Priors

Tianhao Xie1, Eugene Belilovsky1,2, Sudhir Mudur1 and Tiberiu Popa1
1Concordia University, Montréal, Canada 2MILA, Montréal, Canada

###### Abstract

Direct mesh editing and deformation are key components in the geometric
modeling and animation pipeline. Direct mesh editing methods are typically
framed as optimization problems combining user-specified vertex constraints
with a regularizer that determines the position of the rest of the vertices.
The choice of the regularizer is key to the realism and authenticity of the
final result. Physics and geometry-based regularizers are not aware of the
global context and semantics of the object, and the more recent deep learning
priors are limited to a specific class of 3D object deformations. In this
work, our main contribution is a local mesh editing method called DragD3D for
global context-aware realistic deformation through direct manipulation of a
few vertices. DragD3D is not restricted to any class of objects. It achieves
this by combining the classic geometric ARAP (as rigid as possible)
regularizer with 2D priors obtained from a large-scale diffusion model.
Specifically, we render the objects from multiple viewpoints through a
differentiable renderer and use the recently introduced DDS loss which scores
the faithfulness of the rendered image to one from a diffusion model. DragD3D
combines the approximate gradients of the DDS with gradients from the ARAP
loss to modify the mesh vertices via neural Jacobian field, while also
satisfying vertex constraints. We show that our deformations are realistic and
aware of the global context of the objects, and provide better results than
just using geometric regularizers.

<ccs2012> <concept>
<concept_id>10010147.10010371.10010396.10010398</concept_id>
<concept_desc>Computing methodologies Mesh geometry models</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10010147.10010257</concept_id> <concept_desc>Computing
methodologies Machine learning</concept_desc>
<concept_significance>500</concept_significance> </concept> </ccs2012>

[500]Computing methodologies Mesh geometry models [500]Computing methodologies
Machine learning

††volume: 43††issue: 2

##  1 Introduction

Geometric deformation and editing are fundamental operations in the geometric
modeling pipeline that have received a lot of attention over the years. Among
the many varieties of geometric representations and editing modalities,
vertex-based mesh editing through direct manipulation of mesh vertices is
particularly appealing for many applications and this constitutes the main
focus of this work. Classical direct mesh editing methods [BS07, YLW∗21]
employ an optimization framework where the user vertex constraints are
complemented by a regularizer whose main goal is to keep the rest of the shape
realistic. To accomplish this, regularizers either use elasticity or geometric
priors. These regularizers try to minimize locally computed energies and they
often fail for large deformations because: (1) they assume the deformation
behavior is homogeneous across the object, which is not true in practice
[SZGP05, PJS07] and (2), especially for CAD models, there are global semantic
relationships between different parts that are lost when only local
regularizers are considered. Methods such as iWires [GSMCO09] attempted to
address this issue by first analyzing the shape to extract global
relationships between parts using a network of curves, but this is not a
general solution.

One way to address these shortcomings is to use data-driven methods. Data-
driven methods use three main strategies:
(1) learn the shape space of certain classes of objects and use it to guide
the deformations [YM16, WCMN19, YAK∗20, JTM∗21, DYT21, LZZ∗21, YZX∗21, SHM∗22,
HHK∗23]
(2) learn the deformation behavior from a set of example deformations in a
supervised manner [SZGP05, PJS07, TMW∗22]. The biggest challenge in these
first two strategies is the reliance on real 3D data. While the collection of
3D datasets available for research has greatly improved in the last few years,
it is minuscule compared to the richness of 2D images available; it can be
argued that a large gap between these two will always exist due to the
inherent challenges in 3D acquisition compared to 2D acquisition.
(3) Another strategy is to rely on optimizing the shape using 2D priors
obtained from large pre-trained 2D models such as CLIP [MKXBP22, GAG∗23,
SFL∗23] or stable diffusion [PJBM22, MPCOMA23, TMT∗23], and guide the
deformation through differential rendering. This strategy has only been
relatively recently explored in the context of 3D shape synthesis from a text
prompt, but presently none of these can accommodate user-specified geometric
constraints, and it is not obvious how to extend these methods to do so. As a
result, in this category, there is very little work in the context of surface
editing.

In this work, we propose a mesh editing method that allows for direct editing
of vertices and combines a local geometric regularizer with global guidance
based on a large-scale image generative model. This novel formulation supports
realistic shape deformation of dense meshes through vertex level editing of a
small number of vertices. More specifically, given a number of vertex
constraints, a region of influence in the model, and a very brief text
description of the object (could just be the object name, like chair, car,
etc.), we optimize deformation for the 3D mesh using the neural Jacobian field
for mesh representation and a loss which combines user supplied vertex
constraints and ARAP [SA07] for geometric rigidity, and a DDS score obtained
from multiple renderings of the object scored against large scale 2D priors to
provide the global context. Geometry changes are restricted to the region of
influence. We show that our method produces realistic and meaningful
deformations with just a few user constraints yielding better results than
traditional methods. Further, it is not restricted to any specific class or
classes of geometric shapes.

Our main contribution is a direct 3D mesh editing algorithm that yields a
global context-aware realistically deformed mesh by dragging just a few mesh
vertices. More specific contributions are:

  * •

We are the first to integrate the requirement of satisfying user-specified
geometric constraints while utilizing the global context obtained through the
use of 2D priors from a large-scale generative model.

  * •

A novel loss function defined as a weighted combination of user-specified
vertex constraints, ARAP, and a DDS score using the 2D prior, collectively
optimizes for geometry and realism in the deformed shape.

  * •

We introduce a dynamic weighting strategy for the optimization that gradually
accommodates the user constraint in a global context-aware manner.

##  2 Related Work

Figure 1:  Overview of our method. Input: The user specifies constraints (red
vertices moved to new positions in blue). Our approach combines these
constraints with a shape regularizer and DDS score (natural image prior)
applied on multiple views to perform the mesh deformation. Output: resulting
deformation

In basic terms, geometric editing techniques can be divided into data-driven
and non-data-driven methods. Non-data-driven shape deformation techniques
typically frame the editing operation as an optimization problem [YLW∗21]
where user constraints are coupled with a regularizer that guides the rest of
the shape. Shape regularizers come in two flavors: physically-based
regularizers that try to minimize physically-based energy functions, often
based on elastic energies [CZW23] and geometric regularizers that try to
preserve the original shape of the object and are often based on shape
differential operators [BS07, SA07, CPSS10].

With the introduction of neural representations of geometry Yang et al.
[YBHK21] proposed a neural implicit level set representation that supports
editing operations, but the edited shapes are obtained using an optimization
framework with physically-based losses and not data-driven constraints. Yuan
et al. [YSL∗22, YSL∗23] proposed a geometric editing method based on the NERF
representation. However, the editing is performed using the classical as-
rigid-as-possible (ARAP) method applied on a mesh obtained from the NERF and
which is converted back to the NERF representation afterwards.

Regardless of representation, one of the shortcomings of these methods is the
lack of a high-level semantic understanding of the object; so they often
result in unrealistic deformations. To address such shortcomings, some papers
rely on additional heuristics: [KSSCO08, GSMCO09] such as maintaining the
proportion of geometric features taking inspiration from CAD design or say,
editing the object using a set of salient curves on the object. However, a
more general solution is to use a data-driven mechanism to guide the
deformation.

###  2.1 Data-driven Geometric Editing

Data-driven geometric editing methods are tailored to the representation of
the underlying geometry. Implicit neural representations such as DeepSDF
[PFS∗19] or NERF [MST∗21] are popular due to their generative prowess. Yumar
et al. [YM16] proposed a method that uses an implicit occupancy grid to learn
the space of shapes from a set of 3D data. Deng et al. [DYT21] propose a
method that uses a neural signed distance function representation that also
learns a neural shape model from a collection of 3D objects. Both of these
methods are limited to a small and specific set of objects. Liu et al.
[LZZ∗21] introduced a method that learns a conditional radiance field over an
entire object class to guide the deformation behavior. However, all the above
methods are subject to the same limitation of being reliant on the
availability of relevant 3D data which is much less accessible as compared to
2D images.

One way to overcome this limitation is to leverage large-scale image-based
models such as CLIP [RKH∗21] or stable diffusion [PJBM22]. Since models like
CLIP and stable diffusion use text prompts as a controlling mechanism, it has
opened the door to text-based 3D generation methods [PJBM22, MKXBP22]. For
purposes of shape editing, Hyung et al. [HHK∗23] propose a purely text-based
editing framework using a NERF representation and Mikaeili et al. [MPCOMA23]
propose a sketch-based editing method using a NERF representation. The latter
method is close to ours in spirit in that it uses the SDS loss [PJBM22] to
guide the deformation using a sketch-based interface. However, this method
does not allow for direct control over the geometry, an operation that is very
challenging when using an implicit representation and would also necessitate
user interaction in several views.

Despite the popularity of the neural implicit representations mentioned above,
3D triangular meshes are still very much in use in many real life applications
due to their simplicity, efficiency, and downstream hardware processing
support through GPUs. Wang et al. [WCMN19] propose a method that trains an
end-to-end network that deforms the source model to resemble the target.
Because the method infers per-vertex offset displacements, it is not suited
for vertex-based mesh editing say, by specifying only a few vertex
constraints. Wang et al. [YAK∗20] propose a neural cage network that infers
cage coordinates of the points inside. Both these networks are trained by
combining shape-based Laplacian losses and other heuristics tailored towards
generic man-made objects. Therefore despite being neural methods, the
deformation behaviour learned is not driven directly by 3D or deformation
data.

Early data-driven methods [SZGP05, PJS07] focus on learning a deformation
behavior from a set of sample poses in terms of deformation gradients. More
recently Tang et al. [TMW∗22] used supervision to learn the prior deformation
of a specific class of objects, mostly animals from an existing database. All
of the above approaches rely not only on the availability of 3D data, but also
the availability of sample deformations of the same mesh for supervised
learning. Jakab et al. [JTM∗21] discover key feature points in a dataset of
objects and cast the problem as transforming a source 3D object to a target 3D
object from the same object category. The feature points used for deformation
are not user selected. Furthermore, these mesh based methods have similar
limitations to their neural counterparts in that they rely on a 3D database of
objects of a certain class.

To overcome this limitation Gao et al. [GAG∗23] edit a triangular mesh using
only a text prompt using CLIP embeddings. One major challenge with the CLIP
methods is that CLIP embeddings do not encode a viewing direction resulting in
ambiguities denoted as the Janus effect [GAG∗23]. The other major challenge is
that it seems difficult to accommodate the satisfaction of user-specified
geometric constraints.

The use of large scale pre-trained models for 2D image editing has been
explored in several works. Pan et al. [PTL∗23] propose a method for point-
based image editing by optimizing in the latent space of the StyleGAN2
generator [KLA∗20]. Similarly, Shi et al. [SXP∗23] and Mou et al. [MWS∗23]
achieve something similar by optimizing the diffusion latent space. In
[GSW∗22], 3D textured meshes are generated from the learned latent space of
images. Unfortunately, these methods cannot be extended to direct vertex
editing of a given 3D mesh, as they are completely reliant on locating the
deformed representation in the rich latent space of 2D images. As already
mentioned, we do not have that luxury of such a rich latent space for 3D
shapes. In our work, we also use the latent space of images, but we use it
regularize the 3D geometry by enforcing it to result in a natural image when
rendered from an arbitrary view.

In summary, our main goal in this work is to provide a general mesh editing
method with user-specified vertex constraints that do not involve supervision,
hence does not rely on 3D data for training, and produces realistic results
without being restricted to a specific class of objects. We achieve this by
harnessing the rich and vast knowledge about natural and human-made objects
that are represented in today’s pre-trained large-scale models in the image
format. Our method is described in the following section.

Figure 2:  Gallery of results 1 Figure 3:  Gallery of results 2

##  3 Method

###  3.1 Overview

The overview of our proposed method is shown in Figure 1. The user specifies:
a set of mesh vertices (red) as handles to drag, paired with a corresponding
set of target 3D positions (blue), an optional mask to specify the part of the
mesh that is allowed to be modified, and a short (typically one word) prompt
generally describing the object (e.g., a car, a chair, etc.).

We optimize the modifiable part of the mesh in the gradient domain [GAG∗23]
guided by three losses: (1) $l^{2}$ distance of the user constraint loss, (2)
Delta Denoising Score (DDS) loss which makes the deformed mesh have realistic
appearances when rendered from random viewpoints and provides global guidance
to our 3D model, and (3) ARAP loss to control the local geometric behavior.

The remainder of the paper is organized as follows. Sections 3.2\- 3.4 review
the neural jacobian fields, the space in which we do the optimization, the
Delta Denoising Score, and the As-Rigid-as-possible regularizer respectively.
Section 3.5 explains how everything fits together and provides additional
implementation details. Section 4 presents our experiments and analysis of the
method including ablation studies related to our design decisions, and
finally, Section 5 has conclusions, limitations, and future work.

###  3.2 Neural Jacobian field

The neural Jacobian field [AGK∗22] operates in the intrinsic gradient domain
of triangular meshes to learn highly-accurate mappings between meshes. For
triangular mesh vertices $\Phi$, the per-face Jacobian $J_{i}\in\Re^{3\times
3}$ is defined as

$J_{i}=\Phi\nabla^{T}_{i},$ (1)

where $\nabla^{T}_{i}$ is the gradient operator of triangle $t_{i}$. Given
matrices $M_{i}\in\Re^{3\times 3}$ for every triangle, we can compute vertex
positions $\Phi^{*}$ whose Jacobians $J_{i}=\Phi^{*}\nabla^{T}_{i}$ are least-
square closest to $M_{i}$ by solving the Poisson equation. The solution is
obtained by solving the following linear system

$\Phi^{*}=\mathcal{L}^{-1}\mathcal{A}\nabla^{T}M,$ (2)

where $\mathcal{A}$ is the mesh’s mass matrix, $\mathcal{L}$ is the mesh’s
Laplacian, and $M$ is the stacking of the input matrices $M_{i}$. (For a
detailed mathematical definition of the Jacobian field, please refer to
[AGK∗22]). In another result of interest to this work[GAG∗23], Gao et.al found
that when deforming a triangular mesh with guidance from text-to-image models,
such as CLIP[RKH∗21], instead of working directly with vertex positions,
optimizing the Jacobian field of the mesh can produce a more smoothly-
deforming mesh and avoid entanglement.

###  3.3 Delta Denoising Score

Score Distillation Sampling(SDS) was introduced in [PJBM22] to use a $2D$
prior to guide the synthesis of $3D$ shapes. Delta Denoising
Score(DDS)[HACO23] extended the $SDS$ mechanism for guidance in image editing.
Given an input image $x$, a text embedding $z$, a denoising model
$\epsilon_{\phi}$ with parameters $\phi$, a randomly sampled time step
$t\sim\mathcal{U}(0,1)$ drawn from the uniform distribution, and noise
$\epsilon\sim\mathcal{N}(0,I)$ following a normal distribution, the weighted
denoising score can be expressed as

$\mathcal{L}_{Diff}(\phi,x,z,\epsilon,t)=w(t)\|\epsilon_{\phi}(x_{t},z,t)-\epsilon\|^{2}_{2}$
(3)

where $w(t)$ is a weighting function that depends on the time step $t$, and
$x_{t}$ is the $x$ added the noise at time step $t$. For text-conditioned
diffusion models that use classifier-free guidance [HS22], the denoised image
is expressed as

$\hat{\epsilon}_{\phi}(x_{t},z,t)=(1+w)\epsilon_{\phi}(x_{t},z,t)-w\epsilon_{\phi}(x_{t},t),$
(4)

which is a weighted sum of conditioned and unconditioned denoising.

Thus, as shown in [PJBM22], given an arbitrary differentiable parametric
function $g_{\theta}$ that renders images, the gradient of the $SDS$ guidance
is given by

$\nabla_{\theta}\mathcal{L}_{SDS}(x,z,\epsilon,t)=\hat{\epsilon}_{\phi}((x_{t},z,t)-\epsilon){{\partial
x_{t}}\over{\partial\theta}}.$ (5)

Using that gradient to optimize the $g_{\theta}$ can produce images that look
natural. However, for image editing, [HACO23] has shown that $SDS$ can produce
non-detailed and blurry outputs due to the noisy gradient. To overcome this
problem, $DDS$ was introduced which does two $SDS$ processes for the edited
image and reference image respectively, and retrieves the gradient by
subtraction of these two,

$\nabla_{\theta}\mathcal{L}_{DDS}=\nabla_{\theta}\mathcal{L}_{SDS}(x_{edit},z_{edit})-\nabla_{\theta}\mathcal{L}_{SDS}(x_{ref},z_{ref}),$
(6)

where $x_{edit}$ and $x_{ref}$ are the edited image and reference image,
$z_{edit}$ and $z_{ref}$ are the text prompts of the edited image and
reference image. It shows that the $DDS$ has a less noisy gradient and can
produce a higher fidelity image in the image editing task.

###  3.4 As-Rigid-As-Possible Regularizer

The ARAP regularizer [SA07] does not penalize any isometric deformations
allowing local rotations, but it penalizes local stretches. More specifically:

$\mathcal{L}_{reg}=\sum_{i}\sum_{j\in
N(i)}w_{ij}\|(V^{{}^{\prime}}_{i}-V^{{}^{\prime}}_{j})-R_{i}(V_{i}-V_{j})\|^{2},$
(7)

where $V_{i}$ is the initial vertex position, $V^{{}^{\prime}}_{i}$ is the
deformed vertex position, ${N(i)}$ is the set of vertices adjacent to
${V_{i}}$, and ${R_{i}}$ is the local estimation rotation matrix for the one
ring of vertices around vertex ${i}$. ${w_{ij}}$ is the standard cotangent
Laplacian weight [MDSB03]. At every iteration, ${R_{i}}$ can be computed
analytically using SVD decomposition on the local co-variance matrix[Ume91].

Figure 4: Comparison with ARAP shows the more natural results of DragD3D
compared to ARAP: symmetric deformation of the car even when only one anchor
is used, detail preservation of the cat’s face, a more natural extension of
the back of the chair and a more natural deformation of the elephant’s trunk.

###  3.5 Putting it all together: DragD3D

Suppose a user constraint start point is $c_{i}$, and target point is
$c_{i}^{{}^{\prime}}$, then the user constraint loss is given by

$\mathcal{L}_{user}=\sum_{i=0}^{N}\|c_{i}^{{}^{\prime}}-c_{i}\|_{2}^{2},$ (8)

where $N$ is the number of user-specified handles. The $DDS$ gradient is
computed as per equation 6. In our case, suppose the undeformed mesh is
$\Phi$, the deformed mesh is $\hat{\Phi}$, the random viewpoint is $vp$, the
differentiable renderer is $g(\cdot)$, and the text prompt is $z$, the $DDS$
score of the deformed mesh can be expressed as

$\nabla_{\hat{\Phi}}\mathcal{L}_{DDS}=\nabla_{\hat{\Phi}}\mathcal{L}_{SDS}(g(\hat{\Phi},vp),z)-\nabla_{\Phi}\mathcal{L}_{SDS}(g(\Phi,vp),z).$
(9)

For the denoising model $\epsilon_{\phi}$, we use the Stable
Diffusion[RBL∗22]. We also applied the As-Rigid-As-Possible energy on the
deformed mesh as a regularizer by evaluating the equation 7, to preserve the
local structure of the mesh and improve the smoothness, especially when a mask
is specified.

Overall, we optimize for

$\mathcal{L}_{total}=\lambda_{user}\mathcal{L}_{user}+\lambda_{DDS}\mathcal{L}_{DDS}+\lambda_{reg}\mathcal{L}_{reg},$
(10)

where the weight of the user constraint $\lambda_{user}$ was set to
dynamically increase in a linear fashion from $1$ to $50$ over the whole
optimization process.

For DDS guidance, we use Stable Diffusion 2.1[RBL∗22] as the backbone along
with the $Perp-Neg$ algorithm [AZS∗23]. The guidance scale is $100$, and the
gradient scale is $0.00002$. The weight of the $ARAP$ regularizer can be set
in the range of $0.04$ to $0.2$, depending on the deformation and mesh quality
desired.

For all examples in this paper except the $car$, we optimize for $2000$
iterations. With our prototype implementation, it takes around $40$ minutes on
an Nvidia RTX 3090. However, with code optimization and hyperparameter tuning,
we can expect significant improvements in execution times. Since the
deformation of the car is larger, we optimized for $4000$ iterations. We use
Adan[XZL∗22] as the optimizer with a learning rate set to $0.005$.

As noted in [MKXBP22] such a large-scale optimization problem can easily get
stuck in an undesired local minimum, so in order to minimize this possibility
we do the following to help the optimization process. In every iteration, four
camera positions were sampled randomly with azimuth
$\in[-180^{\circ},180^{\circ}]$, elevation angle $\in[0^{\circ},90^{\circ}]$,
camera distance$\in[d_{0},d_{0}+2]$, where $d_{0}$ is the camera distance of
specifying handles.

##  4 Experiments and Discussion

###  4.1 Results

We demonstrate the effectiveness of our method on a variety of meshes
belonging to different object classes, natural and human-made. This clearly
demonstrates that our editing method is not limited to any class or classes of
objects. A gallery of results is shown in the Figures  DragD3D: Vertex-based
Editing for Realistic Mesh Deformations using 2D Diffusion Priors,2 and 3. In
figures 2 and 3, the left three columns are the input to the pipeline, and the
right two columns are the output of the pipeline. A 360 degree view of our
models as well as time lapse video of the optimization is provided in the
accompanying video. All the meshes were rendered with and without texture to
show the geometry clearly. The meshes were obtained from TurboSquid,
Thingi10K[ZJ16], and TEXTure[RMA∗23]. For the meshes without texture (car,
cat, chair, spaceship, castle, and snow monster), we generated the texture
using [RMA∗23]. As can be seen from all result figures, the deformation
results are realistic while satisfying the user constraints and, as shown in
Figure 4, clearly outperform the results from pure geometric deformations that
are not aware of the global structure of the object.

###  4.2 Comparison with optimization-based methods

Our main goal is to provide a general mesh editing method that provides
realistic results by leveraging the large-scale 2D models developed with
millions of images of vast classes of objects from multiple views. Unlike many
others, who are able to do this only for a specific class of objects or a
given set of deformations to provide the guidance to the method. Therefore,
our comparison is with a method that does not use any data-driven guidance for
the deformations. There are a multitude of such methods, we refer to two
state-of-the-art reviews [BS07, YLW∗21]. For practical reasons (i.e.
nonavailability of code, no clear winner in terms of deformation quality),
from this multitude of methods we pick the as rigid as possible as a
representative of the optimization-based class of methods [YLW∗21].

In figure 4, we compare our method with the improved version of As-Rigid-As-
Possible(ARAP)[CPSS10], which was implemented by Libigl[JP∗18]. As noted in
[YLW∗21] the ARAP deformation model is representative of the optimization-
based methods.

Although ARAP is an efficient geometry deformation energy computation, it
introduces some artifacts, e.g. bending, and is not able to deform the mesh
symmetrically when the handle is not symmetric, which makes the deformed shape
look unnatural. For instance, in Figure4, $(a)$, the car roof was only lifted
on the handle side non-smoothly; $(c)$, the nose of the cat was bent and the
mouth was stretched; $(e)$ only the areas around the handles were lifted, the
chair top and the armrest were bent; and $(g)$ the trunk of the elephant was
bent and twisted unnaturally. On the other hand, our method can mitigate
problems, such as bending and tilting, produced by $ARAP$ energy, and as a
result, yielding far more realistic deformations. For example, in Figure4,
$(b)$, the roof of the car was lifted on both sides, and all the other parts
of the car were also deformed to make it look natural; $(d)$, the cat’s nose
and mouth are deformed without bending and stretching, and the left ear isn’t
notched for the silhouette; $(f)$, the top of the chair was deformed to a
headrest which has a small arc and the armrest was kept horizontal; $(h)$,
rather than only deform the trunk as in $(g)$, the head of the elephant was
also lifted, making the deformation look more realistic.

Figure 5: Effectiveness of ARAP regularizer’s weight $\lambda_{reg}$. (a)
$\lambda_{reg}=0.2$, (b) $\lambda_{reg}=0.04$

. Note the geometric artifact in top-b) and the nose and mouth distortion in
bottom-a).

Figure 6: Ablation study. a-f: removing one by one the various features of our
method. g: our method

###  4.3 Design decisions

Weighting of ARAP vs. DDS. One way to understand the effects of these two
losses is to look at them as complementing each other. The DDS loss is very
noisy and it affects the low-frequency part of the shape while the ARAP loss
regularizes better the high-frequency part of the shape. Therefore the ARAP
weight balances the local shape behaviour against the global semantics of the
object and the choice depends on the magnitude of the deformation. Large
deformations have the benefit of a larger influence of the global DDS.
Therefore a smaller ARAP loss works better. Figure 5 illustrates the influence
of the weight. The top result is better with higher ARAP loss Figure 5a) to
keep the geometry locally consistent as small dents like those shown in Figure
5-top-b) are not captured by the noisy DDS loss. Figure 5-bottom shows a
result that works better if we allow a larger influence of the DDS loss
(Figure 5-bottom-b), which translates into a low ARAP weight. A high ARAP loss
will result in unwanted distortions on the face (Figure 5-bottom-a)

Dynamic User Loss Given that the user specified constraint is the more
important one to be satisfied in local mesh editing, the user constraint needs
to be appropriately weighted in our loss function. If a very large weight
value is set for this from the very beginning, then $\mathcal{L}_{user}$ will
dominate the gradient and reduce the effectiveness of the other losses (i.e.
the $ARAP$ regularizer and DDS guidance). This, in turn, will lead to drastic
deformation from the beginning. As shown in Figure 7 $(a)$, in the early
stages of the optimization, this drastic deformation can produce artifacts and
non-smooth mesh which then are difficult to remedy by the $ARAP$ regularizer
and DDS guidance in later optimization stages. The dynamic weighting approach
we have introduced helps mitigate this problem, as can be seen in Figure 7
$(b)$.

Ablation Studies. We perform various experiments to show the effectiveness of
the Jacobian fields, $ARAP$ regularizer, random camera views, guidance,
texture, and dynamic weight strategy. As can be seen in Figure 6, if we
optimize for vertex displacement directly, the mesh structure cannot be kept.
Without $ARAP$ regularizer, the local structure of the connection part between
the mask and unmasked area cannot be maintained, and the mesh is not as smooth
if we do not include the $ARAP$ regularizer. If the camera was fixed on four
canonical views, front, back, and sides, and without any change in the
elevation angle, the deformation was overfitted to these four views, which
produced non-smooth deformation. As shown in Figure 6, if we replace the $DDS$
guidance with the $SDS$ guidance in our pipeline, due to the noisy gradient of
$SDS$, the edited mesh is less smooth and has more distortion. Without
texture, the guidance of the 2D diffusion priors is not as efficient as
before, which makes for lower-quality editing. What’s more, without dynamic
weighting for user constraint, we observed that user constraint loss decreased
so fast that the mesh was deformed in a very distorted way and it made the DDS
guidance incapable of remedying the distortion.

Lastly, we test dependence on the precision of the text prompt given by the
user. As Figure 8 shows, the description does not need to be very precise, as
long as it matches the object. This is clear as the prompt is used only in the
DDS loss and primarily serves to provide the global context for this object’s
shape.

Figure 7: The deformed mesh in the early stage of optimization with and
without the dynamic weighting strategy.

##  5 Conclusions, Limitations and Future work

Dragging 3D mesh vertices in space has been used in practice to provide the
fine control that designers seek in shape design. For dense meshes this poses
the problem of having to automatically determine how all other vertices should
change. Pre-trained large scale image models incorporate vast knowledge about
the appearance of shapes in the real world. Recovering very specific shapes
from this generalized knowledge is very challenging. Tapping into this
knowledge for fine level mesh editing is the challenge we have successfully
addressed in this work.

Figure 8: Same deformation with different prompts. The prompts need not be
very precise.

Our method has a few limitations. As with most mesh based methods, the quality
of the triangulation can affect the result. The run-time of the algorithm is
large compared to traditional mesh deformation methods, but is still aligned
with other optimizations based on similar kinds of losses. In future work, we
will look at ways to speed up the optimization. Our deformation method
requires a simple prompt to accompany the rest of the constraints. This is
easy to provide and in the future, we could use automatic prompting [LRLJ23].
Our method still has some issues dealing with very large deformations. The DDS
loss is very noisy and is not sensitive to small geometric artifacts. These
artifacts are usually corrected by the ARAP loss, but large deformations rely
more on the DDS loss and as a consequence may experience geometric artifacts.

## References

  * [AGK∗22] Aigerman N., Gupta K., Kim V. G., Chaudhuri S., Saito J., Groueix T.:  Neural jacobian fields: Learning intrinsic mappings of arbitrary meshes.  _arXiv preprint arXiv:2205.02904_ (2022).
  * [AZS∗23] Armandpour M., Zheng H., Sadeghian A., Sadeghian A., Zhou M.:  Re-imagine the negative prompt algorithm: Transform 2d diffusion into 3d, alleviate janus problem and beyond.  _arXiv preprint arXiv:2304.04968_ (2023).
  * [BS07] Botsch M., Sorkine O.:  On linear variational surface deformation methods.  _IEEE transactions on visualization and computer graphics 14_ , 1 (2007), 213–230.
  * [CPSS10] Chao I., Pinkall U., Sanan P., Schröder P.:  A simple geometric model for elastic deformations.  _ACM transactions on graphics (TOG) 29_ , 4 (2010), 1–6.
  * [CZW23] Chen H., Zheng C., Wampler K.:  Local deformation for interactive shape editing.  _arXiv preprint arXiv:2306.06550_ (2023).
  * [DYT21] Deng Y., Yang J., Tong X.:  Deformed implicit field: Modeling 3d shapes with learned dense correspondence.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2021), pp. 10286–10296.
  * [GAG∗23] Gao W., Aigerman N., Groueix T., Kim V., Hanocka R.:  Textdeformer: Geometry manipulation using text guidance.  In _ACM SIGGRAPH 2023 Conference Proceedings_ (2023), pp. 1–11.
  * [GSMCO09] Gal R., Sorkine O., Mitra N. J., Cohen-Or D.:  iwires: An analyze-and-edit approach to shape manipulation.  In _ACM SIGGRAPH 2009 papers_. 2009, pp. 1–10.
  * [GSW∗22] Gao J., Shen T., Wang Z., Chen W., Yin K., Li D., Litany O., Gojcic Z., Fidler S.:  Get3d: A generative model of high quality 3d textured shapes learned from images.  In _Advances In Neural Information Processing Systems_ (2022).
  * [HACO23] Hertz A., Aberman K., Cohen-Or D.:  Delta denoising score.  _arXiv preprint arXiv:2304.07090_ (2023).
  * [HHK∗23] Hyung J., Hwang S., Kim D., Lee H., Choo J.:  Local 3d editing via 3d distillation of clip knowledge.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2023), pp. 12674–12684.
  * [HS22] Ho J., Salimans T.:  Classifier-free diffusion guidance.  _arXiv preprint arXiv:2207.12598_ (2022).
  * [JP∗18] Jacobson A., Panozzo D., et al.:  libigl: A simple C++ geometry processing library, 2018.  https://libigl.github.io/.
  * [JTM∗21] Jakab T., Tucker R., Makadia A., Wu J., Snavely N., Kanazawa A.:  Keypointdeformer: Unsupervised 3d keypoint discovery for shape control.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2021), pp. 12783–12792.
  * [KLA∗20] Karras T., Laine S., Aittala M., Hellsten J., Lehtinen J., Aila T.:  Analyzing and improving the image quality of stylegan.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ (2020), pp. 8110–8119.
  * [KSSCO08] Kraevoy V., Sheffer A., Shamir A., Cohen-Or D.:  Non-homogeneous resizing of complex models.  _ACM Transactions on Graphics (TOG) 27_ , 5 (2008), 1–9.
  * [LRLJ23] Luo T., Rockwell C., Lee H., Johnson J.:  Scalable 3d captioning with pretrained models.  _arXiv preprint arXiv:2306.07279_ (2023).
  * [LZZ∗21] Liu S., Zhang X., Zhang Z., Zhang R., Zhu J.-Y., Russell B.:  Editing conditional radiance fields.  In _Proceedings of the IEEE/CVF international conference on computer vision_ (2021), pp. 5773–5783.
  * [MDSB03] Meyer M., Desbrun M., Schröder P., Barr A. H.:  Discrete differential-geometry operators for triangulated 2-manifolds.  In _Visualization and mathematics III_. Springer, 2003, pp. 35–57.
  * [MKXBP22] Mohammad Khalid N., Xie T., Belilovsky E., Popa T.:  Clip-mesh: Generating textured meshes from text using pretrained image-text models.  In _SIGGRAPH Asia 2022 conference papers_ (2022), pp. 1–8.
  * [MPCOMA23] Mikaeili A., Perel O., Cohen-Or D., Mahdavi-Amiri A.:  Sked: Sketch-guided text-based 3d editing.  _arXiv preprint arXiv:2303.10735_ (2023).
  * [MST∗21] Mildenhall B., Srinivasan P. P., Tancik M., Barron J. T., Ramamoorthi R., Ng R.:  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Communications of the ACM 65_ , 1 (2021), 99–106.
  * [MWS∗23] Mou C., Wang X., Song J., Shan Y., Zhang J.:  Dragondiffusion: Enabling drag-style manipulation on diffusion models.  _arXiv preprint arXiv:2307.02421_ (2023).
  * [PFS∗19] Park J. J., Florence P., Straub J., Newcombe R., Lovegrove S.:  Deepsdf: Learning continuous signed distance functions for shape representation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ (2019), pp. 165–174.
  * [PJBM22] Poole B., Jain A., Barron J. T., Mildenhall B.:  Dreamfusion: Text-to-3d using 2d diffusion.  _arXiv preprint arXiv:2209.14988_ (2022).
  * [PJS07] Popa T., Julius D., Sheffer A.:  Interactive and linear material aware deformations.  _International Journal of Shape modeling 13_ , 01 (2007), 73–100.
  * [PTL∗23] Pan X., Tewari A., Leimkühler T., Liu L., Meka A., Theobalt C.:  Drag your gan: Interactive point-based manipulation on the generative image manifold.  In _ACM SIGGRAPH 2023 Conference Proceedings_ (2023), pp. 1–11.
  * [RBL∗22] Rombach R., Blattmann A., Lorenz D., Esser P., Ommer B.:  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ (2022), pp. 10684–10695.
  * [RKH∗21] Radford A., Kim J. W., Hallacy C., Ramesh A., Goh G., Agarwal S., Sastry G., Askell A., Mishkin P., Clark J., et al.:  Learning transferable visual models from natural language supervision.  In _International conference on machine learning_ (2021), PMLR, pp. 8748–8763.
  * [RMA∗23] Richardson E., Metzer G., Alaluf Y., Giryes R., Cohen-Or D.:  Texture: Text-guided texturing of 3d shapes.  _arXiv preprint arXiv:2302.01721_ (2023).
  * [SA07] Sorkine O., Alexa M.:  As-rigid-as-possible surface modeling.  In _Symposium on Geometry processing_ (2007), vol. 4, pp. 109–116.
  * [SFL∗23] Sanghi A., Fu R., Liu V., Willis K. D., Shayani H., Khasahmadi A. H., Sridhar S., Ritchie D.:  Clip-sculptor: Zero-shot generation of high-fidelity and diverse shapes from natural language.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2023), pp. 18339–18348.
  * [SHM∗22] Shechter M., Hanocka R., Metzer G., Giryes R., Cohen-Or D.:  Neuralmls: Geometry-aware control point deformation.  _arXiv preprint arXiv:2201.01873_ (2022).
  * [SXP∗23] Shi Y., Xue C., Pan J., Zhang W., Tan V. Y., Bai S.:  Dragdiffusion: Harnessing diffusion models for interactive point-based image editing.  _arXiv preprint arXiv:2306.14435_ (2023).
  * [SZGP05] Sumner R. W., Zwicker M., Gotsman C., Popović J.:  Mesh-based inverse kinematics.  _ACM transactions on graphics (TOG) 24_ , 3 (2005), 488–495.
  * [TMT∗23] Tsalicoglou C., Manhardt F., Tonioni A., Niemeyer M., Tombari F.:  Textmesh: Generation of realistic 3d meshes from text prompts.  _arXiv preprint arXiv:2304.12439_ (2023).
  * [TMW∗22] Tang J., Markhasin L., Wang B., Thies J., Nießner M.:  Neural shape deformation priors.  _Advances in Neural Information Processing Systems 35_ (2022), 17117–17132.
  * [Ume91] Umeyama S.:  Least-squares estimation of transformation parameters between two point patterns.  _IEEE Transactions on Pattern Analysis & Machine Intelligence 13_, 04 (1991), 376–380.
  * [WCMN19] Wang W., Ceylan D., Mech R., Neumann U.:  3dn: 3d deformation network.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2019), pp. 1038–1046.
  * [XZL∗22] Xie X., Zhou P., Li H., Lin Z., Yan S.:  Adan: Adaptive nesterov momentum algorithm for faster optimizing deep models.  _arXiv preprint arXiv:2208.06677_ (2022).
  * [YAK∗20] Yifan W., Aigerman N., Kim V. G., Chaudhuri S., Sorkine-Hornung O.:  Neural cages for detail-preserving 3d deformations.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2020), pp. 75–83.
  * [YBHK21] Yang G., Belongie S., Hariharan B., Koltun V.:  Geometry processing with neural fields.  _Advances in Neural Information Processing Systems 34_ (2021), 22483–22497.
  * [YLW∗21] Yuan Y.-J., Lai Y.-K., Wu T., Gao L., Liu L.:  A revisit of shape editing techniques: From the geometric to the neural viewpoint.  _Journal of Computer Science and Technology 36_ , 3 (2021), 520–554.
  * [YM16] Yumer M. E., Mitra N. J.:  Learning semantic deformation flows with 3d convolutional networks.  In _European Conference on Computer Vision_ (2016), Springer, pp. 294–311.
  * [YSL∗22] Yuan Y.-J., Sun Y.-T., Lai Y.-K., Ma Y., Jia R., Gao L.:  Nerf-editing: geometry editing of neural radiance fields.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2022), pp. 18353–18364.
  * [YSL∗23] Yuan Y.-J., Sun Y.-T., Lai Y.-K., Ma Y., Jia R., Kobbelt L., Gao L.:  Interactive nerf geometry editing with shape priors.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ (2023).
  * [YZX∗21] Yang B., Zhang Y., Xu Y., Li Y., Zhou H., Bao H., Zhang G., Cui Z.:  Learning object-compositional neural radiance field for editable scene rendering.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ (2021), pp. 13779–13788.
  * [ZJ16] Zhou Q., Jacobson A.:  Thingi10k: A dataset of 10,000 3d-printing models.  _arXiv preprint arXiv:1605.04797_ (2016).