Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (2024)

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 equationsccs: Mathematics of computingIntegral equations

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (1)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (2)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (3)

1. 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 dsuperscript𝑑{\mathbb{R}^{d}}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT.Consider a d𝑑ditalic_d-ball Br(𝐱)subscript𝐵𝑟𝐱B_{r}(\mathbf{x})italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ), centered at 𝐱𝐱\mathbf{x}bold_x with radius r𝑟ritalic_r, fully contained within the domain. The integral equation

(1)u(𝐱)=1|Br(𝐱)|Br(𝐱)u(𝐲)d𝐲+Br(𝐱)f(𝐳)G(𝐱,𝐳)d𝐳𝑢𝐱1subscript𝐵𝑟𝐱subscriptsubscript𝐵𝑟𝐱𝑢𝐲differential-d𝐲subscriptsubscript𝐵𝑟𝐱𝑓𝐳𝐺𝐱𝐳differential-d𝐳u(\mathbf{x})=\frac{1}{\lvert\partial B_{r}(\mathbf{x})\rvert}\int_{\partial B%_{r}(\mathbf{x})}u(\mathbf{y})\,\mathrm{d}\mathbf{y}+\int_{B_{r}(\mathbf{x})}f%(\mathbf{z})G(\mathbf{x},\mathbf{z})\,\mathrm{d}\mathbf{z}italic_u ( bold_x ) = divide start_ARG 1 end_ARG start_ARG | ∂ italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) | end_ARG ∫ start_POSTSUBSCRIPT ∂ italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) end_POSTSUBSCRIPT italic_u ( bold_y ) roman_d bold_y + ∫ start_POSTSUBSCRIPT italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) end_POSTSUBSCRIPT italic_f ( bold_z ) italic_G ( bold_x , bold_z ) roman_d bold_z

holds for the Poisson equation Δu=fΔ𝑢𝑓\Delta u=froman_Δ italic_u = italic_f in general, where |Br(𝐱)|subscript𝐵𝑟𝐱\lvert\partial B_{r}(\mathbf{x})\rvert| ∂ italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) | denotes the surface area of the sphere that bounds the ball Br(𝐱)subscript𝐵𝑟𝐱B_{r}(\mathbf{x})italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) and G𝐺Gitalic_G denotes the Green’s function for the Poisson equation on Br(𝐱)subscript𝐵𝑟𝐱B_{r}(\mathbf{x})italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x )(Sawhney and Crane, 2020). On the right-hand side, the first term is a boundary integral over the (d1)𝑑1(d-1)( italic_d - 1 )-sphere, and the second term is a volume integral over the d𝑑ditalic_d-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 NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT points 𝐳isubscript𝐳𝑖\mathbf{z}_{i}bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT inside the ball with probability density function (PDF) p(𝐳i)𝑝subscript𝐳𝑖p(\mathbf{z}_{i})italic_p ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), we get the recursive relationship used in WoS:

(2)u^(𝐱)=u^(𝐲)+1NVi=1NVG(𝐱,𝐳i)f(𝐳i)p(𝐳i),^𝑢𝐱^𝑢𝐲1subscript𝑁𝑉superscriptsubscript𝑖1subscript𝑁𝑉𝐺𝐱subscript𝐳𝑖𝑓subscript𝐳𝑖𝑝subscript𝐳𝑖\widehat{u}(\mathbf{x})=\widehat{u}(\mathbf{y})+\frac{1}{N_{V}}\sum_{i=1}^{N_{%V}}\frac{G(\mathbf{x},\mathbf{z}_{i})f(\mathbf{z}_{i})}{p(\mathbf{z}_{i})},over^ start_ARG italic_u end_ARG ( bold_x ) = over^ start_ARG italic_u end_ARG ( bold_y ) + divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG italic_G ( bold_x , bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_f ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_ARG start_ARG italic_p ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_ARG ,

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 u(𝐲)𝑢𝐲u(\mathbf{y})italic_u ( bold_y ) is unknown at point 𝐲𝐲\mathbf{y}bold_y 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 𝐱𝐱\mathbf{x}bold_x during the recursion falls within a small distance ϵitalic-ϵ\epsilonitalic_ϵ of the domain boundary, by using the known solution at the closest boundary point 𝐱¯¯𝐱\overline{\mathbf{x}}over¯ start_ARG bold_x end_ARG as the solution estimate: u^(𝐱)=u(𝐱¯)^𝑢𝐱𝑢¯𝐱\widehat{u}(\mathbf{x})=u(\overline{\mathbf{x}})over^ start_ARG italic_u end_ARG ( bold_x ) = italic_u ( over¯ start_ARG bold_x end_ARG ). 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 𝒮𝒮\mathcal{S}caligraphic_S embedded in dsuperscript𝑑{\mathbb{R}^{d}}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT such that dim(𝒮)<ddimension𝒮𝑑\dim(\mathcal{S})<droman_dim ( caligraphic_S ) < italic_d:

(3)Δ𝒮u𝒮(𝐲)=f𝒮(𝐲),𝐲𝒮,formulae-sequencesubscriptΔ𝒮subscript𝑢𝒮𝐲subscript𝑓𝒮𝐲𝐲𝒮\Delta_{\mathcal{S}}u_{\mathcal{S}}(\mathbf{y})=f_{\mathcal{S}}(\mathbf{y}),%\quad\mathbf{y}\in\mathcal{S},roman_Δ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) = italic_f start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) , bold_y ∈ caligraphic_S ,

where Δ𝒮subscriptΔ𝒮\Delta_{\mathcal{S}}roman_Δ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT denotes the Laplace-Beltrami operator. For convenience, we will use the word surface to refer to any nonzero codimension object in dsuperscript𝑑{\mathbb{R}^{d}}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT, 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 𝐲𝒮𝐲𝒮\mathbf{y}\in\mathcal{S}bold_y ∈ caligraphic_S is the solution u𝒮(𝐲)subscript𝑢𝒮𝐲u_{\mathcal{S}}(\mathbf{y})italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) 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 𝒮𝒮\mathcal{S}caligraphic_S is smooth and define the closest point query to the surface, for 𝐱d𝐱superscript𝑑\mathbf{x}\in{\mathbb{R}^{d}}bold_x ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT, as

