Ryusuke Sugimoto0000-0001-5894-0423University of WaterlooCanadarsugimot@uwaterloo.ca,Nathan King0000-0003-4105-0189University of WaterlooCanadan5king@uwaterloo.ca,Toshiya Hachisuka0000-0003-0284-3776University of WaterlooCanadatoshiya.hachisuka@uwaterloo.caandChristopher Batty0000-0003-3830-7772University of WaterlooCanadachristopher.batty@uwaterloo.ca
Abstract.
We present projected walk on spheres (PWoS), a novel pointwise and discretization-free Monte Carlo solver for surface PDEs with Dirichlet boundaries, as a generalization of the walk on spheres method (WoS)(Sawhney and Crane, 2020; Muller, 1956). We adapt the recursive relationship of WoS designed for PDEs in volumetric domains to a volumetric neighborhood around the surface, and at the end of each recursion step, we project the sample point on the sphere back to the surface. We motivate this simple modification to WoS with the theory of the closest point extension used in the closest point method.To define the valid volumetric neighborhood domain for PWoS, we develop strategies to estimate the local feature size of the surface and to compute the distance to the Dirichlet boundaries on the surface extended in their normal directions. We also design a mean value filtering method for PWoS to improve the method’s efficiency when the surface is represented as a polygonal mesh or a point cloud. Finally, we study the convergence of PWoS and demonstrate its application to graphics tasks, including diffusion curves, geodesic distance computation, and wave propagation animation. We show that our method works with various types of surfaces, including a surface of mixed codimension.
projected walk on spheres, walk on spheres, surface partial differential equations, closest point method, Monte Carlo methods
††ccs: Mathematics of computingPartial differential equations††ccs: Mathematics of computingIntegral equations1. Introduction
Methods to solve surface partial differential equations (PDEs) have become ubiquitous tools in computer graphics research and production. They are used for surface editing(Desbrun etal., 1999), texture synthesis(Turk, 1991), surface fluid animation(Stam, 2003), geodesic distance computation(Crane etal., 2013), and diffusion curves on surfaces(Bartsch etal., 2023; DeGoes etal., 2022) today. A common approach to tackle these problems is to discretize the surface and solve a globally coupled linear system using discrete surface differential operators.
For volumetric PDEs, Monte Carlo methods have recently garnered significant attention in the graphics community due to their unique advantages over traditional discretization-based PDE solvers, including the ability to estimate solution values in a pointwise, spatial discretization-free manner. One such method is the walk on spheres (WoS) method (Muller, 1956) introduced to graphics by Sawhney and Crane (2020). They primarily focused on the (constant coefficient) Poisson and screened Poisson equations in a volumetric domain, and follow-up work likewise emphasizes volumetric problems. In the present work, we consider instead the problem of solving surface PDEs.
Sawhney and Seyb etal.[(2022)] proposed an extension of WoS to support second-order linear elliptic PDEs with spatially varying coefficients, and as one application, they demonstrated a method to solve the Laplace equation on a 2D surface embedded in 3D. However, this approach requires that the conformal parametrization of the surface be readily available, limiting the method’s applicability.
We propose a simpler generalization of WoS for surface PDEs, the projected walk on spheres (PWoS) method, which only assumes the availability of a closest surface point query and an unoriented surface normal direction query. PWoS supports Dirichlet boundary conditions and inherits the advantages of WoS: PWoS is a pointwise, discretization-free Monte Carlo method.Since our method does not require the meshing of the surface, it is particularly advantageous over traditional approaches, such as the finite element method, when the computation can be localized and when the surface is given as an implicit representation, such as a signed distance function.The resulting solution is also free of mesh-dependent discretization artifacts, such as from linear interpolation,as we show in fig.1.Compared to WoS, which performs walks on spheres inside the domain, PWoS performs walks on spheres inside a Cartesian embedding neighborhood domain around the surface. After each step of the walk, it projects the sampled point on the sphere onto the surface.We motivate this simple modification to the original WoS through its connection to the closest point method (CPM) (Ruuth and Merriman, 2008; März and Macdonald, 2012).
Furthermore, inspired by the mean value filtering method for WoS by Bakbouk and Peers (2023), we design a mean value filtering method with a discrete basis function to allow more efficient estimation of solutions when the surface is represented as a polygonal mesh or a point cloud. To confirm PWoS’s accuracy, we perform convergence studies of the method applied to the surface Poisson and screened Poisson equations. Finally, we demonstrate its use in several representative graphics applications, including diffusion curves, geodesic distance computation, and wave propagation animation.
In summary, our key contributions are:
- •
Our novel PWoS algorithm that generalizes the WoS method to solve surface PDEs with Dirichlet boundaries, supported by the theory of the closest point extension.
- •
A mean value filtering method for PWoS with a discrete basis for efficiency improvement.
- •
Evaluation of PWoS with convergence studies and graphics applications.
2. Background
This section briefly reviews the two core mathematical ideas on which our method is based.
2.1. Walk on Spheres
WoS solves volumetric PDEs such as the Poisson equation over a Cartesian domain in .Consider a -ball , centered at with radius , fully contained within the domain. The integral equation
(1) |
holds for the Poisson equation in general, where denotes the surface area of the sphere that bounds the ball and denotes the Green’s function for the Poisson equation on (Sawhney and Crane, 2020). On the right-hand side, the first term is a boundary integral over the -sphere, and the second term is a volume integral over the -ball.If we perform Monte Carlo integration of the first term by uniformly sampling a point on the sphere and of the second term by sampling points inside the ball with probability density function (PDF) , we get the recursive relationship used in WoS:
(2) |
where the hat notation indicates that a variable is a Monte Carlo estimate.The first term on the right-hand side is also a Monte Carlo estimate because the solution is unknown at point in general.At each recursion step, WoS applies this recursive relationship to the largest ball inside the domain bounded by Dirichlet boundaries. It terminates the recursion when the sample point during the recursion falls within a small distance of the domain boundary, by using the known solution at the closest boundary point as the solution estimate: . WoS thereby estimates the solution at each evaluation point independently, offering intrinsic parallelism. Our method generalizes WoS, originally proposed for volumetric PDEs, by incorporating the closest point extension theory of CPM.
2.2. Surface PDEs and Closest Point Extension
Consider the Poisson equation defined on a surface embedded in such that :
(3) |
where denotes the Laplace-Beltrami operator. For convenience, we will use the word surface to refer to any nonzero codimension object in , including mixed-codimension objects. One typically solves such a surface PDE by discretizing the differential operator and solving a sparse linear system. For triangle meshes one can use the cotangent Laplacian(MacNeal, 1949); for other surface representations, a corresponding discrete Laplacian must be defined(Bunge and Botsch, 2023; Sharp and Crane, 2020; Belkin etal., 2009).The closest point method(Ruuth and Merriman, 2008) instead addresses the surface PDE (eq.3) in a more general way by changing the domain to an embedding space, which is a Cartesian neighborhood surrounding the original surface. CPM then solves an embedding PDE defined on the embedding space, whose solution when restricted to points is the solution to the original surface PDE. We briefly summarize CPM theory and refer readers to work by King etal. (2023) for an in-depth review of CPM.
We first assume that is smooth and define the closest point query to the surface, for , as
(4) |
In general, the mapping may not be unique: there may exist more than one closest point for a given . We define the neighborhood where the closest point function is unique as,where is the local feature size at point defined as the minimum Euclidean distance from to the medial axis (Amenta and Bern, 1999). The medial axis is defined as the set of points in where there is more than one closest point. Note that when is a watertight surface the definition of the medial axis that we use contains both the interior part that is bounded by and the exterior part that lies outside the bounded domain.
Within the neighborhood , or a subset of it, surface differential operators can be replaced by Cartesian differential operators with closest point extensions (see (Ruuth and Merriman, 2008; März and Macdonald, 2012)). The closest point extension operator extends surface functions onto to be constant in the normal direction of and is defined as . For functions the extension acts on the restriction of to the surface, i.e., The Laplace-Beltrami operator in eq.3 is equivalent to the following:
(5) |
To define the embedding PDE on , we also extend as . The equation , for , is ill-posed because is constant in the normal direction of but is not guaranteed to be. Therefore, the embedding PDE for eq.3 becomes
(6) |
where is a function that compensates for not being constant in the normal direction of . The function is nonzero for and to ensure eq.6 is consistent with the surface PDE (eq.3) on . Any function with these conditions has the form , where and .
The Macdonald-Brandman-Ruuth approach (see (Chen and Macdonald, 2015, Section 2.3)) takes to allow eq.6 to be written as an equation in one unknown, , since (but importantly except on ). We instead do not restrict to be and interpret as a modification to the source term , then solve for the unknown solution in eq.6. As proven by von Glehn etal. (2013, Section 3.2), since is the extension of a surface function. The property that allows our projection step during the walk in PWoS, detailed in section3.
We show through our numerical examples that taking for all ,provides qualitatively correct results for graphics applications and quantitatively convergent results in most examples in section5. However, the choice of causes eq.6 to be ill-posed as discussed above, and we observe some bias in some convergence studies in section5.1 when is complex. Interesting future work would involve constructing a more accurate function to improve convergence.
In the traditional CPM, one solves the embedding PDE inside a narrow tubular subset of that is within a constant distance to .For the typical grid-based variant, the tubular subset is spatially discretized with a grid of uniform spacing. Interpolation and finite differences are applied on the grid to approximate the closest point extension and the spatial Cartesian differential operators, respectively, and then the resulting linear system is solved. Other variants of CPM(Petras and Ruuth, 2016; Piret, 2012; Cheung etal., 2015) also require some discretization within .Thus, while traditional CPM is agnostic to the specific surface representation, it still discretizes the embedding space and solves a globally coupled system.Moreover, imposing exterior or interior boundary conditions requires tedious grid operations(King etal., 2023).By contrast, we develop a spatial discretization-free, pointwise Monte Carlo estimator for surface PDEs by incorporating the closest point extension concept into WoS.
When there are Dirichlet boundaries , on which the solution is given, one can geometrically extend the boundary itself out into the embedding space in the normal directions, assigning it the same boundary value in accordance with the closest point extension. Note that such boundaries need not coincide with the geometric boundaries of the surface itself.In the context of grid-based CPM, King etal. (2023) discuss how to impose such boundary conditions by duplicating degrees of freedom near the extended boundary.In our work, we devise a method that uses only the closest point function to the (pre-extension) boundary , without the need to construct the extended boundary geometry or perform any complex duplication of degrees of freedom.
3. Method
Input
While our algorithm is generalizable to other configurations, we describe our method for the case when is embedded in and . Recall that we use the word surface to refer both to “surfaces” with (curves) as well as surfaces. We also allow mixed codimension where parts of the surface have and the rest of the surface has . We assume that we can query the closest point function for .Additionally, for surfaces with , we assume that we can query the unoriented normal direction of the surface for . For surfaces with , we assume that we can query the tangential direction of the surface for .These assumptions are valid for common surface representations, including, but not limited to, polygonal meshes, oriented point clouds, and implicit functions. The theory discussed in section2.2 is based on the assumption that is smooth; in practice, we observe that applying our technique on discretized surfaces with sharp features behaves well as we show in section5.1.Additionally, we assume the Dirichlet boundaries have a lower dimension than the dimension of and support the closest point query . When solving a two-sided boundary value problem for boundaries with , we also assume that we can query the tangent direction of .
Overview
The core idea of our method is to apply the WoS recursive relationship within while utilizing the closest point extension constraint that . To do so, we modify the walk process to use spheres contained within and to project the walk position at each recursion step, as illustrated in fig.2.
The problem we solve is the embedding PDE within . The Monte Carlo estimate of eq.2 holds by assuming because the embedding PDE is defined with the Cartesian differential operator.To estimate the surface PDE’s solution at point , we consider a 3D ball centered at and fully contained inside . Theoretically, it should be the largest ball fully contained inside that does not cross the extended Dirichlet boundaries , to minimize the number of steps needed to reach the boundary. We determine the radius of such a ball by taking the minimum of a conservative (under-)estimate of the local feature size at point (section3.1) and the distance to the (extended) Dirichlet boundaries (section3.2).In eq.2, the Monte Carlo estimate of the solution on the sphere, , needs to be evaluated at point , which does not lie on in general.The closest point extension constraint of provides a convenient relationship here:the embedding PDE’s solution at point should coincide with the surface PDE’s solution at the projected point, .We can therefore project the point onto at the end of each recursion step hence .After this projection at the end of each step, we continue the recursion.The source term similarly uses the closest point projection for the closest point extension, replacing the recursive relationship of WoS (eq.2) with
(7) |
We choose in our implementation.
Analogous to the original WoS method, we terminate the recursion when the point falls within a distance of the (extended) Dirichlet boundary by taking the boundary value . We provide pseudocode for an instance of our algorithm in algorithm1, where we highlight the difference between our proposed method and WoS. We also provide a visualization of a potential path of our algorithm when is embedded in in fig.2. Notably, PWoS is a generalization of the WoS algorithm: when (i.e., the codimension-zero case), the local feature size is infinite, the distance to the Dirichlet boundary can be computed with the closest point query , and the last closest point projection of has no effect since .When , in addition to the closest point projection at the end of each recursion step, our algorithm utilizes two nontrivial steps: the local feature size estimation using a medial axis point cloud and the computation of the distance to the (extended) Dirichlet boundary. We discuss these in the following two subsections.
Local Feature Size | Avg. Step Count |
---|---|
0.99 | 31.1398 |
0.5 | 47.84 |
0.25 | 128.981 |
0.125 | 462.781 |
0.0625 | 1818.73 |
3.1. Local Feature Size Estimation
To determine the radius of the sphere centered at that is fully contained inside at each recursion step of the walk, we need a conservative lower bound estimate for the local feature size at .One naive approach would be to use a small enough positive constant value for all , similar to the grid-based CPM (Ruuth and Merriman, 2008). This is a valid strategy, but it would often yield a sphere radius smaller than necessary, requiring more recursion steps for the walks to reach a Dirichlet boundary and making the method inefficient. fig.3 illustrates a result of our algorithm on a unit sphere using different (artificial) local feature size estimates. The analytical local feature size of this surface is everywhere. Although using any smaller value would still give consistent results, it significantly increases the average step count for the walks. For more complex shapes, one cannot compute the analytical local feature size in general and using a small constant value in its place is inefficient. This issue motivates our need for a better local feature size estimate.
To estimate the local feature size, we compute a point cloud approximation of the medial axis and estimate the local feature size as the distance from any query point to its nearest point in the medial axis point cloud.One could use any local feature size estimation algorithm and/or medial axis extraction algorithm (see e.g., (Tagliasacchi etal., 2016)), such as one that outputs or uses the medial axis’ connectivity.Since such methods are often comparatively costly, we employed a simple point cloud-based method, which we describe below.
Medial ball extraction
We first densely scatter points inside a ball in having a radius equal to half the length of the diagonal of the bounding box of , so the entire surface is fully enclosed.For each point , we use the closest point query to compute its two opposing normal directions at .Specifically, we normalize the vector to get the first direction and invert its direction to get the second one. We then use the method of Ma etal. (2012) to extract a point cloud that represents the medial axis, as follows. For each side of the surface (i.e., each normal), we start with a large sphere tangent to , whose center necessarily lies on the normal ray. The initial radius of the ball is set to the length of the diagonal of the bounding box of . Then, we iteratively shrink the size of the sphere, moving its center to maintain tangency at the surface point , until the closest surface point from the center of the sphere does not change. This algorithm gives medial balls, i.e., balls with their centers on the medial axis.As the number of initial scattered points increases, the extracted point cloud balls tend to cover the entire medial axis.While Ma etal. (2012) assumed that the surface is represented as an oriented point cloud, we observed that the algorithm works well with other surface representations by adjusting its termination criteria. See the paper by Ma etal. (2012) and our implementation for details.
Scale axis pruning
Directly using the medial ball centers as the medial axis point cloud does not work well in general when contains any noise or artificial sharp corners introduced by the discretization of a smooth surface. We therefore prefer to estimate (only) the stable part of the medial axis, which is not affected by any small perturbation of . A common solution is, therefore, to prune unstable components of the medial axis, which itself remains an active research topic(Tagliasacchi etal., 2016). We take inspiration from the scale axis transform (SAT)(Giesen etal., 2009; Miklos etal., 2010), but design an alternative since SAT considers topology information of the medial axis, which is unnecessary for local feature size estimation. Our alternative is simpler and faster since topology information is omitted.We first scale all the medial balls by a constant factor . Then, for any pair of medial balls and , we remove the smaller ball from the set of medial balls if it is fully contained in the other. That is, if , we remove the ball .
Note that some of the medial balls may have a very large radius before pruning. For example, an exterior medial ball for a surface of a convex shape would have an infinite radius in theory (but our algorithm returns at most the length of the diagonal of the bounding box of ).When such large medial balls are used in our pruning algorithm, they can easily (and undesirably) remove important, stable parts of the medial axis. The original SAT approach was applied only to the interior medial axis of closed surfaces. Therefore, this issue was not observed since the interior medial ball radii are always bounded and proportional to the size of local features of .To address this problem, we consider each pair of tangent balls generated at the same surface point, and replace the larger one with a tangent ball having the radius of the smaller one. In other words, before we prune the medial axis, we shift the medial ball centers of the larger medial ball in the pair and shrink its radius. In fig.2, we visualize the medial ball centers after this shifting operation.
After this pruning, the set of medial ball centers gives the medial axis point cloud we use to estimate the local feature size.Adjusting the scaling parameter allows us to control the pruning strength. Unless otherwise stated, we use the value for our results.
Conservative and nonzero local feature size
As the medial axis is represented as a point cloud and the nearest point distance may give a larger value than the actual local feature size, we multiply by a small constant ( in our implementation) to ensure a conservative estimate of the local feature size.When there are sharp corners in the geometry, the analytical local feature size is zero, and the walk will become stuck. To prevent this problem, when the estimated local feature size is smaller than a positive constant threshold , we return as the local feature size estimate instead. This process essentially rounds sharp corners with rounding radius . The uniform grid size adopted in grid-based CPM has a similar effect. We do not observe any significant error due to this rounding, as we show in section5.1.
3.2. Distance to Extended Dirichlet Boundaries
Dirichlet boundaries are extended in the normal direction of and the solution in the embedding space on this extended boundary is determined by the closest point extension of Dirichlet values on .Therefore, we need to compute the minimum distance to the extended Dirichlet boundaries, and limit the sphere radius in PWoS further if it is less than the local feature size.To determine the distance to the extended Dirichlet boundary from point , we first find the closest point that lies on the boundary before the extension, .The subset of the extended boundary that is extended from takes the shape of a line segment when and a disk when .We set the line segment’s half-length or the disk radius to the local feature size at using the algorithm in section3.1.We can compute the distance from the point to this line segment or disk without explicitly constructing the extended boundary geometry.When , the distance to the extended boundary is given by
(8) |
where and are the surface normal and the local feature size estimate at , respectively.When , the normal direction is not uniquely defined, so we instead use a similar method based on the surface tangent :
(9) |
3.3. Generalizations
3.3.1. Screened Poisson Equation
We have so far considered only the Poisson equation. For WoS, Sawhney and Crane (2020) proposed a generalization of WoS to the screened Poisson equation , where is a positive constant. The embedding PDE for the screened Poisson equation is constructed similarly to eq.6 using closest point extensions(Chen and Macdonald, 2015, Section 2.3). Therefore, similar to WoS(Sawhney and Crane, 2020), we replace the integral equation (eq.1) used in our recursive relationship with
(10) |
where is the Yukawa potential and is a positive number smaller than . To evaluate the first term, instead of multiplying the contribution from the recursion by as suggested by Sawhney and Crane (2020), we use a Russian Roulette strategy: we terminate the path early with probability with zero contribution and otherwise we use the estimate of the solution without multiplying by . As PWoS is a generalization of WoS, we can apply this Russian Roulette strategy to WoS as well.This scheme allows us to terminate the PWoS path early without waiting for it to reach the boundary and without introducing additional bias. It also makes it possible to apply PWoS to problems without Dirichlet boundaries. We can even use this estimator for the screened Poisson equation as a Tikhonov regularization of the Poisson equation without boundaries: solving the screened Poisson equation with small yields an approximate solution to the Poisson equation. Similar regularization ideas appear in multiple contexts[Sabelfeld and Simonov(1994), Section 6.3; Sawhney and Miller etal.(2023)].
3.3.2. Divergence Source Term and Gradient Estimation
Many graphics applications, such as the heat method for geodesic distance computation(Crane etal., 2013) and the projection step of fluid simulation(Foster and Metaxas, 1996; Stam, 1999), give rise to a Poisson equation with a source term expressed as the divergence of a vector field, , where is a vector field defined over the surface.These applications also require the estimation of the gradient of the solution instead of the solution itself. With grid-based CPM, the differential operators defined in the embedding Cartesian domain can be used to solve such problems(Auer etal., 2012; King etal., 2023). For our PWoS, however, we do not use any embedding grid structure, and we do not assume any specific surface representation, so we cannot use such discrete differential operators.
To solve a problem with a divergence of vector field as the source term, Sugimoto etal. (2024) proposed to use integration by parts to convert the volume integral arising from the source term to a form that does not explicitly depend on the divergence of the vector field. This was done in the context of the walk on boundary method(Sugimoto etal., 2023), which is a Monte Carlo volumetric PDE solver based on a different integral equation formulation, and we adapt this approach to (projected) WoS.To estimate the gradient with PWoS, we can use the gradient estimator for WoS(Sawhney and Crane, 2020, Section 3.1), because the solution’s gradient is zero in the surface normal directions due to the closest point extension constraint . The gradient estimator replaces the first step of recursion based on eq.1 with one we can derive by taking the gradient of eq.1.We discuss the details of these estimators in the supplemental note and demonstrate an application for the geodesic distance computation in section5.2.2.
4. Mean Value Filtering with Discrete Basis
The method we described in section3 is suitable for evaluating the solution at a few evaluation points independently. To improve the efficiency of the method for many evaluation points (i.e., mesh vertices), it is often desirable to utilize the spatial consistency of the solution. For WoS, a few such methodshave been proposed, but those approaches are not trivially applicable to our setup, where we have additional closest point extension constraints between the surface and embedding PDEs.Specifically, the boundary value caching approach[Miller and Sawhney etal.(2023)] is based on a boundary integral equation defined in the Cartesian coordinates and is not applicable to surface PDEs. Adaptation of reverse WoS(Qi etal., 2022) or mean value caching(Bakbouk and Peers, 2023) to our setting would require a mapping of a PDF on the surface to a PDF on the spheres used in our walk process, and we found it difficult to design such an algorithm for general surfaces.
Inspired by the filtering method of Bakbouk and Peers (2023), we designed a filtering method that uses discrete basis functions defined over a surface represented as a polygonal mesh or point cloud. When we apply eq.7 at an evaluation point, if is a precomputed Monte Carlo estimate, we do not need to continue the recursion and can simply use the precomputed estimate in its place. In general, we do not have the estimate available at . Thus, we get the estimate by interpolating the estimate of solutions already computed at a discrete set of points by accepting the bias due to interpolation. For example, for a triangle mesh, we use barycentric coordinates as the interpolation basis.This can be interpreted as a PDE-aware smoothing filter of the solution when we consider this process as a weighted average of the solution estimates at nearby evaluation points.
To estimate the solution at all mesh vertices or all points in a point cloud, we first compute an approximate solution to the problem with a low sample count using the method described in section3 and apply this filtering step to get an updated estimate. We can also precompute the filtering weights per evaluation point first and apply the same filter a few times to achieve an even smoother solution without too much additional cost. We can similarly design a gradient estimation filter based on section3.3.2. While this filtering approach utilizes a discrete basis defined over the surface, we do not use any explicit discrete differential operators and do not need to solve a linear system, which is in contrast to the traditional methods based on discrete Laplacians.
The method of Bakbouk and Peers (2023) similarly designed a smoothing filter for volumetric PDEs. Their method can keep the estimate unbiased by assuming that the evaluation points are sampled with a known PDF and that the original estimates are unbiased, but our method introduces bias due to the use of a discrete basis. We nevertheless observe that, within a reasonable runtime budget,our biased filtering method can reduce the error compared to the PWoS algorithm without filtering. We leave the development of an unbiased variant to future work.
5. Results
We implemented PWoS in Houdini 20.0(Side Effects Software, Inc., 2023) without GPU acceleration, using its built-in closest point queries. Our implementation is provided as supplemental material.
5.1. Error Convergence
We conducted error convergence studies of our method using discretized surfaces. In each scene in fig.4, we change the number of sample paths to plot the root mean squared error measured against the analytical or reference solution. We do not apply the mean value filtering method in these studies. The scene setup includes Laplace, Poisson, and screened Poisson equations, and both smooth surfaces and surfaces with sharp corners. In all scenes, we use and except for Laplace equations, for which we use . While the method shows the expected convergence rate of in most examples, we observed a relatively large bias with problems having more complex source term functions, indicating the need for future work to investigate estimating from eq.6 for such problems.
fig.5 compares the convergence of our method for problems with different local feature sizes. We observe that larger local feature sizecorresponds to faster convergence with lower bias.
Mean value filtering.
We run the mean value filtering algorithm on a Laplace equation on a triangulated curved disk surface in fig.6.In this setup, we observe that applying the filter multiple times with a filter constructed with a sufficiently high sample count can significantly reduce the error, even when the initial estimate is computed with a small number of samples. As expected, the filter is more effective when constructed with more samples, but the error decreases slower than the rate , where is the number of samples used to construct the filter. As no recursive estimation is required with the mean value filtering, it significantly improves the efficiency of PWoS.
5.2. Applications
5.2.1. Diffusion curves
Diffusion curves (Orzan etal., 2008) succinctly represent an image as a collection of curves with associated colors. The final image, exhibiting smooth color gradients, is recovered by solving a Laplace equation with the curves dictating boundary conditions. In our application, we solve the surface Laplace equation using PWoS.With our approach, the surface need not have a boundary curve conforming to the discretized mesh, which contrasts with the common approach(DeGoes etal., 2022).fig.7 shows the reconstruction of color at each point on the discretized surface, represented as a quadrilateral mesh and a combination of triangles, polylines, and point clouds.Our method naturally supports two-sided boundaries, with different colors specified on each side of a curve, and surface geometries with mixed-codimension.Additionally, the pointwise nature of PWoS allows it to be applied in a view-dependent manner. For example, given a camera configuration, for antialiasing, we sample points within each pixel to generate rays. We then generate PWoS samples at the ray-surface intersection points. No computational resources are wasted on surface points that are invisible to the camera(fig.1), and we can obtain clean results without relying on a fine discretization of the surface.Boundary integral-based approaches (Bang etal., 2023; Sun etal., 2012) would similarly allow domain discretization-free evaluation of diffusion curves, but they first require a global linear system solve. Moreover, such methods are not applicable to general curved surfaces, and would need to map the results in the 2D domain to the surface via UV coordinates, for example.
5.2.2. Geodesic distance.
Crane etal. (2013) proposed the heat method, which solves two standard surface PDEs in series to compute the geodesic distance from the boundary .The steps are summarized as
- (1)
where ,
- (2)
, and
- (3)
where ,
where is a small positive constant and is the geodesic distance.Step 2 uses the gradient of the solution to the screened Poisson equation found in Step 1. With a discretization-based method, a discrete gradient operator is used to estimate this gradient; in our method, we directly evaluate the gradient of during Step 1 using the method in section3.3.2, without needing itself. We evaluate the gradient at mesh vertices and normalize it to get at mesh vertices. In Step 3, again, we do not rely on a discrete divergence operator to solve the Poisson equation, but instead use the method described in section3.3.2. When our Poisson solver requires the evaluation of at a point, we interpolate from the mesh vertices and (re)normalize it. We can similarly compute the geodesic distance on a surface represented as a point cloud.fig.8 compares our PWoS-based heat method on surfaces represented as polygonal meshes or oriented point clouds against the heat method with grid-based CPM(King etal., 2023) and the exact geodesic distance computed with Geometry Central(Sharp etal., 2019).Our results are consistent with the reference implementation, albeit with minor deviations.
5.2.3. Surface wave animation
Using our screened Poisson equation solver, we can solve some classes of time-dependent problems. We discretize the wave equation in time with implicit Euler to get a screened Poisson equation and solve it with time stepping (fig.9). At each time step, we store the solution at the vertices of the mesh and query the solution from previous frames by interpolating the values. In contrast to grid-based CPM(Auer etal., 2012), our method directly deals with surface geometry without defining an embedding grid.
5.3. Performance
The performance of our method depends on several factors; we report timings for two representative examples. We measured these timings using a workstation with two Intel(R) Xeon(R) Silver 4316 CPUs, each with 20 CPU cores.For the scene in fig.1 bottom left, the image resolution is 640 by 480 and the number of samples per pixel was 1024. The precomputation step, including medial axis computation, took less than 1 minute, and the rest of the main parts of PWoS took 2 hours and 11 minutes. We did not apply mean value filtering.For the scene in fig.7 left, we have 28,119 evaluation points. The medial axis point cloud extraction took 2.4 seconds, the initial solution estimate with 1 sample took 1.3 seconds, and the application of 10 mean value filtering steps with a filter constructed with 128 samples took 13 seconds.Optimizing the implementation with GPU acceleration may further improve performance.
6. Conclusion and Discussion
We have developed a Monte Carlo method for surface PDEs by augmenting the formulation of the walk on spheres method with a closest point projection step. Our algorithm is justified through its connection to the theory of closest point extensions drawn from the CPM literature. To accelerate its convergence, we have developed a practical mean value filtering method that utilizes a discrete basis defined over the surface. We have further analyzed the method’s convergence on representative analytical tests and demonstrated its application to graphics problems.
PWoS currently supports only Dirichlet boundary conditions; efficient Neumann or Robin boundary handling similar to the walk on stars method for volumetric PDEs[Miller and Sawhney etal.(2024); Sawhney and Miller etal.(2023)] would require the availability of a few more queries, such as a ray intersection query against the (extended) boundaries.
While we used a local feature estimation algorithm to allow walks with larger step sizes, the local feature size estimation itself imposes additional smoothness assumptions on the surface. To respect small-scale local features, the walk can require many iterations to reach a Dirichlet boundary. This effect is partly due to our algorithm (like WoS) being based on an integral equation that holds only locally inside a ball. Revisiting this choice using an integral equation based on a global relationship, such as the one underpinning the walk on boundary method(Sugimoto etal., 2023), could lead to a more efficient alternative for surface PDEs.
Lastly, our method relied on the assumption that the closest point extension compensation term (i.e., ineq.6) in the embedding PDE is negligible. We empirically showed that the algorithm designed with this assumption works well when the source term has a relatively simple expression, but we do not yet have a complete understanding of when this assumption is strictly valid.However, since tends to zero continuously as approaches the surface, the influence of ignoring this term is expected to decrease as we shrink the embedding space (i.e., shrink the spheresize). One can always take a smaller sphere size, albeit at a higher computational cost, as we show in fig.10.Extending our method to consider the effect of the compensation term would further improve the reliability and broaden the applicability of our method.
Acknowledgements.
This research was partially funded by NSERC Discovery Grants (RGPIN-2021-02524 & RGPIN-2020-03918), CFI-JELF (Grant 40132), and a grant from Autodesk.The first author was partially funded by the David R. Cheriton Graduate Scholarship.The second author was partially funded by the Ontario Graduate Scholarship.
References
- (1)
- Amenta and Bern (1999)N. Amenta and M. Bern.1999.Surface Reconstruction by Voronoi Filtering.Discrete & Computational Geometry22, 4 (1999),481–504.https://doi.org/10.1007/PL00009475
- Auer etal. (2012)S. Auer, C.B. Macdonald,M. Treib, J. Schneider, andR. Westermann. 2012.Real-Time Fluid Effects on Surfaces using theClosest Point Method.Computer Graphics Forum31, 6 (2012),1909–1923.https://doi.org/10.1111/j.1467-8659.2012.03071.x
- Bakbouk and Peers (2023)Ghada Bakbouk and PieterPeers. 2023.Mean Value Caching for Walk on Spheres. InEurographics Symposium on Rendering,Tobias Ritschel andAndrea Weidlich (Eds.). TheEurographics Association, Delft, Netherlands,10pages.https://doi.org/10.2312/sr.20231120
- Bang etal. (2023)Seungbae Bang, KirillSerkh, Oded Stein, and Alec Jacobson.2023.An Adaptive Fast-Multipole-Accelerated HybridBoundary Integral Equation Method for Accurate Diffusion Curves.ACM Trans. Graph. 42,6, Article 215 (dec2023), 28pages.https://doi.org/10.1145/3618374
- Bartsch etal. (2023)Alec Bartsch, ColinThompson, and Fernando de Goes.2023.A Procedural Approach for Stylized Bark Shading.In ACM SIGGRAPH 2023 Talks (¡conf-loc¿, ¡city¿LosAngeles¡/city¿, ¡state¿CA¡/state¿, ¡country¿USA¡/country¿, ¡/conf-loc¿)(SIGGRAPH ’23). Association forComputing Machinery, New York, NY, USA, Article10, 2pages.https://doi.org/10.1145/3587421.3595419
- Belkin etal. (2009)Mikhail Belkin, Jian Sun,and Yusu Wang. 2009.Constructing Laplace operator from point clouds in. In Proceedings of the TwentiethAnnual ACM-SIAM Symposium on Discrete Algorithms (New York, New York)(SODA ’09). Society forIndustrial and Applied Mathematics, USA,1031–1040.https://doi.org/10.1137/1.9781611973068.112
- Bunge and Botsch (2023)A. Bunge and M.Botsch. 2023.A Survey on Discrete Laplacians for GeneralPolygonal Meshes.Computer Graphics Forum42, 2 (2023),521–544.https://doi.org/10.1111/cgf.14777
- Chen and Macdonald (2015)Yujia Chen and ColinB.Macdonald. 2015.The Closest Point Method and Multigrid Solvers forElliptic Equations on Surfaces.SIAM Journal on Scientific Computing37, 1 (2015),A134–A155.https://doi.org/10.1137/130929497
- Cheung etal. (2015)KaChun Cheung, LeevanLing, and StevenJ. Ruuth.2015.A localized meshless method for diffusion on foldedsurfaces.J. Comput. Phys. 297(2015), 194–206.https://doi.org/10.1016/j.jcp.2015.05.021
- Crane etal. (2013)Keenan Crane, ClarisseWeischedel, and Max Wardetzky.2013.Geodesics in heat: A new approach to computingdistance based on heat flow.ACM Trans. Graph. 32,5, Article 152 (oct2013), 11pages.https://doi.org/10.1145/2516971.2516977
- DeGoes etal. (2022)Fernando DeGoes, WilliamSheffler, and Kurt Fleischer.2022.Character articulation through profile curves.ACM Trans. Graph. 41,4, Article 139 (jul2022), 14pages.https://doi.org/10.1145/3528223.3530060
- Desbrun etal. (1999)Mathieu Desbrun, MarkMeyer, Peter Schröder, and AlanH.Barr. 1999.Implicit fairing of irregular meshes usingdiffusion and curvature flow. In Proceedings ofthe 26th Annual Conference on Computer Graphics and Interactive Techniques(SIGGRAPH ’99). ACMPress/Addison-Wesley Publishing Co., USA,317–324.https://doi.org/10.1145/311535.311576
- Foster and Metaxas (1996)Nick Foster and DimitriMetaxas. 1996.Realistic Animation of Liquids.Graphical Models and Image Processing58, 5 (1996),471–483.https://doi.org/10.1006/gmip.1996.0039
- Giesen etal. (2009)Joachim Giesen, BalintMiklos, Mark Pauly, and CamilleWormser. 2009.The scale axis transform. InProceedings of the Twenty-Fifth Annual Symposium onComputational Geometry (Aarhus, Denmark) (SCG’09). Association for Computing Machinery,New York, NY, USA, 106–115.https://doi.org/10.1145/1542362.1542388
- King etal. (2023)Nathan King, Haozhe Su,Mridul Aanjaneya, Steven Ruuth, andChristopher Batty. 2023.A Closest Point Method for Surface PDEs with InteriorBoundary Conditions for Geometry Processing.arXiv:2305.04711[cs.GR]
- Ma etal. (2012)Jaehwan Ma, SangWon Bae,and Sunghee Choi. 2012.3D medial axis point approximation using nearestneighbors and the normal field.The Visual Computer 28,1 (2012), 7–19.https://doi.org/10.1007/s00371-011-0594-7
- MacNeal (1949)RichardHenri MacNeal.1949.The solution of partial differential equationsby means of electrical networks.Ph. D. Dissertation.California Institute of Technology.https://doi.org/10.7907/PZ04-5290
- März and Macdonald (2012)Thomas März andColinB. Macdonald. 2012.Calculus on Surfaces with General Closest PointFunctions.SIAM J. Numer. Anal. 50,6 (2012), 3303–3328.https://doi.org/10.1137/120865537
- Miklos etal. (2010)Balint Miklos, JoachimGiesen, and Mark Pauly.2010.Discrete scale axis representations for 3Dgeometry. In ACM SIGGRAPH 2010 Papers (LosAngeles, California) (SIGGRAPH ’10).Association for Computing Machinery,New York, NY, USA, Article 101,10pages.https://doi.org/10.1145/1833349.1778838
- Miller etal. (2023)Bailey Miller, RohanSawhney, Keenan Crane, and IoannisGkioulekas. 2023.Boundary Value Caching for Walk on Spheres.ACM Trans. Graph. 42,4, Article 82 (jul2023), 11pages.https://doi.org/10.1145/3592400
- Miller etal. (2024)Bailey Miller, RohanSawhney, Keenan Crane, and IoannisGkioulekas. 2024.Walkin’ Robin: Walk on Stars with Robin BoundaryConditions.ACM Trans. Graph. 43,4, Article 41 (jul2024), 18pages.https://doi.org/10.1145/3658153
- Muller (1956)MervinE. Muller.1956.Some Continuous Monte Carlo Methods for theDirichlet Problem.The Annals of Mathematical Statistics27, 3 (1956),569 – 589.https://doi.org/10.1214/aoms/1177728169
- Orzan etal. (2008)Alexandrina Orzan, AdrienBousseau, Holger Winnemöller, PascalBarla, Joëlle Thollot, and DavidSalesin. 2008.Diffusion curves: a vector representation forsmooth-shaded images. In ACM SIGGRAPH 2008Papers (Los Angeles, California) (SIGGRAPH ’08).Association for Computing Machinery,New York, NY, USA, Article 92,8pages.https://doi.org/10.1145/1399504.1360691
- Petras and Ruuth (2016)A. Petras and S.J.Ruuth. 2016.PDEs on moving surfaces via the closest pointmethod and a modified grid based particle method.J. Comput. Phys. 312(2016), 139–156.https://doi.org/10.1016/j.jcp.2016.02.024
- Piret (2012)Cécile Piret.2012.The orthogonal gradients method: A radial basisfunctions method for solving partial differential equations on arbitrarysurfaces.J. Comput. Phys. 231,14 (2012), 4662–4675.https://doi.org/10.1016/j.jcp.2012.03.007
- Qi etal. (2022)Yang Qi, Dario Seyb,Benedikt Bitterli, and WojciechJarosz. 2022.A bidirectional formulation for Walk on Spheres.Computer Graphics Forum41, 4 (2022),51–62.https://doi.org/10.1111/cgf.14586
- Ruuth and Merriman (2008)StevenJ. Ruuth andBarry Merriman. 2008.A simple embedding method for solving partialdifferential equations on surfaces.J. Comput. Phys. 227,3 (jan 2008),1943–1961.https://doi.org/10.1016/j.jcp.2007.10.009
- Sabelfeld and Simonov (1994)KarlK. Sabelfeld andNikolaiA. Simonov. 1994.Random Walks on Boundary for SolvingPDEs.De Gruyter, Berlin.https://doi.org/10.1515/9783110942026
- Sawhney and Crane (2020)Rohan Sawhney and KeenanCrane. 2020.Monte Carlo geometry processing: a grid-freeapproach to PDE-based methods on volumetric domains.ACM Trans. Graph. 39,4, Article 123 (aug2020), 18pages.https://doi.org/10.1145/3386569.3392374
- Sawhney etal. (2023)Rohan Sawhney, BaileyMiller, Ioannis Gkioulekas, and KeenanCrane. 2023.Walk on Stars: A Grid-Free Monte Carlo Method forPDEs with Neumann Boundary Conditions.ACM Trans. Graph. 42,4 (aug 2023),22pages.https://doi.org/10.1145/3592398
- Sawhney etal. (2022)Rohan Sawhney, DarioSeyb, Wojciech Jarosz, and KeenanCrane. 2022.Grid-Free Monte Carlo for PDEs with SpatiallyVarying Coefficients.ACM Trans. Graph. 41,4, Article 53 (jul2022), 17pages.https://doi.org/10.1145/3528223.3530134
- Sharp and Crane (2020)Nicholas Sharp andKeenan Crane. 2020.A Laplacian for Nonmanifold Triangle Meshes.Computer Graphics Forum39, 5 (2020),69–80.https://doi.org/10.1111/cgf.14069
- Sharp etal. (2019)Nicholas Sharp, KeenanCrane, etal. 2019.GeometryCentral: A modern C++ library ofdata structures and algorithms for geometry processing.https://geometry-central.net/
- Side Effects Software, Inc. (2023)Side Effects Software, Inc.2023.Houdini.https://www.sidefx.com/products/houdini/Computer Software.
- Stam (1999)Jos Stam. 1999.Stable Fluids. InProceedings of the 26th Annual Conference onComputer Graphics and Interactive Techniques(SIGGRAPH ’99). ACMPress/Addison-Wesley Publishing Co., USA,121–128.https://doi.org/10.1145/311535.311548
- Stam (2003)Jos Stam. 2003.Flows on surfaces of arbitrary topology.ACM Trans. Graph. 22,3 (jul 2003),724–731.https://doi.org/10.1145/882262.882338
- Sugimoto etal. (2024)Ryusuke Sugimoto,Christopher Batty, and ToshiyaHachisuka. 2024.Velocity-Based Monte Carlo Fluids. InACM SIGGRAPH 2024 Conference Papers (Denver, CO,USA) (SIGGRAPH ’24). Associationfor Computing Machinery, New York, NY, USA, Article8, 11pages.https://doi.org/10.1145/3641519.3657405
- Sugimoto etal. (2023)Ryusuke Sugimoto, TerryChen, Yiti Jiang, Christopher Batty,and Toshiya Hachisuka. 2023.A Practical Walk-on-Boundary Method for BoundaryValue Problems.ACM Trans. Graph. 42,4, Article 81 (jul2023), 16pages.https://doi.org/10.1145/3592109
- Sun etal. (2012)Xin Sun, Guofu Xie,Yue Dong, Stephen Lin,Weiwei Xu, Wencheng Wang,Xin Tong, and Baining Guo.2012.Diffusion curve textures for resolution independenttexture mapping.ACM Trans. Graph. 31,4, Article 74 (jul2012), 9pages.https://doi.org/10.1145/2185520.2185570
- Tagliasacchi etal. (2016)Andrea Tagliasacchi,Thomas Delame, Michela Spagnuolo,Nina Amenta, and Alexandru Telea.2016.3D Skeletons: A State-of-the-Art Report.Computer Graphics Forum35, 2 (2016),573–597.https://doi.org/10.1111/cgf.12865
- Turk (1991)Greg Turk.1991.Generating textures on arbitrary surfaces usingreaction-diffusion.SIGGRAPH Comput. Graph.25, 4 (jul1991), 289–298.https://doi.org/10.1145/127719.122749
- von Glehn etal. (2013)Ingrid von Glehn, ThomasMärz, and ColinB. Macdonald.2013.An embedded method-of-lines approach to solvingpartial differential equations on surfaces.arXiv:1307.5657[math.NA]