(4)cp𝒮(𝐱)=argmin𝐲𝒮𝐱𝐲2.subscriptcp𝒮𝐱𝐲𝒮argminsubscriptdelimited-∥∥𝐱𝐲2\mathrm{cp}_{\mathcal{S}}(\mathbf{x})=\underset{\mathbf{y}\in\mathcal{S}}{%\operatorname{argmin}}\lVert\mathbf{x}-\mathbf{y}\rVert_{2}.roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) = start_UNDERACCENT bold_y ∈ caligraphic_S end_UNDERACCENT start_ARG roman_argmin end_ARG ∥ bold_x - bold_y ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT .

In general, the mapping cp𝒮(𝐱)subscriptcp𝒮𝐱\mathrm{cp}_{\mathcal{S}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) may not be unique: there may exist more than one closest point for a given 𝐱𝐱\mathbf{x}bold_x. We define the neighborhood 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) where the closest point function is unique as𝒩(𝒮)={𝐱d|𝐱cp𝒮(𝐱)2<LFS(cp𝒮(𝐱))}𝒩𝒮conditional-set𝐱superscript𝑑subscriptdelimited-∥∥𝐱subscriptcp𝒮𝐱2LFSsubscriptcp𝒮𝐱\mathcal{N}(\mathcal{S})=\left\{\mathbf{x}\in{\mathbb{R}^{d}}\;\big{|}\;\lVert%\mathbf{x}-\mathrm{cp}_{\mathcal{S}}(\mathbf{x})\rVert_{2}<\mathrm{LFS}\left(%\mathrm{cp}_{\mathcal{S}}(\mathbf{x})\right)\right\}caligraphic_N ( caligraphic_S ) = { bold_x ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT | ∥ bold_x - roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT < roman_LFS ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ) },where LFS(𝐲)LFS𝐲\mathrm{LFS}(\mathbf{y})roman_LFS ( bold_y ) is the local feature size at point 𝐲𝒮𝐲𝒮\mathbf{y}\in\mathcal{S}bold_y ∈ caligraphic_S defined as the minimum Euclidean distance from 𝐲𝐲\mathbf{y}bold_y to the medial axis med(𝒮)med𝒮\mathrm{med}(\mathcal{S})roman_med ( caligraphic_S )(Amenta and Bern, 1999). The medial axis med(𝒮)med𝒮\mathrm{med}(\mathcal{S})roman_med ( caligraphic_S ) is defined as the set of points in dsuperscript𝑑{\mathbb{R}^{d}}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT where there is more than one closest point. Note that when 𝒮𝒮\mathcal{S}caligraphic_S is a watertight surface the definition of the medial axis that we use contains both the interior part that is bounded by 𝒮𝒮\mathcal{S}caligraphic_S and the exterior part that lies outside the bounded domain.

Within the neighborhood 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ), 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 E𝐸Eitalic_E extends surface functions onto 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) to be constant in the normal direction of 𝒮𝒮\mathcal{S}caligraphic_S and is defined as Eu𝒮(𝐱)=u𝒮(cp𝒮(𝐱))𝐸subscript𝑢𝒮𝐱subscript𝑢𝒮subscriptcp𝒮𝐱Eu_{\mathcal{S}}(\mathbf{x})=u_{\mathcal{S}}(\mathrm{cp}_{\mathcal{S}}(\mathbf%{x}))italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) = italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ). For functions u𝒩(𝒮)𝑢𝒩𝒮u\in\mathcal{N}(\mathcal{S})italic_u ∈ caligraphic_N ( caligraphic_S ) the extension E𝐸Eitalic_E acts on the restriction of u𝑢uitalic_u to the surface, i.e., Eu=E(u|𝒮).𝐸𝑢𝐸evaluated-at𝑢𝒮Eu=E(u|_{\mathcal{S}}).italic_E italic_u = italic_E ( italic_u | start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ) . The Laplace-Beltrami operator in eq.3 is equivalent to the following:

(5)Δ𝒮u𝒮(𝐲)=Δ[Eu𝒮](𝐲),𝐲𝒮.formulae-sequencesubscriptΔ𝒮subscript𝑢𝒮𝐲Δdelimited-[]𝐸subscript𝑢𝒮𝐲𝐲𝒮\Delta_{\mathcal{S}}u_{\mathcal{S}}(\mathbf{y})=\Delta[Eu_{\mathcal{S}}](%\mathbf{y}),\quad\mathbf{y}\in\mathcal{S}.roman_Δ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) = roman_Δ [ italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ] ( bold_y ) , bold_y ∈ caligraphic_S .

To define the embedding PDE on 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ), we also extend f𝒮subscript𝑓𝒮f_{\mathcal{S}}italic_f start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT as f(𝐱)=Ef𝒮(𝐱)𝑓𝐱𝐸subscript𝑓𝒮𝐱f(\mathbf{x})=Ef_{\mathcal{S}}(\mathbf{x})italic_f ( bold_x ) = italic_E italic_f start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ). The equation Δ[Eu𝒮](𝐱)=f(𝐱)Δdelimited-[]𝐸subscript𝑢𝒮𝐱𝑓𝐱\Delta[Eu_{\mathcal{S}}](\mathbf{x})=f(\mathbf{x})roman_Δ [ italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ] ( bold_x ) = italic_f ( bold_x ), for 𝐱𝒩(𝒮)𝐱𝒩𝒮\mathbf{x}\in\mathcal{N}(\mathcal{S})bold_x ∈ caligraphic_N ( caligraphic_S ), is ill-posed because f𝑓fitalic_f is constant in the normal direction of 𝒮𝒮\mathcal{S}caligraphic_S but Δ[Eu𝒮]Δdelimited-[]𝐸subscript𝑢𝒮\Delta[Eu_{\mathcal{S}}]roman_Δ [ italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ] is not guaranteed to be. Therefore, the embedding PDE for eq.3 becomes

(6)Δ[Eu𝒮](𝐱)=f(𝐱)+g(𝐱),𝐱𝒩(𝒮),formulae-sequenceΔdelimited-[]𝐸subscript𝑢𝒮𝐱𝑓𝐱𝑔𝐱𝐱𝒩𝒮\Delta[Eu_{\mathcal{S}}](\mathbf{x})=f(\mathbf{x})+g(\mathbf{x}),\quad\mathbf{%x}\in\mathcal{N}(\mathcal{S}),roman_Δ [ italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ] ( bold_x ) = italic_f ( bold_x ) + italic_g ( bold_x ) , bold_x ∈ caligraphic_N ( caligraphic_S ) ,

where g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) is a function that compensates for Δ[Eu𝒮]Δdelimited-[]𝐸subscript𝑢𝒮\Delta[Eu_{\mathcal{S}}]roman_Δ [ italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ] not being constant in the normal direction of 𝒮𝒮\mathcal{S}caligraphic_S. The function g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) is nonzero for 𝐱𝒩(𝒮)𝒮𝐱𝒩𝒮𝒮\mathbf{x}\in\mathcal{N}(\mathcal{S})\setminus\mathcal{S}bold_x ∈ caligraphic_N ( caligraphic_S ) ∖ caligraphic_S and g|𝒮=0evaluated-at𝑔𝒮0g|_{\mathcal{S}}=0italic_g | start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = 0 to ensure eq.6 is consistent with the surface PDE (eq.3) on 𝒮𝒮\mathcal{S}caligraphic_S. Any function g𝑔gitalic_g with these conditions has the form g(𝐱)=γ(v(𝐱)Ev(𝐱))𝑔𝐱𝛾𝑣𝐱𝐸𝑣𝐱g(\mathbf{x})=\gamma(v(\mathbf{x})-Ev(\mathbf{x}))italic_g ( bold_x ) = italic_γ ( italic_v ( bold_x ) - italic_E italic_v ( bold_x ) ), where γ𝛾\gamma\in\mathbb{R}italic_γ ∈ blackboard_R and γ0𝛾0\gamma\neq 0italic_γ ≠ 0.

The Macdonald-Brandman-Ruuth approach (see (Chen and Macdonald, 2015, Section 2.3)) takes v|𝒮=u𝒮evaluated-at𝑣𝒮subscript𝑢𝒮v|_{\mathcal{S}}=u_{\mathcal{S}}italic_v | start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT to allow eq.6 to be written as an equation in one unknown, v(𝐱)𝑣𝐱v(\mathbf{x})italic_v ( bold_x ), since Eu𝒮=Ev𝐸subscript𝑢𝒮𝐸𝑣Eu_{\mathcal{S}}=Evitalic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = italic_E italic_v (but importantly vEv𝑣𝐸𝑣v\neq Evitalic_v ≠ italic_E italic_v except on 𝒮𝒮\mathcal{S}caligraphic_S). We instead do not restrict v|𝒮evaluated-at𝑣𝒮v|_{\mathcal{S}}italic_v | start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT to be u𝒮subscript𝑢𝒮u_{\mathcal{S}}italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT and interpret g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) as a modification to the source term f(𝐱)𝑓𝐱f(\mathbf{x})italic_f ( bold_x ), then solve for the unknown solution u(𝐱)=Eu𝒮𝑢𝐱𝐸subscript𝑢𝒮u(\mathbf{x})=Eu_{\mathcal{S}}italic_u ( bold_x ) = italic_E italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT in eq.6. As proven by von Glehn etal. (2013, Section 3.2), u(𝐱)=Eu(𝐱)𝑢𝐱𝐸𝑢𝐱u(\mathbf{x})=Eu(\mathbf{x})italic_u ( bold_x ) = italic_E italic_u ( bold_x ) since u(𝐱)𝑢𝐱u(\mathbf{x})italic_u ( bold_x ) is the extension of a surface function. The property that u(𝐱)=Eu(x)=u(cp𝒮(𝐱))𝑢𝐱𝐸𝑢𝑥𝑢subscriptcp𝒮𝐱u(\mathbf{x})=Eu(x)=u(\mathrm{cp}_{\mathcal{S}}(\mathbf{x}))italic_u ( bold_x ) = italic_E italic_u ( italic_x ) = italic_u ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ) allows our projection step during the walk in PWoS, detailed in section3.

We show through our numerical examples that taking g(𝐱)=0𝑔𝐱0g(\mathbf{x})=0italic_g ( bold_x ) = 0 for all 𝐱𝒩(𝒮)𝐱𝒩𝒮\mathbf{x}\in\mathcal{N}(\mathcal{S})bold_x ∈ caligraphic_N ( caligraphic_S ),provides qualitatively correct results for graphics applications and quantitatively convergent results in most examples in section5. However, the choice of g(𝐱)=0𝑔𝐱0g(\mathbf{x})=0italic_g ( bold_x ) = 0 causes eq.6 to be ill-posed as discussed above, and we observe some bias in some convergence studies in section5.1 when f𝑓fitalic_f is complex. Interesting future work would involve constructing a more accurate g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) function to improve convergence.

In the traditional CPM, one solves the embedding PDE inside a narrow tubular subset of 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) that is within a constant distance to 𝒮𝒮\mathcal{S}caligraphic_S.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 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ).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 𝒞𝒮𝒞𝒮\mathcal{C}\subset\mathcal{S}caligraphic_C ⊂ caligraphic_S, on which the solution u𝒮subscript𝑢𝒮u_{\mathcal{S}}italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT 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 cp𝒞(𝐱)subscriptcp𝒞𝐱\mathrm{cp}_{\mathcal{C}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ) to the (pre-extension) boundary 𝒞𝒞\mathcal{C}caligraphic_C, 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 𝒮𝒮\mathcal{S}caligraphic_S is embedded in 3superscript3{\mathbb{R}^{3}}blackboard_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT and dim(𝒮)=1,2dimension𝒮12\dim(\mathcal{S})=1,2roman_dim ( caligraphic_S ) = 1 , 2. Recall that we use the word surface to refer both to “surfaces” with dim(𝒮)=1dimension𝒮1\dim(\mathcal{S})=1roman_dim ( caligraphic_S ) = 1 (curves) as well as dim(𝒮)=2dimension𝒮2\dim(\mathcal{S})=2roman_dim ( caligraphic_S ) = 2 surfaces. We also allow mixed codimension where parts of the surface have dim(𝒮)=1dimension𝒮1\dim(\mathcal{S})=1roman_dim ( caligraphic_S ) = 1 and the rest of the surface has dim(𝒮)=2dimension𝒮2\dim(\mathcal{S})=2roman_dim ( caligraphic_S ) = 2. We assume that we can query the closest point function cp𝒮(𝐱)subscriptcp𝒮𝐱\mathrm{cp}_{\mathcal{S}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) for 𝐱3𝐱superscript3\mathbf{x}\in{\mathbb{R}^{3}}bold_x ∈ blackboard_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT.Additionally, for surfaces with dim(𝒮)=2dimension𝒮2\dim(\mathcal{S})=2roman_dim ( caligraphic_S ) = 2, we assume that we can query the unoriented normal direction of the surface 𝐧(𝐱)𝐧𝐱\mathbf{n}(\mathbf{x})bold_n ( bold_x ) for 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S. For surfaces with dim(𝒮)=1dimension𝒮1\dim(\mathcal{S})=1roman_dim ( caligraphic_S ) = 1, we assume that we can query the tangential direction of the surface 𝐭(𝐱)𝐭𝐱\mathbf{t}(\mathbf{x})bold_t ( bold_x ) for 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S.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 𝒮𝒮\mathcal{S}caligraphic_S 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 𝒞𝒞\mathcal{C}caligraphic_C have a lower dimension than the dimension of 𝒮𝒮\mathcal{S}caligraphic_S and support the closest point query cp𝒞subscriptcp𝒞\mathrm{cp}_{\mathcal{C}}roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT. When solving a two-sided boundary value problem for boundaries with dim(𝒞)=1dimension𝒞1\dim(\mathcal{C})=1roman_dim ( caligraphic_C ) = 1, we also assume that we can query the tangent direction of 𝒞𝒞\mathcal{C}caligraphic_C.

Input :surface 𝒮𝒮\mathcal{S}caligraphic_S, boundary 𝒞𝒞\mathcal{C}caligraphic_C, evaluation point 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S, sample walk count NPsubscript𝑁𝑃N_{P}italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT, volume sample count NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT, tolerance ϵitalic-ϵ\epsilonitalic_ϵ

FunctionEstimateSolution(𝒮𝒮\mathcal{S}caligraphic_S, 𝒞𝒞\mathcal{C}caligraphic_C, NPsubscript𝑁𝑃N_{P}italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT, NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT, 𝐱𝐱\mathbf{x}bold_x, ϵitalic-ϵ\epsilonitalic_ϵ):

absent\mathcal{M}\leftarrowcaligraphic_M ← medialAxisPointCloud(𝒮𝒮\mathcal{S}caligraphic_S)

// section3.1

u^sum0subscript^𝑢sum0\widehat{u}_{\mathrm{sum}}\leftarrow 0over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sum end_POSTSUBSCRIPT ← 0

forn1𝑛1n\leftarrow 1italic_n ← 1 to NPsubscript𝑁𝑃N_{P}italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPTdo

u^^𝑢\widehat{u}over^ start_ARG italic_u end_ARG \leftarrow RecursiveEstimate(𝒮𝒮\mathcal{S}caligraphic_S, \mathcal{M}caligraphic_M, 𝒞𝒞\mathcal{C}caligraphic_C, NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT, 𝐱𝐱\mathbf{x}bold_x, ϵitalic-ϵ\epsilonitalic_ϵ)

u^sumsubscript^𝑢sum\widehat{u}_{\mathrm{sum}}over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sum end_POSTSUBSCRIPT \leftarrowu^sum+u^subscript^𝑢sum^𝑢\widehat{u}_{\mathrm{sum}}+\widehat{u}over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sum end_POSTSUBSCRIPT + over^ start_ARG italic_u end_ARG

end for

return u^sum/NPsubscript^𝑢sumsubscript𝑁𝑃\widehat{u}_{\mathrm{sum}}/N_{P}over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sum end_POSTSUBSCRIPT / italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT

FunctionRecursiveEstimate(𝒮𝒮\mathcal{S}caligraphic_S, \mathcal{M}caligraphic_M, 𝒞𝒞\mathcal{C}caligraphic_C, NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT, 𝐱𝐱\mathbf{x}bold_x, ϵitalic-ϵ\epsilonitalic_ϵ):

δ𝛿absent\delta\leftarrowitalic_δ ← distanceToBoundary(𝒮𝒮\mathcal{S}caligraphic_S, \mathcal{M}caligraphic_M, 𝒞𝒞\mathcal{C}caligraphic_C, 𝐱𝐱\mathbf{x}bold_x)

// section3.2

ifδ<ϵ𝛿italic-ϵ\delta<\epsilonitalic_δ < italic_ϵthen

return u(cp𝒞(𝐱))𝑢subscriptcp𝒞𝐱u(\mathrm{cp}_{\mathcal{C}}(\mathbf{x}))italic_u ( roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ) )

l𝑙absentl\leftarrowitalic_l ← localFeatureSize(𝒮𝒮\mathcal{S}caligraphic_S, \mathcal{M}caligraphic_M, 𝐱𝐱\mathbf{x}bold_x)

// section3.1

rmin(l,δ)𝑟𝑙𝛿r\leftarrow\min(l,\delta)italic_r ← roman_min ( italic_l , italic_δ )

𝐲𝐲\mathbf{y}bold_y \leftarrow uniformSphereSample(center=𝐱𝐱\mathbf{x}bold_x , radius=r𝑟ritalic_r)

u^spheresubscript^𝑢sphereabsent\widehat{u}_{\mathrm{sphere}}\leftarrowover^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sphere end_POSTSUBSCRIPT ← RecursiveEstimate(𝒮𝒮\mathcal{S}caligraphic_S, medmed\mathrm{med}roman_med, 𝒞𝒞\mathcal{C}caligraphic_C, NVsubscript𝑁𝑉N_{V}italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT, cp𝒮(𝐲)subscriptcp𝒮𝐲\mathrm{cp}_{\mathcal{S}}(\mathbf{y})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ), ϵitalic-ϵ\epsilonitalic_ϵ )

{𝐳1,𝐳NV}subscript𝐳1subscript𝐳subscript𝑁𝑉absent\{\mathbf{z}_{1},...\mathbf{z}_{N_{V}}\}\leftarrow{ bold_z start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … bold_z start_POSTSUBSCRIPT italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_POSTSUBSCRIPT } ← ballSample(center=𝐱𝐱\mathbf{x}bold_x , radius=r𝑟ritalic_r)

u^ball1NVi=1NVG(𝐱,𝐳i)f(cp𝒮(𝐳i))p(𝐳i)subscript^𝑢ball1subscript𝑁𝑉superscriptsubscript𝑖1subscript𝑁𝑉𝐺𝐱subscript𝐳𝑖𝑓subscriptcp𝒮subscript𝐳𝑖𝑝subscript𝐳𝑖\widehat{u}_{\mathrm{ball}}\leftarrow\frac{1}{N_{V}}\sum_{i=1}^{N_{V}}\frac{G(%\mathbf{x},\mathbf{z}_{i})f({\color[rgb]{0,0.443,0.737}\mathrm{cp}_{\mathcal{S%}}(\mathbf{z}_{i})})}{p(\mathbf{z}_{i})}over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_ball end_POSTSUBSCRIPT ← divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG italic_G ( bold_x , bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_f ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ) end_ARG start_ARG italic_p ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_ARG

return u^sphere+u^ballsubscript^𝑢spheresubscript^𝑢ball\widehat{u}_{\mathrm{sphere}}+\widehat{u}_{\mathrm{ball}}over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_sphere end_POSTSUBSCRIPT + over^ start_ARG italic_u end_ARG start_POSTSUBSCRIPT roman_ball end_POSTSUBSCRIPT

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (4)

Overview

The core idea of our method is to apply the WoS recursive relationship within 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) while utilizing the closest point extension constraint that u(𝐱)=u(cp𝒮(𝐱))𝑢𝐱𝑢subscriptcp𝒮𝐱u(\mathbf{x})=u(\mathrm{cp}_{\mathcal{S}}(\mathbf{x}))italic_u ( bold_x ) = italic_u ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ). To do so, we modify the walk process to use spheres contained within 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) and to project the walk position at each recursion step, as illustrated in fig.2.

The problem we solve is the embedding PDE Δu(𝐱)=f(𝐱)+g(𝐱)Δ𝑢𝐱𝑓𝐱𝑔𝐱\Delta u(\mathbf{x})=f(\mathbf{x})+g(\mathbf{x})roman_Δ italic_u ( bold_x ) = italic_f ( bold_x ) + italic_g ( bold_x ) within 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ). The Monte Carlo estimate of eq.2 holds by assuming g(𝐱)=0𝑔𝐱0g(\mathbf{x})=0italic_g ( bold_x ) = 0 because the embedding PDE is defined with the Cartesian differential operator.To estimate the surface PDE’s solution at point 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S, we consider a 3D ball centered at 𝐱𝐱\mathbf{x}bold_x and fully contained inside 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ). Theoretically, it should be the largest ball fully contained inside 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) that does not cross the extended Dirichlet boundaries 𝒞𝒞\mathcal{C}caligraphic_C, 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 𝐱𝐱\mathbf{x}bold_x (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, u^(𝐲)^𝑢𝐲\widehat{u}(\mathbf{y})over^ start_ARG italic_u end_ARG ( bold_y ), needs to be evaluated at point 𝐲𝐲\mathbf{y}bold_y, which does not lie on 𝒮𝒮\mathcal{S}caligraphic_S in general.The closest point extension constraint of u𝑢uitalic_u provides a convenient relationship here:the embedding PDE’s solution at point 𝐲𝐲\mathbf{y}bold_y should coincide with the surface PDE’s solution at the projected point, cp𝒮(𝐲)subscriptcp𝒮𝐲\mathrm{cp}_{\mathcal{S}}(\mathbf{y})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ).We can therefore project the point 𝐲𝐲\mathbf{y}bold_y onto 𝒮𝒮\mathcal{S}caligraphic_S at the end of each recursion step hence u^(𝐲)=u^(cp𝒮(𝐲))^𝑢𝐲^𝑢subscriptcp𝒮𝐲\widehat{u}(\mathbf{y})=\widehat{u}(\mathrm{cp}_{\mathcal{S}}(\mathbf{y}))over^ start_ARG italic_u end_ARG ( bold_y ) = over^ start_ARG italic_u end_ARG ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) ).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)u^(𝐱)=u^(cp𝒮(𝐲))+1NVi=1NVG(𝐱,𝐳i)f(cp𝒮(𝐳i))p(𝐳i).^𝑢𝐱^𝑢subscriptcp𝒮𝐲1subscript𝑁𝑉superscriptsubscript𝑖1subscript𝑁𝑉𝐺𝐱subscript𝐳𝑖𝑓subscriptcp𝒮subscript𝐳𝑖𝑝subscript𝐳𝑖\widehat{u}(\mathbf{x})=\widehat{u}(\mathrm{cp}_{\mathcal{S}}(\mathbf{y}))+%\frac{1}{N_{V}}\sum_{i=1}^{N_{V}}\frac{G(\mathbf{x},\mathbf{z}_{i})f(\mathrm{%cp}_{\mathcal{S}}(\mathbf{z}_{i}))}{p(\mathbf{z}_{i})}.over^ start_ARG italic_u end_ARG ( bold_x ) = over^ start_ARG italic_u end_ARG ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) ) + divide start_ARG 1 end_ARG start_ARG italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG italic_G ( bold_x , bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_f ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ) end_ARG start_ARG italic_p ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_ARG .

We choose p(𝐳i)1/𝐱𝐳i2proportional-to𝑝subscript𝐳𝑖1subscriptdelimited-∥∥𝐱subscript𝐳𝑖2p(\mathbf{z}_{i})\propto 1/\lVert\mathbf{x}-\mathbf{z}_{i}\rVert_{2}italic_p ( bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ∝ 1 / ∥ bold_x - bold_z start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in our implementation.

Analogous to the original WoS method, we terminate the recursion when the point 𝐱𝐱\mathbf{x}bold_x falls within a distance ϵitalic-ϵ\epsilonitalic_ϵ of the (extended) Dirichlet boundary by taking the boundary value u(cp𝒮(𝐱))𝑢subscriptcp𝒮𝐱u(\mathrm{cp}_{\mathcal{S}}(\mathbf{x}))italic_u ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ). 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 𝒮𝒮\mathcal{S}caligraphic_S is embedded in 2superscript2{\mathbb{R}^{2}}blackboard_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT in fig.2. Notably, PWoS is a generalization of the WoS algorithm: when dim(𝒮)=ddimension𝒮𝑑\dim(\mathcal{S})=droman_dim ( caligraphic_S ) = italic_d (i.e., the codimension-zero case), the local feature size is infinite, the distance to the Dirichlet boundary 𝒞𝒞\mathcal{C}caligraphic_C can be computed with the closest point query cp𝒞subscriptcp𝒞\mathrm{cp}_{\mathcal{C}}roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT, and the last closest point projection of 𝐲𝐲\mathbf{y}bold_y has no effect since cp𝒮(𝐲)=𝐲subscriptcp𝒮𝐲𝐲\mathrm{cp}_{\mathcal{S}}(\mathbf{y})=\mathbf{y}roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) = bold_y.When dim(𝒮)<ddimension𝒮𝑑\dim(\mathcal{S})<droman_dim ( caligraphic_S ) < italic_d, 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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (5)

Local Feature SizeAvg. Step Count
0.9931.1398
0.547.84
0.25128.981
0.125462.781
0.06251818.73

3.1. Local Feature Size Estimation

To determine the radius of the sphere centered at 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S that is fully contained inside 𝒩(𝒮)𝒩𝒮\mathcal{N}(\mathcal{S})caligraphic_N ( caligraphic_S ) at each recursion step of the walk, we need a conservative lower bound estimate for the local feature size at 𝐱𝐱\mathbf{x}bold_x.One naive approach would be to use a small enough positive constant value for all 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S, 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 1111 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 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S 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 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT inside a ball in 3superscript3{\mathbb{R}^{3}}blackboard_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT having a radius equal to half the length of the diagonal of the bounding box of 𝒮𝒮\mathcal{S}caligraphic_S, so the entire surface is fully enclosed.For each point 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, we use the closest point query cp𝒮(𝐱i)subscriptcp𝒮subscript𝐱𝑖\mathrm{cp}_{\mathcal{S}}(\mathbf{x}_{i})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) to compute its two opposing normal directions at cp𝒮(𝐱i)subscriptcp𝒮subscript𝐱𝑖\mathrm{cp}_{\mathcal{S}}(\mathbf{x}_{i})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ).Specifically, we normalize the vector 𝐱icp𝒮(𝐱i)subscript𝐱𝑖subscriptcp𝒮subscript𝐱𝑖\mathbf{x}_{i}-\mathrm{cp}_{\mathcal{S}}(\mathbf{x}_{i})bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) 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 cp𝒮(𝐱i)subscriptcp𝒮subscript𝐱𝑖\mathrm{cp}_{\mathcal{S}}(\mathbf{x}_{i})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), 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 𝒮𝒮\mathcal{S}caligraphic_S. Then, we iteratively shrink the size of the sphere, moving its center to maintain tangency at the surface point cp𝒮(𝐱i)subscriptcp𝒮subscript𝐱𝑖\mathrm{cp}_{\mathcal{S}}(\mathbf{x}_{i})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ), 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 𝒮𝒮\mathcal{S}caligraphic_S 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 𝒮𝒮\mathcal{S}caligraphic_S. 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 s>1𝑠1s>1italic_s > 1. Then, for any pair of medial balls Br1(𝐱1)subscript𝐵subscript𝑟1subscript𝐱1B_{r_{1}}(\mathbf{x}_{1})italic_B start_POSTSUBSCRIPT italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) and Br2(𝐱2)subscript𝐵subscript𝑟2subscript𝐱2B_{r_{2}}(\mathbf{x}_{2})italic_B start_POSTSUBSCRIPT italic_r start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ), we remove the smaller ball from the set of medial balls if it is fully contained in the other. That is, if sr1<sr2+𝐱1𝐱22𝑠subscript𝑟1𝑠subscript𝑟2subscriptdelimited-∥∥subscript𝐱1subscript𝐱22s\cdot r_{1}<s\cdot r_{2}+\lVert\mathbf{x}_{1}-\mathbf{x}_{2}\rVert_{2}italic_s ⋅ italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT < italic_s ⋅ italic_r start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + ∥ bold_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - bold_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, we remove the ball Br1(𝐱1)subscript𝐵subscript𝑟1subscript𝐱1B_{r_{1}}(\mathbf{x}_{1})italic_B start_POSTSUBSCRIPT italic_r start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ).

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 𝒮𝒮\mathcal{S}caligraphic_S).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 𝒮𝒮\mathcal{S}caligraphic_S.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 s𝑠sitalic_s allows us to control the pruning strength. Unless otherwise stated, we use the value s=1.15𝑠1.15s=1.15italic_s = 1.15 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 (0.90.90.90.9 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 λ𝜆\lambdaitalic_λ, we return λ𝜆\lambdaitalic_λ as the local feature size estimate instead. This process essentially rounds sharp corners with rounding radius λ𝜆\lambdaitalic_λ. 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 𝒞𝒞\mathcal{C}caligraphic_C are extended in the normal direction of 𝒮𝒮\mathcal{S}caligraphic_S and the solution in the embedding space on this extended boundary is determined by the closest point extension of Dirichlet values on 𝒞𝒞\mathcal{C}caligraphic_C.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 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S, we first find the closest point that lies on the boundary before the extension, cp𝒞(𝐱)subscriptcp𝒞𝐱\mathrm{cp}_{\mathcal{C}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ).The subset of the extended boundary that is extended from cp𝒞(𝐱)subscriptcp𝒞𝐱\mathrm{cp}_{\mathcal{C}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ) takes the shape of a line segment when dim(𝒮)=2dimension𝒮2\dim(\mathcal{S})=2roman_dim ( caligraphic_S ) = 2 and a disk when dim(𝒮)=1dimension𝒮1\dim(\mathcal{S})=1roman_dim ( caligraphic_S ) = 1.We set the line segment’s half-length or the disk radius to the local feature size at cp𝒞(𝐱)subscriptcp𝒞𝐱\mathrm{cp}_{\mathcal{C}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ) using the algorithm in section3.1.We can compute the distance from the point 𝐱𝒮𝐱𝒮\mathbf{x}\in\mathcal{S}bold_x ∈ caligraphic_S to this line segment or disk without explicitly constructing the extended boundary geometry.When dim(𝒮)=2dimension𝒮2\dim(\mathcal{S})=2roman_dim ( caligraphic_S ) = 2, the distance δ𝛿\deltaitalic_δ to the extended boundary is given by

(8)𝐫=cp𝒞(𝐱)𝐱,δ=𝐫clamp(𝐫𝐧,l,l)𝐧2,formulae-sequence𝐫subscriptcp𝒞𝐱𝐱𝛿subscriptdelimited-∥∥𝐫clamp𝐫𝐧𝑙𝑙𝐧2\begin{split}\mathbf{r}&=\mathrm{cp}_{\mathcal{C}}(\mathbf{x})-\mathbf{x},\\\delta&=\lVert\mathbf{r}-\mathrm{clamp}(\mathbf{r}\cdot\mathbf{n},-l,l)\cdot%\mathbf{n}\rVert_{2},\end{split}start_ROW start_CELL bold_r end_CELL start_CELL = roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ) - bold_x , end_CELL end_ROW start_ROW start_CELL italic_δ end_CELL start_CELL = ∥ bold_r - roman_clamp ( bold_r ⋅ bold_n , - italic_l , italic_l ) ⋅ bold_n ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , end_CELL end_ROW

where 𝐧𝐧\mathbf{n}bold_n and l𝑙litalic_l are the surface normal and the local feature size estimate at cp𝒞(𝐱)subscriptcp𝒞𝐱\mathrm{cp}_{\mathcal{C}}(\mathbf{x})roman_cp start_POSTSUBSCRIPT caligraphic_C end_POSTSUBSCRIPT ( bold_x ), respectively.When dim(𝒮)=1dimension𝒮1\dim(\mathcal{S})=1roman_dim ( caligraphic_S ) = 1, the normal direction is not uniquely defined, so we instead use a similar method based on the surface tangent 𝐭𝐭\mathbf{t}bold_t:

(9)𝐪=𝐫(𝐫𝐭)𝐭,δ={|𝐫𝐭|,if𝐪2<l,𝐫l(𝐪/𝐪2)2,otherwise.formulae-sequence𝐪𝐫𝐫𝐭𝐭𝛿cases𝐫𝐭ifsubscriptdelimited-∥∥𝐪2𝑙subscriptdelimited-∥∥𝐫𝑙𝐪subscriptdelimited-∥∥𝐪22otherwise\begin{split}\mathbf{q}&=\mathbf{r}-(\mathbf{r}\cdot\mathbf{t})\mathbf{t},\\\delta&=\begin{cases}\lvert\mathbf{r}\cdot\mathbf{t}\rvert,&\mathrm{if}\;%\lVert\mathbf{q}\rVert_{2}<l,\\\lVert\mathbf{r}-l\cdot(\mathbf{q}/\lVert\mathbf{q}\rVert_{2})\rVert_{2},&%\mathrm{otherwise}.\end{cases}\end{split}start_ROW start_CELL bold_q end_CELL start_CELL = bold_r - ( bold_r ⋅ bold_t ) bold_t , end_CELL end_ROW start_ROW start_CELL italic_δ end_CELL start_CELL = { start_ROW start_CELL | bold_r ⋅ bold_t | , end_CELL start_CELL roman_if ∥ bold_q ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT < italic_l , end_CELL end_ROW start_ROW start_CELL ∥ bold_r - italic_l ⋅ ( bold_q / ∥ bold_q ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , end_CELL start_CELL roman_otherwise . end_CELL end_ROW end_CELL end_ROW

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 Δuσu=fΔ𝑢𝜎𝑢𝑓\Delta u-\sigma u=froman_Δ italic_u - italic_σ italic_u = italic_f, where σ𝜎\sigmaitalic_σ 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)u(𝐱)=cr,σ|Br(𝐱)|Br(𝐱)u(𝐲)d𝐲+Br(𝐱)f(𝐳)Gσ(𝐱,𝐳)d𝐳,𝑢𝐱subscript𝑐𝑟𝜎subscript𝐵𝑟𝐱subscriptsubscript𝐵𝑟𝐱𝑢𝐲differential-d𝐲subscriptsubscript𝐵𝑟𝐱𝑓𝐳subscript𝐺𝜎𝐱𝐳differential-d𝐳u(\mathbf{x})=\frac{c_{r,\sigma}}{\lvert\partial B_{r}(\mathbf{x})\rvert}\int_%{\partial B_{r}(\mathbf{x})}u(\mathbf{y})\,\mathrm{d}\mathbf{y}+\int_{B_{r}(%\mathbf{x})}f(\mathbf{z})G_{\sigma}(\mathbf{x},\mathbf{z})\,\mathrm{d}\mathbf{%z},italic_u ( bold_x ) = divide start_ARG italic_c start_POSTSUBSCRIPT italic_r , italic_σ end_POSTSUBSCRIPT end_ARG start_ARG | ∂ italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) | end_ARG ∫ start_POSTSUBSCRIPT ∂ italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) end_POSTSUBSCRIPT italic_u ( bold_y ) roman_d bold_y + ∫ start_POSTSUBSCRIPT italic_B start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT ( bold_x ) end_POSTSUBSCRIPT italic_f ( bold_z ) italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT ( bold_x , bold_z ) roman_d bold_z ,

where Gσsubscript𝐺𝜎G_{\sigma}italic_G start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT is the Yukawa potential and cr,σsubscript𝑐𝑟𝜎c_{r,\sigma}italic_c start_POSTSUBSCRIPT italic_r , italic_σ end_POSTSUBSCRIPT is a positive number smaller than 1111. To evaluate the first term, instead of multiplying the contribution from the recursion by cr,σsubscript𝑐𝑟𝜎c_{r,\sigma}italic_c start_POSTSUBSCRIPT italic_r , italic_σ end_POSTSUBSCRIPT as suggested by Sawhney and Crane (2020), we use a Russian Roulette strategy: we terminate the path early with probability 1cr,σ1subscript𝑐𝑟𝜎1-c_{r,\sigma}1 - italic_c start_POSTSUBSCRIPT italic_r , italic_σ end_POSTSUBSCRIPT with zero contribution and otherwise we use the estimate of the solution without multiplying by cr,σsubscript𝑐𝑟𝜎c_{r,\sigma}italic_c start_POSTSUBSCRIPT italic_r , italic_σ end_POSTSUBSCRIPT. 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 σ𝜎\sigmaitalic_σ 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, f𝒮=𝒮𝐡𝒮subscript𝑓𝒮subscript𝒮subscript𝐡𝒮f_{\mathcal{S}}=\nabla_{\mathcal{S}}\cdot\mathbf{h}_{\mathcal{S}}italic_f start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = ∇ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ⋅ bold_h start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT, where 𝐡𝒮subscript𝐡𝒮\mathbf{h}_{\mathcal{S}}bold_h start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT 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 u(𝐱)=u(cp𝒮(𝐱))𝑢𝐱𝑢subscriptcp𝒮𝐱u(\mathbf{x})=u(\mathrm{cp}_{\mathcal{S}}(\mathbf{x}))italic_u ( bold_x ) = italic_u ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_x ) ). 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 u^(cp𝒮(𝐲))^𝑢subscriptcp𝒮𝐲\widehat{u}(\mathrm{cp}_{\mathcal{S}}(\mathbf{y}))over^ start_ARG italic_u end_ARG ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) ) 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 u^(cp𝒮(𝐲))^𝑢subscriptcp𝒮𝐲\widehat{u}(\mathrm{cp}_{\mathcal{S}}(\mathbf{y}))over^ start_ARG italic_u end_ARG ( roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ) ) available at cp𝒮(𝐲)subscriptcp𝒮𝐲\mathrm{cp}_{\mathcal{S}}(\mathbf{y})roman_cp start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ( bold_y ). 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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (6)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (7)

(a) Helix
Ends
mylaplace

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (8)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (9)

(b) Helix
Ends
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (10)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (11)

(c) Z-Order (s=1.15𝑠1.15s=1.15italic_s = 1.15)
Ends
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (12)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (13)

(d) Z-Order (s=1.05𝑠1.05s=1.05italic_s = 1.05)
Ends
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (14)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (15)

(e) Circle
Two-sided
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (16)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (17)

(f) Torus
Torus Knot
mylaplace

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (18)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (19)

(g) Dziuk Surface
Closed
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (20)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (21)

(h) Dziuk Surface
No Boundary
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (22)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (23)

(i) Sphere
Closed
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (24)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (25)

(j) Sphere
Open
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (26)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (27)

(k) Sphere
Closed
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (28)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (29)

(l) Sphere
No Boundary
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (30)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (31)

(m) Punched
Closed
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (32)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (33)

(n) Punched
Open
Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (34)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (35)

(o) Punched
Closed
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (36)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (37)

(p) Punched
No Boundary
Screened Poisson

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (38)Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (39)

High

Mid

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (40)

Low

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (41)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (42)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (43)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (44)

-No mean   value filtering
-1x mean   value filtering
-10x mean   value filtering

Initial Samples

Filter Samples

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 NPsubscript𝑁𝑃N_{P}italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT 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 ϵ=0.001italic-ϵ0.001\epsilon=0.001italic_ϵ = 0.001 and NV=32subscript𝑁𝑉32N_{V}=32italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT = 32 except for Laplace equations, for which we use NV=0subscript𝑁𝑉0N_{V}=0italic_N start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT = 0. While the method shows the expected convergence rate of 𝒪(1/NP)𝒪1subscript𝑁𝑃\mathcal{O}(1/\sqrt{N_{P}})caligraphic_O ( 1 / square-root start_ARG italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT end_ARG ) 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 g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) 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 𝒪(1/NP)𝒪1subscript𝑁𝑃\mathcal{O}(1/\sqrt{N_{P}})caligraphic_O ( 1 / square-root start_ARG italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT end_ARG ), where NPsubscript𝑁𝑃N_{P}italic_N start_POSTSUBSCRIPT italic_P end_POSTSUBSCRIPT 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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (45)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (46)

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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (47)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (48)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (49)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (50)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (51)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (52)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (53)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (54)

Triangle Mesh
(Ours)

Point Cloud
(Ours)

Grid-Based CPM
(King etal., 2023)

Exact
(Sharp etal., 2019)

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 𝒞𝒞\mathcal{C}caligraphic_C.The steps are summarized as

  1. (1)

    Δ𝒮u𝒮(1/t)u𝒮=0subscriptΔ𝒮subscript𝑢𝒮1𝑡subscript𝑢𝒮0\Delta_{\mathcal{S}}u_{\mathcal{S}}-(1/t)u_{\mathcal{S}}=0roman_Δ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT - ( 1 / italic_t ) italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = 0 where u𝒮=1on𝒞subscript𝑢𝒮1on𝒞u_{\mathcal{S}}=1\;\mathrm{on}\;\mathcal{C}italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = 1 roman_on caligraphic_C,

  2. (2)

    𝐗=(𝒮u𝒮)/𝒮u𝒮2𝐗subscript𝒮subscript𝑢𝒮subscriptdelimited-∥∥subscript𝒮subscript𝑢𝒮2\mathbf{X}=-(\nabla_{\mathcal{S}}u_{\mathcal{S}})/\lVert\nabla_{\mathcal{S}}u_%{\mathcal{S}}\rVert_{2}bold_X = - ( ∇ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ) / ∥ ∇ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, and

  3. (3)

    Δ𝒮ϕ𝒮=𝒮𝐗subscriptΔ𝒮subscriptitalic-ϕ𝒮subscript𝒮𝐗\Delta_{\mathcal{S}}\phi_{\mathcal{S}}=\nabla_{\mathcal{S}}\cdot\mathbf{X}roman_Δ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT italic_ϕ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = ∇ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT ⋅ bold_X where ϕ𝒮=0on𝒞subscriptitalic-ϕ𝒮0on𝒞\phi_{\mathcal{S}}=0\;\mathrm{on}\;\mathcal{C}italic_ϕ start_POSTSUBSCRIPT caligraphic_S end_POSTSUBSCRIPT = 0 roman_on caligraphic_C,

where t𝑡titalic_t is a small positive constant and ϕitalic-ϕ\phiitalic_ϕ 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 u𝑢uitalic_u during Step 1 using the method in section3.3.2, without needing u𝑢uitalic_u itself. We evaluate the gradient at mesh vertices and normalize it to get 𝐗𝐗\mathbf{X}bold_X 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 𝐗𝐗\mathbf{X}bold_X at a point, we interpolate 𝐗𝐗\mathbf{X}bold_X 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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (55)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (56)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (57)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (58)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (59)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (60)

Frame 7

Frame 19

Frame 31

Frame 43

Frame 55

Frame 200

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.

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (61)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (62)

Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (63)

(e)

(g)

(h)

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., g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) 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 g(𝐱)𝑔𝐱g(\mathbf{x})italic_g ( bold_x ) tends to zero continuously as 𝐱𝐱\mathbf{x}bold_x 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 indsuperscript𝑑\mathbb{R}^{d}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT. 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]
Projected Walk on Spheres: A Monte Carlo Closest Point Method for Surface PDEs (2024)

References

Top Articles
Latest Posts
Recommended Articles
Article information

Author: Nathanael Baumbach

Last Updated:

Views: 6589

Rating: 4.4 / 5 (55 voted)

Reviews: 86% of readers found this page helpful

Author information

Name: Nathanael Baumbach

Birthday: 1998-12-02

Address: Apt. 829 751 Glover View, West Orlando, IN 22436

Phone: +901025288581

Job: Internal IT Coordinator

Hobby: Gunsmithing, Motor sports, Flying, Skiing, Hooping, Lego building, Ice skating

Introduction: My name is Nathanael Baumbach, I am a fantastic, nice, victorious, brave, healthy, cute, glorious person who loves writing and wants to share my knowledge and understanding with you.