Anugrah Ryan C Dunn
email: rcdunn@ucsd.edu

Jo Joshy

Jui-Te Lin

Cédric Girerd

Tania K Morimoto

John T Hwang

Jo Anugrah

Joshy

Scalable Enforcement of Geometric Non-interference Constraints for Gradient-Based Optimization

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Accurate detection of physical interference between two or more bodies is crucial in the design of many engineering systems. Modeling interference between physical bodies is, therefore, an important problem in computational design.

Non-interference constraints appear in numerical optimization problems that manipulate an object within an environment containing other objects such that there is no collision. Efficient and accurate modeling of the non-interference constraints is critical for fast and reliable solutions in the overarching optimization problem. Prior literature on these problems describe these constraints using inconsistent terminology, e.g., anatomical constraints [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF][START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF], spatial integration constraints [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF], boundary constraints [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF][START_REF] Andrew | Coupled wind turbine design and layout optimization with nonhomogeneous wind turbines[END_REF], and interference checks [START_REF] Fadel | Packing Optimization of Free-Form Objects in Engineering Design[END_REF]. We observe that these terms represent the same underlying constraint. We propose to call these constraints geometric non-interference constraints since they are employed in design optimization to ensure a design where there exists no interference between two or more geometric shapes or paths of motion.

In our study, a geometric shape is associated with the design configuration of an engineering system at a particular instance of time. The geometric shapes of interest in this paper are curves in two dimensions, or orientable surfaces in three dimensions. We assume that the geometric shapes are non-self-intersecting but make no assumptions on whether they are open or closed. A path of motion or trajectory is the set of points that traces the motion of a point on the engineering system as the system changes configuration over time. The paths considered in this paper are simply curves in two or three dimensions. We use the term layout to refer to a set of geometric shapes.

Based on the definitions above, we identify three major classes of optimization problems with geometric non-interference constraints: layout optimization, shape optimization, and optimal path planning. All three classes are parts of the scope of problems we address in this paper. Layout optimization optimizes the positions of geometric shapes via translation subject to geometric noninterference, with or without additional boundary constraints. For example, the wind farm layout optimization problem (WFLOP) consists in positioning the wind turbines within a wind farm in an optimal way such that interference between turbines and the boundary of the farm is avoided [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF][START_REF] Andrew | Coupled wind turbine design and layout optimization with nonhomogeneous wind turbines[END_REF][START_REF] Cazzaro | Variable neighborhood search for large offshore wind farm layout optimization[END_REF][START_REF] Thomas | A Comparison of Eight Optimization Methods Applied to a Wind Farm Layout Optimization Problem[END_REF][START_REF] Baker | Best Practices for Wake Model and Optimization Algorithm Selection in Wind Farm Layout Optimization[END_REF]. Another example of a layout optimization problem is the packing problem. Packing problems consist of positioning objects within a space to ensure the minimum space is occupied or the maximum number of objects are placed without geometric interference [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF][START_REF] Fadel | Packing Optimization of Free-Form Objects in Engineering Design[END_REF][START_REF] Brandt | The air cargo load planning problem[END_REF].

Shape optimization seeks to optimize geometric shapes subject to geometric non-interference, with or without additional boundary constraints. For example, shape optimization of an aircraft fuselage optimizes the shape of a fuselage with constraints ensuring that the passengers, crew, payload, and all the subsystems fit inside the fuselage [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF].

Optimal path planning optimizes the trajectory of a point or a set of points subject to geometric non-interference, with or without additional boundary constraints. Robot motion planning is a class of problems that falls under optimal path planning and is widely researched [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF]. The design optimization of surgical robots is an example of a problem involving robot motion planning that has had recent attention [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF][START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF]. In the design optimization of surgical robots, non-interference constraints are imposed such that the robot does not collide with the anatomy of a patient during operation. Additionally, it is desirable for the robot to maintain a safe distance from the anatomy, motivating the use of a distance-based non-interference constraint formulation in such problems.

The problems just mentioned are solved using numerical optimization algorithms. Historically, gradient-free algorithms have been more commonly used to solve such problems, e.g., in layout optimization [START_REF] Lodi | Two-dimensional packing problems: A survey[END_REF][START_REF] Cagan | A survey of computational approaches to three-dimensional layout problems[END_REF][START_REF] Fasano | Solving non-standard packing problems by global optimization and heuristics[END_REF] and in robot motion planning [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF]. A major reason behind this was the difficulty in efficiently computing the derivatives for a complex model. As models become more complex, that is, with more disciplines and design variables, solutions become impracticable with gradient-free algorithms since these algorithms scale poorly with the number of design variables. However, the recent emergence of modeling frameworks such as OpenMDAO [START_REF] Gray | OpenMDAO: An Open-Source Framework for Multidisciplinary Design, Analysis, and Optimization[END_REF] has enabled efficient design of large-scale and multidisciplinary systems using gradient-based optimization, including some of the aforementioned problems with geometric non-interference constraints [2-4, 8, 9].

Geometric non-interference constraint functions for gradient-based optimization require special consideration. These functions must be continuously differentiable or smooth in order to be used with a gradient-based optimization algorithm. They should also be efficient to compute because optimization algorithms evaluate constraint functions and their derivatives repeatedly over many optimization iterations. During some iterations, the optimizer may violate an interference constraint, and useful gradient information on such iterations is still required despite it being infeasible. Consequently, any non-interference constraint function must be defined in the event of an overlap between objects and provide necessary gradient information.

Design body

Feasible space Ω 𝜙 > 0

Boundary Γ 𝜙 = 0 𝜙 < 0

Infeasible space

Figure 1: An ideal constraint function ϕ indicates geometric interferences using signed distances of representative points x (i) on a body from the boundary Γ defining the feasible space Ω.

Figure 1 shows a diagram with two iterations of a design body in an optimization problem. One of the designs shown is feasible while the other is not. The feasible design is the one where the design body is completely inside the feasible space whereas the infeasible design has at least one point on the design body lying outside the feasible space. For the ϕ defined in Fig. 1, enforcing the optimization constraint ϕ(x (i) ) ≥ ϵ for certain representative points x (i) chosen on the surface of the design body guarantees non-interference by ensuring that all x (i) stay within the feasible region for the final optimized design. The constant ϵ can be any small positive value appropriate for a given problem.

Existing non-interference constraint formulations suffer from various limitations. The formulation of quasi-phi-functions by Stoyan et al. [START_REF] Stoyan | Optimized Object Packings Using Quasi-Phi-Functions[END_REF] provides an analytical form to represent an interference for simple geometric shapes. Quasi-phi-functions are continuous but only piecewise continuously differentiable. These functions are also not generalized to represent any arbitrary shape. The formulation by Brelje et al. [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF] is generalized to any triangulated 3D geometric shape, but has computational limitations. The computational complexity of their method is O(N Γ ), where N Γ is the number of elements in the triangulation. They are able to overcome this scaling issue by making use of graphics processing units (GPUs) but demonstrate their formulation on a geometric shape with only 626 elements in the triangulation. In their recent work on the WFLOP, Risco et al. [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF] formulate a generic explicit method for geometric shapes in 2D, but the method suffers from the same scaling issues as in [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF] and contains discontinuous derivatives. The formulation by Bergeles et al. [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF] employs a distance potential function that is calculated with the k-nearest neighbors. With the use of a k-d tree structure, the computational complexity of the k-nearest neighbor search scales better than linearly, O(k log(N Γ )) on average, but the structure is not suitable for gradient-based optimization because the derivatives are discontinuous when the set of k-nearest neighbors switches.

Outside the domain of non-interference constraint formulations currently employed in optimization, we discovered a significant body of research conducted on a remarkably similar problem by the computer graphics community. Surface reconstruction in the field of computer graphics is the process of converting a set of points into a surface for graphical representation. A popular approach for surface reconstruction is the representation of surfaces by an implicit function. Implicit surface reconstruction methods such as Poisson [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF], Multi-level Partition of Unity (MPU) [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF], and Smooth Signed Distance (SSD) [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF], to name a few, construct an implicit function from a point cloud to represent a surface. We observed that some of these distance-based formulations can be applied to overcome prior limitations in enforcing geometric non-interference constraints in gradient-based optimization.

The objective of this paper is to devise a general methodology based on an appropriate surface reconstruction method to generate a smooth and fast-to-evaluate geometric non-interference constraint function from an oriented point cloud. It is desired that the function locally approximates the signed distance to a geometric shape and that its evaluation time scales independently of the number of points sampled over the shape N Γ . The function must also be an accurate implicit representation of the surface implied by the given point cloud. The contribution of this paper is a new formulation for representing geometric non-interference constraints in gradient-based optimization. We investigate various properties of the proposed formulation, its efficiency compared to existing non-interference constraint formulations, and its accuracy compared to state-of-the-art surface reconstruction methods. Additionally, we demonstrate the computational speedup of our formulation in an experiment with a path planning optimization and shape optimization problem.

The remainder of this paper proceeds as follows. Section 2 reviews existing geometric noninterference constraint formulations for optimization. In this section, we also provide a thorough survey of implicit surface reconstruction methods in order to identify methodologies to be brought into our formulation. Section 3 presents our methodology to generate an implicit level set function for representing geometric non-interference constraints. Section 4 provides numerical results that quantify the accuracy and efficiency of our formulation. We demonstrate our accuracy on common benchmarking models from the computer graphics community and also on geometries associated with aircraft design. We finally demonstrate the application of our method using a medical robot design optimization problem. In Sec. 5, we summarize our approach, its impact, and discuss future work.

Related Work

This section presents an overview of related work to the defined research problem. We begin by reviewing prior methods for enforcing non-interference constraints in gradient-based optimization in subsection 2.1. We then review the problem of surface reconstruction and its complexities in subsection 2.2. In subsection 2.3, we present a literature review for methods that approximate the signed distance function to contextualize our formulation within the field of surface reconstruction.

Previous methods for enforcing non-interference constraints in gradient-based optimization

We identify two preexisting methods for enforcing geometric non-interference constraints in gradient-based optimization that are both continuous and differentiable. Previous constraint formulations that utilize the nearest neighbor distance, e.g., Risco et al. [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF] and Bergeles et al. [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF], have been used in optimization, but we note the that they are non-differentiable and may incur numerical difficulties in gradient-based optimization.

Brelje et al. [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF] implement a general mesh-based constraint formulation for non-interference constraints between two triangulations of objects. Two nonlinear constraints define their formulation. The first constraint is that the minimum distance of the design shape to the geometric shape is greater than zero, and the second constraint is that the intersection length between the two bodies is zero, i.e., there is no intersection. A binary check, e.g., ray tracing, must be used to reject optimization iterations where the design shape is entirely in the infeasible region, where the previous two constraints are satisfied. As noted by Brelje et al., this formulation may be susceptible to representing very thin objects, where the intersection length is very sensitive to the step size of the optimizer. Additionally, the constraint function has a computational complexity of O(N Γ ), which may be addressed by the use of graphics processing units (GPUs).

Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF] implement a modified signed distance function, making it differentiable throughout. Using an oriented set of points to represent the bounds of the feasible region, the constraint function is a distance-based weighted sum of signed distances between the points and a set of points on the design shape. This representation is inexact and is found to compromise accuracy for a smoothness in the constraint representation in practice. Additionally, their formulation has a computational complexity of O(N Γ ).

Surface reconstruction

Our research goal-to derive a smooth level set function from a set of oriented points-closely aligns with the problem of surface reconstruction in computer graphics. Surface reconstruction is done in many ways, and we refer the reader to [START_REF] Berger | A Survey of Surface Reconstruction from Point Clouds[END_REF][START_REF] Huang | Surface Reconstruction from Point Clouds: A Survey and a Benchmark[END_REF] for a full survey on surface reconstruction methods from point clouds. We, in particular, focus on surface reconstruction with implicit function representations from point clouds. Implicit surface reconstruction is done by constructing an indicator function between the interior and exterior of a surface, whose isocontour represents a smooth surface implied by the point cloud. The methodologies for surface reconstruction use implicit functions as a means to an end; however, the focus of our investigation is on the implicit function itself for enforcing non-interference constraints. We identify that the direct connection between non-interference constraints and implicit functions in surface reconstruction is that the reconstructed surface represents the boundary between the feasible and infeasible region in a continuous and differentiable way.

The surface reconstruction problem begins with a representation of a geometric shape. Geometric non-interference constraints may be represented by geometric shapes using scanned samples of the surface of an anatomy [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF][START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF], outer mold line (OML) meshes [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF], user defined polygons [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF], and a sampled set of points of seabed depths [START_REF] Cazzaro | Variable neighborhood search for large offshore wind farm layout optimization[END_REF]. Many geometric shape representations, including those mentioned, can be sampled and readily converted into an oriented point cloud and posed as a surface reconstruction problem.

The construction of any point cloud comes with additional complexities. For example, machine tolerance of scanners introduce error into scans, and meshing algorithms produce different point cloud representations for the same geometric shape. As a result, implicit surface reconstruction methods often take into consideration nonuniform sampling, noise, outliers, misalignment between scans, and missing data in point clouds. Implicit surface reconstruction methods have been shown to address these issues well, including hole-filling [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF][START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF][START_REF] Davis | Filling holes in complex surfaces using volumetric diffusion[END_REF], reconstructing surfaces from noisy samples [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF][START_REF] Quynh Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF][START_REF] Pan | Phase-Field Guided Surface Reconstruction Based on Implicit Hierarchical B-Splines[END_REF], reconstructing sharp corners and edges [START_REF] Quynh Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF], and reconstructing surfaces without normal vectors in the point cloud [START_REF] Pan | Phase-Field Guided Surface Reconstruction Based on Implicit Hierarchical B-Splines[END_REF][START_REF] Hoppe | Surface Reconstruction from Unorganized Points[END_REF].

Basis functions are commonly used to define the space of implicit functions for implicit surface reconstruction. Basis functions are constructed from a discrete set of points scattered throughout the domain, whose distribution and locations play an important role to defining the implicit function. Examples of these points include control points for B-splines, centers for radial basis functions, and shifts for wavelets. Implicit surface reconstruction methods distribute these points in various ways.

One approach is to adaptively subdivide the implicit function's domain using an octree structure. Octrees, as used by [17-19, 26, 28, 29], recursively subdivide the domain into octants using various heuristics in order to form neighborhoods of control points near the surface. Heuristics include point density [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF], error-controlled [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF], and curvature-based [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF] subdivisions. Octrees are notable because the error of the surface reconstruction decays with the sampling width between control points, which decreases exponentially with respect to the octree depth [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF]. Additionally, the neighborhoods of control points from octrees can be solved for and evaluated in parallel using graphics processing units (GPUs), which allows for on-demand surface reconstruction as demonstrated in [START_REF] Zhou | Data-Parallel Octrees for Surface Reconstruction[END_REF].

Another approach for distributing the points that control the implicit function is to locate them directly on the points in the point cloud. In the formulation by Carr et al. [START_REF] Carr | Reconstruction and Representation of 3D Objects with Radial Basis Functions[END_REF], a chosen subset of points in the point cloud and points projected in the direction of the normal vectors are used to place the radial basis function (RBF) centers, resulting in fewer centers than octrees that are still distributed near the surface. The explicit formulation by Hicken and Kaur [START_REF] Hicken | An Explicit Level-Set Formula to Approximate Geometries[END_REF] uses all points in the point cloud to define the implicit function and shows favorable decay in surface reconstruction error as the number of points in the point cloud N Γ increases. This structure has been used in combination with RBFs for hole-filling in [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF] and anisotropic basis functions for representing sharp corners in [START_REF] Quynh Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF].

Another approach is to construct a uniform grid of points to control the implicit function.

Unlike the aforementioned approaches, the distribution of points is decoupled from the resolution of the point cloud. As a result, deformations to the geometric shape can be represented without loss in accuracy near the surface as shown by Zhao et al. [START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF]. This makes it a popular structure in partial differential equation (PDE) based reconstruction methods that evolve the surface during reconstruction, such as in [START_REF] Tasdizen | Geometric Surface Processing via Normal Maps[END_REF][START_REF] Jakobsen | Variational Volumetric Surface Reconstruction from Unorganized Points[END_REF]. In general, more points representing the implicit function are required to achieve the same level of accuracy to other approaches. As a result, implicit functions defined by a uniform grid are more computationally expensive to solve for in both time and memory usage than the aforementioned approaches, as experienced by Sibley and Taubin [START_REF] Sibley | Vectorfield Isosurface-Based Reconstruction from Oriented Points[END_REF], but can be reduced by a GPU-based multigrid approach as implemented by Jakobsen et al. [START_REF] Jakobsen | Variational Volumetric Surface Reconstruction from Unorganized Points[END_REF].

Formulations for signed distance function approximation

The signed distance function (SDF) presents an ideal candidate for implicit surface reconstruction and geometric non-interference constraints. It is known that the zero level set of the SDF is a smooth representation of the points in a point cloud, and its gradient field is a smooth representation of the normal vector field from the normal vectors in a point cloud. As a result, many formulations to approximate the SDF have been researched for implicit surface reconstruction. We note that there exists other methodologies, such as wavelets [START_REF] Manson | Streaming Surface Reconstruction Using Wavelets[END_REF] and a Fast Fourier Transform (FFT) based method [START_REF] Kazhdan | Reconstruction of Solid Models from Oriented Point Sets[END_REF], that fit a smooth indicator function instead, but are less applicable for non-interfernce constraints where a measurement of distance is desired. We identify four categories that approximate the SDF in some way: explicit formulations, interpolation formulations with RBFs, PDE-based formulations, and energy minimization formulations.

Explicit formulations

Explicit formulations use the data defined in the point cloud to define linear approximations to the SDF. These formulations then apply smoothing to these linear approximations in order to define the level set function. Risco et al. [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF] present the simplest approach which uses the nearest edge and normal vector to define the function explicitly. The resultant constraint function is piecewise continuous but non-differentiable at points where the nearest edge switches. Belyaev et al. [START_REF] Belyaev | Signed Lp-distance fields[END_REF] derive a special smoothing method for defining signed L p -distance functions, which is a continuous and smooth transition between piecewise functions. Hicken and Kaur [START_REF] Hicken | An Explicit Level-Set Formula to Approximate Geometries[END_REF] use modified constraint aggregation methods to define the function in a smooth and differentiable way. Upon the investigation of Hicken and Kaur, the signed L p -distance functions give poor approximations of the surface. Additionally, Hicken and Kaur's formulation is shown to increase in accuracy as the data in the point cloud, number of points N Γ , increases. We identify Hicken and Kaur's explicit formulation as a good candidate for enforcing non-interference constraints, as it is continuous and differentiable with good accuracy.

Given an oriented point cloud is a set of ordered pairs

P = {(p i , ⃗ n i ) : i = 1, . . . , N Γ }
, where p i is the location of the points sampled over the geometric shape, and ⃗ n i are the unit normal vectors at p i , the level set function by Hicken and Kaur [START_REF] Hicken | An Explicit Level-Set Formula to Approximate Geometries[END_REF] is given as

ϕ H (x) = N Γ i=1 d i (x)e -ρ(∆ i (x)-∆ min ) N Γ j=1 e -ρ(∆ j (x)-∆ min ) , ( 1
)
where d i (x) is the signed distance to the hyperplane defined by the point and normal vector pair in the point cloud (p i , ⃗ n i ), ∆ i (x) is the Euclidean distance to the point p i , ∆ min is the Euclidean distance to the nearest neighbor, and ρ is a smoothing parameter. To improve accuracy, Hicken and Kaur suggest modifications to make the linear approximation in to a quadratic approximation by using the principal curvatures of the surface. Unless readily provided by a smooth geometric representation, the principal curvatures must be approximated from the point cloud, similar to the approximation method by Tang and Feng [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF]. To reduce the computational complexity, Hicken and Kaur suggest only evaluating the k-nearest neighbors, since the basis weights exponentially decay with distance. However, using the k-nearest neighbors will remove the function's differentiability as the set of k-nearest neighbors changes. As a result, their formulation must scale in computational complexity by O(N Γ ) to be differentiable. While not originally purposed for geometric non-interference constraints, the formulation by Hicken and Kaur is on par in computational complexity with current non-interference constraint formulations [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF][START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF][START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF]. No calculations are required to set up the function because the equation is explicitly defined by the data in the point cloud, but it is susceptible to the various complexities of point clouds, such as noise and outliers.

Interpolation formulations with radial basis functions

Another method to construct the level set function is to solve an interpolation problem given an oriented point cloud P. Because the data points of P always lie on the zero contour, nonzero interpolation points for the implicit function can be defined on the interior and exterior, as originally done by Turk and O'Brien [START_REF] Turk | Modelling with Implicit Surfaces That Interpolate[END_REF]. Radial basis functions (RBFs) are then formulated to interpolate the data. To avoid overfitting, thin-plate splines can be used to formulate the smoothest interpolator for the data, as noted in [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF][START_REF] Carr | Reconstruction and Representation of 3D Objects with Radial Basis Functions[END_REF]. Solving for the weights of a RBF involves solving a linear system, which is often dense and very computationally expensive due to their global support. Turk and O'Brien [START_REF] Turk | Modelling with Implicit Surfaces That Interpolate[END_REF] solve up to 3,000 RBF centers, and improvements by Carr et al. [START_REF] Carr | Reconstruction and Representation of 3D Objects with Radial Basis Functions[END_REF] allow up to 594,000 RBF centers to be constructed in reasonable time (hours). On top of the significant computational expense, interpolating RBFs have been criticized for having blobby reconstructions [START_REF] Quynh Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF][START_REF] Turk | Modelling with Implicit Surfaces That Interpolate[END_REF] which poorly represent sharp features in the geometric shapes.

PDE-based formulations

Another approach is to construct the level set function as a smooth vector field that smoothly approximates the normal vectors ⃗ n i given by the point cloud P. The vector field is then integrated and fit, usually by a least squares fitting, to make the zero level set fit the point cloud. We classify the methods that solve for the vector field as a solution to a partial differential equations (PDEs) as PDE-based methods. Poisson's method [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF] uses variational techniques to Poisson's equation to construct a vector field. Improvements to this method add penalization weights to better fit the zero contour to the point cloud in [START_REF] Kazhdan | Screened Poisson Surface Reconstruction[END_REF]. Tasdizen et al. [START_REF] Tasdizen | Geometric Surface Processing via Normal Maps[END_REF] prioritize minimal curvature and minimal error in the vector field by solving a set of coupled second order PDEs to derive their level set function. Zhao et al. [START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF] use the level set method, originally introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], for surface reconstruction, with the advantage of modeling deformable shapes. In the aformentioned PDE-based methods, the setup for the implicit function reduces to solving a PDE by time-stepping [START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF][START_REF] Tasdizen | Geometric Surface Processing via Normal Maps[END_REF] or a sparse linear system [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF][START_REF] Kazhdan | Reconstruction of Solid Models from Oriented Point Sets[END_REF] in the case of Poisson's equation. Kazhdan et al. [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF] note that care should be taken when choosing a smoothing filter for the normal field defined by ⃗ n i , especially for nonuniformly sampled points. In the analysis done by Calakli and Taubin [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF], they found that Poisson's method often over-smooths some surfaces. We also note that solutions to PDEs are more difficult to implement than other methods in practice.

Energy minimization formulations

Another methodology is to solve an optimization problem that minimizes some energy function with respect to the values of the basis function directly. The smooth signed distance (SSD) surface reconstruction method [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF] minimizes an energy function with three terms. Minimizing these three terms maximizes smoothness and minimizes the approximation error of the zero level set and the gradient field to the data in P, all in a least squares sense. Alternative forms, such as in [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF][START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF], propose a different energy term to this formulation, which does a direct least squares fit to the approximate signed distance function. We perform a more thorough discussion of the four energy terms in Section 3, as our method also poses an energy minimization problem.

The energy minimization problem posed by these papers is a well-posed unconstrained quadratic programming (QP) problem. The solution to these unconstrained QP problems reduces to the solution of a linear system. Making use of hierarchical structures, such as octrees, and compactly supported basis functions, the linear system is sparse and recursively solved at increasing depths of the structure. These advantages allow for fast solutions on the order of minutes as reported by [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF][START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF]. It should be noted that the time and space (memory) consumed by hierarchical approaches grows exponentially with the depth of the octree, so many implementations limit the depth up to 11. The resultant number of control points in Tang and Feng [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF] are on the order of 10 6 .

Summary

We note that interpolation formulations with RBFs, PDE-based formulations, and energy minimization formulations are different approaches to the same problem of approximating the SDF. The primary differences lie within the derivation and implementation of such methods. The energy minimization formulation by Calakli and Taubin [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF] performs a least squares fit to the data in the point cloud. Thin-plate spline RBFs are an exact solution to the same energy minimization formulation to interpolate the data and maximize smoothness, as derived by [START_REF] Buhmann | Radial Basis Functions: Theory and Implementations[END_REF]. The two-step energy minimization formulation by Sibley and Taubin [START_REF] Sibley | Vectorfield Isosurface-Based Reconstruction from Oriented Points[END_REF] follows the same approach as PDE-based methods, where a vector field is solved for and then a least squares fit is done to fit the surface. We recommend the interested reader to Calakli and Taubin [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF] who discuss the similarities and differences between SSD and Poisson surface reconstruction methods.

We summarize the context for all the methods in Table 1, highlighting the main differences in their formulation, basis function representation, and distribution of points controlling the function. We note our method is an energy minimization formulation, which uses the same energy terms as Calakli and Taubin [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF], but with a different basis function and different distribution of control points.

Methodology

Signed distance function

We now present our methodology to compute an implicit level set function ϕ by approximating the signed distance function (SDF) of a geometric shape. We assume that the geometric shape partitions its neighboring space into a feasible region and an infeasible region as shown in Fig. 2. Our goal is to generate a level set function such that the zero contour of the function approximates the boundary between the feasible and infeasible regions, i.e., the geometric shape. We also require that evaluating the implicit function at any point in the domain of interest will determine if the point Table 1: Signed distance approximations and their basis function representations.

Formulation

Basis Function Point Distribution [START_REF] Hicken | An Explicit Level-Set Formula to Approximate Geometries[END_REF] Explicit Exponential On-surface [START_REF] Turk | Modelling with Implicit Surfaces That Interpolate[END_REF] Interpolating Thin-plate RBF On-& off-surface [START_REF] Carr | Reconstruction and Representation of 3D Objects with Radial Basis Functions[END_REF] Interpolating Polyharmonic RBF On-& off-surface [START_REF] Quynh Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF] Interpolating Anisotropic RBF On-& off-surface [START_REF] Tasdizen | Geometric Surface Processing via Normal Maps[END_REF] PDE-based Polynomial Uniform [START_REF] Zhao | Fast surface reconstruction using the level set method[END_REF][START_REF] Jakobsen | Variational Volumetric Surface Reconstruction from Unorganized Points[END_REF] PDE-based Linear Uniform [START_REF] Michael Kazhdan | Poisson Surface Reconstruction[END_REF] PDE-based Quadratic Octree [START_REF] Kazhdan | Screened Poisson Surface Reconstruction[END_REF] PDE-based Quadratic B-splines Octree [START_REF] Sibley | Vectorfield Isosurface-Based Reconstruction from Oriented Points[END_REF] Energy minimization Linear Uniform [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF] Energy minimization Linear Octree [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF][START_REF] Pan | Phase-Field Guided Surface Reconstruction Based on Implicit Hierarchical B-Splines[END_REF][START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF] Energy minimization Quadratic B-splines Octree Our method Energy minimization Cubic B-splines Uniform is located on the boundary or within one of the two regions, as indicated by the signed distance of the point from the boundary. We follow the convention of denoting distances in the feasible region as positive and those in the infeasible region as negative. The signed distance function in the neighborhood N ⊂ R n of a point on the geometric shape can then be defined as

d Γ (x) = +D(x) if x ∈ Ω -D(x) if x ∈ N \ Ω (2)
where D : R n → R ≥0 measures the shortest distance of a point x to the boundary Γ. The local feasible region Ω ⊂ N and the local infeasible region N \ Ω are separated by the local boundary Γ l = Γ ∩ N within the neighborhood, and Γ l ⊂ Ω. Note that n = 2 or n = 3 for geometric shapes: n = 2 implies Γ is a curve in two dimensions while n = 3 implies Γ is an orientable surface in three dimensions. We assume that Γ is always a connected set. We make no assumptions on the surface or curve being open or closed. However, we assume that Γ does not contain the boundary points or curves if the curve or surface is open to ensure the existence of a neighborhood where the definition of d Γ is valid. We also note that for closed geometric shapes, the definition of a local neighborhood is not necessary, as the feasible and infeasible regions can be simply defined as the inside and outside of the closed boundary Γ (see Fig. 1), or vice versa. This is identical to the standard definition of the signed distance function.

B-spline functions

Our desired level set function is a smooth approximation to the signed distance function of a geometric shape and is defined as a mapping ϕ : V ⊂ R n → R, where V is the space where we wish to evaluate the level set function as a non-interference constraint function during optimization. This means that the zero level set S = {x : ϕ(x) = 0} implicitly approximates the given geometric shape. To achieve such an approximation, we utilize tensor product B-splines. A B-spline volume

Γ ! Γ Ω 𝒩 ∖ Ω + - 𝒩(𝒑 " ) 𝒑 "
Figure 2: A 2D open funnel partitions the neighborhood of a point p i into a feasible region Ω and an infeasible region N \ Ω.

P : [0, 1] 3 → R 3 is defined as

P(u, v, w) = n i ,n j ,n k i,j,k=0 C i,j,k B i,d 1 (u)B j,d 2 (v)B k,d 3 (w), (3)
where (u, v, w) ∈ [0, 1] 3 are the normalized parametric coordinates, B i,d 1 (u), B j,d 2 (v), B k,d 3 (w) are the B-spline basis functions of degrees d 1 , d 2 , d 3 in the i, j, k directions respectively, and C i,j,k are the control points that form the (n i + 1) × (n j + 1) × (n k + 1) control net in the physical coordinate system. Equation ( 3) is essentially a tensor product, and hence P is also called a tensor product B-spline. A B-spline volume maps a volumetric space in the parametric coordinate system (u, v, w) to the physical coordinate system (x, y, z) by translating and deforming the volumetric space according to the control net C i,j,k . The volumetric space we wish to represent in the physical coordinate system is a rectangular prism V. This is the physical space where geometric non-interference constraints need to be evaluated; therefore, we refer to V as the domain of interest. We require that V encompasses the minimum bounding box for a given point cloud. The minimum bounding box is the smallest closed box in R n that contains the input point cloud representing a geometric shape. We uniformly space the control points C i,j,k across V and consider them to be constant. The domain of interest is not always a cube; therefore, parametric coordinates may be scaled differently along different directions. To compensate for this, some directions may contain more control points than others depending on the dimensions of V.

The B-spline basis functions B i,d 1 (u), B j,d 2 (v), and B k,d 3 (w) are generated by the de Boor's recursion formula [START_REF] De | On calculating with B-splines[END_REF]. The formulas for all three directions are identical, and B i,d 1 (u) along the i direction is computed by recursion

B i,0 (u) = 1 if t i ≤ u < t i+1 0 otherwise , B i,k (u) = u -t i t i+k -t i B i,k-1 (u) + t i+k+1 -u t i+k+1 -t i+1 B i+1,k-1 (u), (4)
where t i denotes the knots in the i direction. The basis functions corresponding to C i,j,k provides support only for (u, v, w)

∈ [u i , u i+d 1 +1 ] × [v j , v j+d 2 +1 ] × [w k , w k+d 3 +1
]; thus the basis functions are sparse. This also means that the number of nonzero terms in the summation of Eq. ( 3) is proportional to the degrees d 1 , d 2 , and d 3 .

We define our B-splines with standard uniform knot vectors and uniformly spaced control points C i,j,k across V. This makes the mapping (u, v, w) → (x, y, z) a linear, one-to-one relationship with (u, v, w) ∈ [0, 1] 3 spanning the entire volumetric space of V. Thus, ∂u ∂x , ∂v ∂y , and ∂w ∂z are constants that depend only on the dimensions of the rectangular prism V. Since we are using a standard uniform knot vector, it should be noted that the control points must lie beyond V in order for the mapping to hold for all points in V. With these assumptions, we apply our target function ϕ to both sides of Eq. ( 3) to obtain

ϕ(x, y, z) = n i ,n j ,n k i,j,k=0 C (ϕ) i,j,k B i,d 1 (u(x))B j,d 2 (v(y))B k,d 3 (w(z)), (5)
where u(x), v(y), and w(z) map the physical coordinates back to normalized parametric coordinates, and C

(ϕ) i,j,k are the values of the function ϕ at the control points. The derivatives with respect to the spatial coordinates and derivatives with respect to the control points are easily derived from this form. The sparsity of the basis functions over the entire domain and the use of standard uniform knot vectors make the computation of ϕ at any given point (x, y, z) in the domain highly efficient. Moreover, ϕ becomes continuously differentiable with easily computed higher order derivatives for appropriately defined B-spline degrees.

Note that while computing the level set function ϕ for a curve in two dimensions, V is a rectangle, P reduces to a B-spline surface, and we omit terms along the k direction in Eq. ( 5). In all of the remaining discussion, we assume n = 3 with the geometric shape being a surface in three dimensions.

Energies for B-spline fitting

The core of our methodology lies in computing appropriate C (ϕ) i,j,k values on the control net so that the level set function ϕ approximates the signed distance function with favorable properties for gradient-based optimization. Note that since our approximation uses B-splines, the resulting function will already be smooth and fast-to-evaluate. Therefore, the remaining task is to formulate an approach for reliably estimating C (ϕ) i,j,k using the information available from an oriented point cloud.

An oriented point cloud is a set of ordered pairs P = {(p i , ⃗ n i ) : i = 1, . . . , N Γ }, where p i are the physical coordinates of the points sampled over the geometric shape, and ⃗ n i are the unit surface (or curve) normal vectors at p i oriented towards the infeasible region. Our method always requires an oriented point cloud as its input. However, we note that in cases where only a point cloud without normal information is available, Principal Component Analysis (PCA) along with a Minimum Spanning Tree (MST) algorithm can be used for estimating normals and their orientation [START_REF] Hoppe | Surface Reconstruction from Unorganized Points[END_REF]. Edge-Aware Resampling (EAR) [START_REF] Huang | Edgeaware point set resampling[END_REF] is another method that can be used for generating noise-free normals that also preserves sharp features.

We calculate C (ϕ) i,j,k values by minimizing an energy function consisting of multiple energies. The terms in the energy function are adopted from existing surface reconstruction methods [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF][START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF][START_REF] Pan | Phase-Field Guided Surface Reconstruction Based on Implicit Hierarchical B-Splines[END_REF][START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF]. Since the zero contour of our desired level set function ϕ should approximate the geometric shape represented by the point cloud, it is straightforward to see that we should minimize energies to approximately satisfy ϕ(p i ) = 0 and ∇ϕ(p i ) = -⃗ n i . Hence we first define energies

E p = 1 N Γ N Γ i=1 ϕ(p i ) 2 and E n = 1 N Γ N Γ i=1 ∥∇ϕ(p i ) + ⃗ n i ∥ 2 , ( 6
)
where E p estimates the approximation error as the average of squared distances of the point cloud from the zero contour of ϕ, and E n measures the average of squared alignment errors of the level set function's gradient when compared to the negative of the unit normal vectors in the point cloud.

Note that we take the negative of the normals (oriented toward the infeasible region) from the point cloud since we want distances given by ϕ to be positive inside the feasible region. Minimizing E p forces the zero contour of ϕ to pass through all the points in the point cloud, and minimizing E n tries to orient the function's direction of steepest increase ∇ϕ along the normal to the geometric shape while pointing toward the feasible direction, both in the least squares sense. Minimizing E n is important since the derivatives of the exact signed distance function d Γ on the boundary of a geometric shape is along the normal to the boundary, and d Γ always satisfies the eikonal equation, i.e., ∥∇d Γ ∥ = 1.

If we perform a direct minimization of energies E p and E n , the resulting function attempts to accurately fit the point data on the geometric shape, and since these energies do not control the behavior of ϕ away from the geometric shape, it could create superfluous zero contours away from the point cloud as reported in previous studies [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF]. To overcome this issue, we define the regularization energy

E r = 1 |V | V ∇ 2 ϕ(x) 2 F dV, (7)
where ∇ 2 ϕ(x) is the Hessian matrix of ϕ evaluated at x, ∥•∥ F represents the Frobenius norm, and |V | = V dV is the total volume of V. The regularization energy E r is interpreted as the aggregate curvature of ϕ over the entire volumetric space of V. The minimization of E r smooths the function ϕ since forcing the Hessian to be zero forces the variations in the gradient field ∇ϕ to a minimum. Since the gradient of ϕ is approximately aligned with the unit normals on the point cloud when minimizing E n , trying to maintain a constant ∇ϕ by minimizing E r also helps satisfy the eikonal equation ∥∇ϕ∥ = 1 for points further away from the point cloud. We evaluate the integral in E r using the B-spline control points C i,j,k lying inside V as quadrature points with unit quadrature weights. Therefore, the regularization energy is approximated as a discrete sum is given by

E r = 1 |V | V ∇ 2 ϕ(x) 2 F dV ≈ 1 N N i=1 ∇ 2 ϕ(x i ) 2 F , ( 8
)
where N is total number of quadrature points, typically about the same as the number of control points N cp = (n i + 1) × (n j + 1) × (n k + 1). Some surface reconstruction techniques employ another energy term E d , which attempts to fit the signed distance function over the entire domain V. However, minimizing this energy was found to create overfitting issues and produce high frequency oscillations in the level set function ϕ in our investigation and previous studies [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF]. As a result, we neglect this energy in our formulation. Nevertheless, we present it here for the sake of completeness. The signed distance energy is given by

E d = 1 N N i=1 (ϕ(x i ) -d Γ (x i )) 2 , ( 9
)
where signed distances d Γ (x i ) are evaluated at the control points within V (same as quadrature points in E r ). The signed distances d Γ (x i ) can be approximated using distances to the nearest neighbor in the point cloud and its normal [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF], or by evaluating the explicit equation ( 1). Note that minimizing E r can act as a regularization to avoid overfitting caused by E d but careful weighting of the four energies according to the geometric shape is necessary.

Final energy minimization problem

Finally, we define the total energy function f as

f = 1 N Γ N Γ i=1 ϕ(p i ) 2 Ep + λ n N Γ N Γ i=1 ∥∇ϕ(p i ) + ⃗ n i ∥ 2 λnEn + λ r N N i=1 ∥∇ 2 ϕ(x i )∥ 2 F λrEr , ( 10
)
where λ n and λ r are the relative penalization weights for E n and E r with respect to E p . The energy minimization problem that yields the desired level set function ϕ is then given by

minimize f = E p + λ n E n + λ r E r with respect to C (ϕ) i,j,k . (11)
Note that the function values at the control points C (ϕ) i,j,k directly affect E p , E n , and E r through the definition of ϕ using B-splines (see Eq. ( 5)). If the geometric shape is a curve in two dimensions, then the optimization variables are C (ϕ) i,j . The choice of penalization weights is not obvious. Penalization weights may require tuning on a case-by-case basis depending on the geometric shape. In general, we recommend λ n ∼ 10 -2 and λ r ∼ 10 -4 based on the parameter study presented in Sec. [START_REF] Criado Risco | Gradient-based Wind Farm Layout Optimization With Inclusion And Exclusion Zones[END_REF].

We provide a summary of our methodology in Algorithm 1 for geometric shapes that are surfaces in three dimensions. The algorithm is easily adapted for curves in two dimensions by simply omitting terms along the k direction.

Algorithm 1 A scalable and differentiable geometric non-interference constraint formulation

1: Discretize the given geometric shape into an oriented point cloud P. 2: Define the domain V for non-interference constraint evaluation, where V is a closed box in R 3 and contains the minimum bounding box of P. 3: Select the appropriate numbers of control points n i , n j , and n k along each direction, and define the corresponding standard uniform knot vectors and spatially uniform control points C i,j,k . Note that some of the control points may lie outside V. 4: Select appropriate weights λ n and λ r , and solve the energy minimization problem [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF] to obtain C (ϕ) i,j,k . 5: Evaluate the geometric non-interference constraint(s) during each optimization iteration using the ϕ given by Eq. ( 5) and optimized C (ϕ) i,j,k from Step 4.

Implementation details

We initialize C (ϕ) i,j,k for the energy minimization problem [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF] by evaluating the explicit equation ( 1) at each control point C i,j,k . This initialization gives an overall good initial guess for the signed distance near the geometric shape, hence it results in smaller E p and E n values. While the initialized B-spline function ϕ is always differentiable for degree two or more, points of non-differentiability in the exact signed distance function will create regions of high curvatures in ϕ. Hence, this initialization may result in a large E r term, even for smooth geometric shapes.

The minimization of E r is thus necessary to smooth these regions of high curvatures, although compromising the accuracy at representing the exact SDF in these regions. Thus, λ r should be large enough to enable smoothing of high curvatures that exist in the exact signed distance initialization but small enough so that it does not induce a large error in representing signed distances near the geometric shape. This is essential for better convergence in the overarching optimization problem when solved using a gradient-based algorithm.

In our implementation, we define the B-spline domain V by extending the minimum bounding box for the point cloud along its diagonal by 15%. The choice of 15% is purely empirical, and it can be lower or higher depending on the optimization problem requirements. The additional margin allows optimization algorithms to evaluate the constraint function at locations that are away from the geometric shape.

Higher order B-splines are computationally expensive. However, they improve the local degrees of freedom and the ability to control the function. In our experiments with the Stanford Bunny model shown in Fig. [START_REF] Andrew | Coupled wind turbine design and layout optimization with nonhomogeneous wind turbines[END_REF], we found no significant reduction in error for B-spline degrees higher than three. Hence, we recommend cubic B-splines for reasonable accuracy and computational efficiency. It is also often computationally expensive to express a non-interference constraint for a body as separate optimization constraints ϕ(x (i) ) ≥ 0 for a large number of points x (i) on the body. Instead, we recommend implementing a smooth minimum or maximum function such as KS-aggregation [START_REF] Kreisselmeier | SYSTEMATIC CONTROL DESIGN BY OPTIMIZING A VECTOR PERFORMANCE INDEX[END_REF] to reduce the number of constraints in the optimization problem.

Lastly, we note that the energy minimization problem [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF] is inherently an unconstrained quadratic programming (QP) problem, and it can be reduced to the solution of a sparse linear system in practice. Since our focus is on creating a scalable and smooth formulation for representing geometric non-interference constraints for gradient-based optimization, we do not delve deeper into the numerical formulation of this linear system. We model our energy minimization problem using the open source OpenMDAO [START_REF] Gray | OpenMDAO: An Open-Source Framework for Multidisciplinary Design, Analysis, and Optimization[END_REF] optimization framework and utilize the general-purpose optimization software SNOPT [START_REF] Gill | SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization[END_REF] to perform energy minimization for numerical studies presented in the next section.

Numerical Study

This section presents the results of various numerical studies using our formulation. We begin by studying our method using simple two-dimensional geometric shapes in subsection 4.1. In subsection 4.2, we investigate the dependence of our method on various weights using the Stanford Bunny dataset. Subsection 4.3 then compares our method to previous non-interference constraint formulations using the Stanford Bunny, and other surface reconstruction methods using three datasets from the Stanford 3D Scanning Repository. We demonstrate the accuracy of our method in representing geometric shapes for aircraft design optimization in subsection 4.4. We conclude the section by demonstrating the application of our method by solving a medical robot design optimization problem with non-interference constraints in subsection 4.5.

We implement the proposed method in a Python environment, and run all experiments on a desktop with an 8-core Ryzen 7 @ 3.6 GHz processor and 32 GB of RAM. We do not implement multi-threading or parallelization with GPUs in any of our numerical experiments.

Investigations using simple geometric shapes in two dimensions

We begin by applying our formulation to curves in two dimensions. Because our formulation is generic, no modifications are required to Eq. ( 10), and the terms in the k direction are simply ignored for 2D geometric shapes.

For 2D curves, the isocontours from the level set function (LSF) ϕ can be readily visualized to facilitate a better understanding of the function both near and far from the curve. Figure 3 visualizes the isocontours of the initialized and energy minimized LSF for a rectangle using our formulation. We initialize C (ϕ) i,j using the explicit equation [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF]. Neglecting the sharp corners in this example, the contours of the initialized function closely match the exact signed distance function (SDF). Thus, the explicit method provides an excellent approximation of the SDF. We observe that the contours of the LSF are more rounded near the corners after energy minimization. Minimizing E r smooths sharp corners on all isocontours, however, not to a degree that compromises E n and E p near the zero contour. We note that E n and E p have less influence compared to E r on the isocontours corresponding to 1 and 2, hence these contours are even more rounded. A LSF representing multiple geometric shapes may also be obtained using a single B-spline. Figure 4 shows the exact SDF and a one-dimensional slice of our energy minimized LSF ϕ along the x axis for a domain containing multiple circles. Notice that the non-differentiable points in the exact SDF can lie inside or outside of a geometric shape. This example illustrates how our energy minimization formulation balances the tradeoff between minimizing the curvature of ϕ and maximizing the accuracy at representing the SDF near points of non-differentiability and high curvature. We see that our LSF ϕ poorly approximates the SDF near points of non-differentiability and high curvature. However, in regions without any non-differentiabilities or high curvatures, the zero level set preserves a good approximation to the exact SDF. For the remainder of the numerical results section, we only consider a single geometric shape within the domain of interest V, because the error of our formulation increases with the minimum bounding box diagonal.

Investigations using a complex geometric shape in three dimensions

We use the well known Stanford Bunny scanned dataset (shown in Fig. 5) to analyze our formulation's performance on three-dimensional geometric shapes. This dataset contains a large oriented point set which we consider as an exact surface. We coarsely sample this point set and apply our method, measuring the accuracy of the resultant energy minimized LSF. The Stanford Bunny contains small scale features, sharp corners, flat surfaces, and smooth surfaces which will test the accuracy of our formulation in representing different geometric features. The sampled point set is free of noise, missing data, and nonuniformity, which are challenges not investigated in this paper. We use the root-mean-squared (RMS) error and max error to evaluate the accuracy of our energy minimized LSF in approximating the signed distance function. The errors are normalized by the minimum bounding box diagonal L to ensure that they are independent of the size of a geometric shape. This allows for a common metric for comparing accuracies across different geometric shapes. The errors are defined as

RMS error = 1 L Ne i=1 (ϕ(x i ) -d Γ (x i )) 2 N e , and (12)
max error = max i=1,2,...,Ne

1 L |ϕ(x i ) -d Γ (x i )|, (13)
where N e is the number of points x i used to calculate the error, and the signed distances are approximated using the explicit equation ( 1). The on-surface error is evaluated on the geometric shape Γ where the true value is zero. The off-surface error is evaluated on points that are near but do not lie on Γ. To acquire these sample points, we take the original sample points and move them in the direction of the normal vectors by specified distances.

The energy minimization problem has two different resolution scales: (1) resolution of the input dataset, and (2) resolution of the B-spline control grid. The energy minimized LSF's ability to approximate the signed distance is best when both resolutions are very fine. Unlike hierarchical structures or explicit methods, the control grid resolution for our method is independent of the input dataset resolution. As a result, an important consideration in our method is the decision of how many control points with which to represent the implicit function. Table 2 tabulates the on-surface errors for our energy minimized function, varying the two resolution scales. A total of nine experiments were ran, using three different resolutions for each resolution scale. It was found that only increasing N Γ was not correlated to any increase in the computation time. Thus, we only consider the average time to solve across the three experiments with the same N cp . In terms of the on-surface error, increasing both resolutions correlates to a decrease in both RMS and maximum error. In terms of the time to solve the energy minimization problem (setup time), increasing the number of control points N cp increases the number of design variables in the minimization problem [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF], consequently increasing the time to solve for C (ϕ) i,j,k .

Table 2: The relative on-surface error for our method applied to the Stanford Bunny model. Penalization weights used were λ n = 10 -2 , and λ r = 5 × 10

N Γ = 5k RMS 1.0 × 10 -3 7.5 × 10 -4 6.7 × 10 -4 Max 4.7 × 10 -3 3.7 × 10 -3 3.3 × 10 -3 N Γ = 25k RMS 8.2 × 10 -4 5.4 × 10 -4 4.6 × 10 -4 Max 3.7 × 10 -3 3.3 × 10 -3 2.7 × 10 -3 N Γ = 64k RMS 7.8 × 10 -4 5.1 × 10 -4 4.2 × 10 -4 Max 6.1 × 10 -3 5.1 × 10 -3
2.9 × 10 -3 Average setup time 11:19 27:47 42:13 (min:sec) While we do not propose an exact method for selecting the penalization weights λ n and λ r , we provide a fixed point parameter study on each weight. Figure 6 shows the resulting errors from varying each penalization weight about the fixed point λ n = λ r = 1 using the Stanford Bunny dataset. The study on λ r shows that small values (λ r < 1) have very little effect on the RMS error, and large values (λ r > 1) significantly reduce the energy minimized function's accuracy. The study on λ n suggests that for a given geometric shape, there exists an optimum value of λ n that minimizes the energy minimized function's error. In all studies, we observed an increase in the minimization time as the corresponding weight increased (not visualized in the figure). These observations lead us to recommend the use of λ n ∼ 10 -2 and λ r ∼ 10 -4 for reasonable accuracy and fast setup time for the energy minimized function.

RMS Error

Off-surface (±0.01) On-surface The ability of our function to represent nonzero level sets of the Stanford Bunny is visualized in Fig. 7. The level sets form good approximations of the offset surfaces, with a maximum relative distance error of 9.8 × 10 -3 . In these visualizations, we observe the region of highest error to be near the neck and feet of the model, where edges and corners exist. Most notably, the 0.005 and 0.01 level sets remove the ears of the model, despite them being in the exact SDF representation. As a thin feature, the removal of the ears in the 0.005 and 0.01 level sets is consistent with similar observations by Tang and Feng [START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF].

Comparison to other methods using complex geometric shapes in three dimensions

We show the computation time and accuracy of our method compared to alternative noninterference constraint formulations for gradient-based methods in Fig. 8, varying the sample size N Γ of the Stanford Bunny. We observe that the method presented by Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF] and the explicit method presented by Hicken and Kaur [START_REF] Hicken | An Explicit Level-Set Formula to Approximate Geometries[END_REF] scale in computational complexity with O(N Γ ), while our method scales independently of N Γ . We note that formulation from Lin et al. is not an attempt at approximating the signed distance function, thus is neglected from the RMS error comparisons. In terms of on-surface error, the explicit method has a steady decay in RMS error with respect to increasing N Γ , suggesting a power law relationship. Our method has a similar decay up to N Γ = 10 4 , where the RMS error decays more slowly for larger N Γ > 10 4 . Similarly, the off-surface RMS error of the explicit method steadily decays for both the ±0.005 and ±0.01 contours, and our method decays until N Γ = 10 4 . For N Γ > 10 4 , the off-surface error of our method decays slowly. Our method's ±0.01 contours have significantly more error than the ±0.005 contours, while the explicit method has similar error for both sets of isocontours. For both on-surface and off-surface error, our method performs better in terms of accuracy up until the ±0.005 contours and N Γ < 2 × 10 4 . From this information, we conclude that the explicit method will outperform in terms of accuracy and underperform in terms evaluation time compared to our method for most very finely sampled geometries. We note that our method can achieve better accuracy than shown in Fig. 8 by a tradeoff in setup time as shown in Table 2, and the explicit method requires a noise-free, uniform sampling to achieve the presented results, which is not always feasible. We apply our method using two additional scanned datasets from the Stanford 3D Scanning Repository. Table 3 records the results of the on-surface error of our method, as well as the reported on-surface error from four notable surface reconstruction methods for necessary context. The methods are smooth signed distance (SSD) reconstruction [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF], Multi-Level Partition of Unity (MPU) [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF], wavelets [START_REF] Manson | Streaming Surface Reconstruction Using Wavelets[END_REF], and screened Poisson (SP) [START_REF] Kazhdan | Screened Poisson Surface Reconstruction[END_REF]. The results for the surface reconstruction methods were obtained from [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF][START_REF] Pan | Phase-Field Guided Surface Reconstruction Based on Implicit Hierarchical B-Splines[END_REF][START_REF] Tang | Multi-scale surface reconstruction based on a curvatureadaptive signed distance field[END_REF][START_REF] Kazhdan | Screened Poisson Surface Reconstruction[END_REF] and were not reproduced in our investigation, resulting in missing data in the table. Of the three scanned datasets, our energy minimized LSF maintains on-surface RMS error and max error on the same order of magnitude compared to the four other methods.

Table 3: Reported on-surface error of surface reconstruction methods and our method for three benchmarking datasets. We reconstruct using N cp = 25, 000, λ n = 10 -2 , and λ r = 5 × 10 -4 .

Model

SSD [START_REF] Calakli | SSD: Smooth Signed Distance Surface Reconstruction[END_REF] MPU [START_REF] Ohtake | Multi-Level Partition of Unity Implicits[END_REF] Wavelets [START_REF] Manson | Streaming Surface Reconstruction Using Wavelets[END_REF] SP [START_REF] Kazhdan | Screened Poisson Surface Reconstruction[END_REF] Our method

Accuracy of our method for aircraft design optimization

We now apply our formulation to a number of geometric shapes involved in novel aircraft design. Aircraft design optimization is a long standing problem and has been the subject of recent interest in problems involving geometric non-interference constraints, e.g., the layout optimization of air cargo [START_REF] Brandt | The air cargo load planning problem[END_REF] and aerodynamic shape optimization [START_REF] Brelje | Flexible Formulation of Spatial Integration Constraints in Aerodynamic Shape Optimization[END_REF]. To enable gradient-based design optimization involving these constraints, a new generic method is required to represent numerous components within an aircraft's design. We recognize the potential for our formulation and demonstrate its capabilities by conducting an experiment.

In this experiment, we apply our formulation and quantify the resultant errors of five geometric shapes commonly associated with aircraft design. The geometries we model include a fuselage and a wing from a novel electric vertical take-off and landing (eVTOL) concept vehicle [START_REF] Silva | VTOL Urban Air Mobility Concept Vehicles for Technology Development[END_REF], a human avatar [START_REF] Matthew P Reed | Developing and implementing parametric human body shape models in ergonomics software[END_REF], a luggage case, and a rectangular prism representing a battery pack within the wing. A visualization of these components in a feasible design configuration is illustrated in Fig. 9.

Table 4 tabulates the on-surface error of the energy minimized LSF for each geometry. We observe that the smallest relative on-surface error is of the smooth fuselage shape, while the largest relative error is of the human avatar. We note from this example that geometries with features reasonably proportioned to their minimum bounding box diagonal are easier to represent using our method, hence our formulation poorly represents small scale features (e.g. hands and feet) of the human avatar while it can represent the smooth fuselage very well. We observe that the bounding boxes of the fuselage, wing, and battery pack are poorly proportioned, yet they do not result in an increase of relative error compared to other geometries. However, their longer minimum bounding box diagonals will result in larger absolute errors.

Application to a medical robot design problem

We now apply our method for enforcing geometric non-interference constraints to a medical robot design problem involving concentric tube robots (CTRs). CTRs are composed of two or more long and slender pre-curved tubes made of superelastic materials. They can be designed to reach points in a large region of interest by rotating and translating the tubes relative to each other at their bases. These characteristics make them ideal for minimally invasive surgeries where a surgeon can operate on a small region of interest with high dexterity through actuation at the base.

In the foundational works of Sears and Dupont [START_REF] Sears | A Steerable Needle Technology Using Curved Concentric Tubes[END_REF] and Webster et al. [START_REF] Webster | Toward Active Cannulas: Miniature Snake-Like Surgical Robots[END_REF], expressions for the shape and tip position of the CTR are derived with respect to the robot's geometric and control variables. Bergeles et al. [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF] use these expressions to perform gradient-free optimization of the CTR's geometric and control variables with anatomical constraints. These anatomical constraints, i.e., geometric non-interference constraints, enforce that the CTR does not interfere with the anatomy (e.g., the right ventricle of the heart shown in Fig. 10) during operation. Recent work by Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF] shows that gradient-based optimization enables an efficient and scalable solution to simultaneously optimize the large set of the tube's geometric and control variables while enforcing anatomical constraints. The experiment we now present follows the workflow of Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF], however, using our new formulation for representing the anatomical constraint function.

The presented workflow involves the solution of multiple optimization problems, including an initial path planning problem, and the geometric design and control of the CTR (the 'simultaneous optimization problem' described by Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF]). The path planning problem solves for a parametric 3D curve that represents an optimal collision-free path to the surgical site within the anatomy. Then, points along this path serve as inputs to the geometric design and control optimization of the CTR, which involve a kinematic model of the robot. In both subproblems, the non-interference constraints are enforced by evaluating a discrete set of points along the path or physical CTR to ensure that no points lie outside of the anatomy.

We begin our experiment with an investigation in the heart anatomy which represents the noninterference constraint of the problem. The initial oriented point cloud of the heart is obtained from segmentation and 3D reconstruction by magnetic resonance imaging (MRI) scans. Due to the limited machine accuracy, error introduced by aligning multiple scans, and normal approximation, the oriented point cloud is noisy, nonuniform, and contains poorly oriented normals. We perform a simple and necessary smoothing step on this point cloud as illustrated in Fig. 10. Although less precise at capturing small scale features, the smoothing step assists our method in reconstructing a smooth zero contour for constraint representation. The smooth representation has relative errors 3.1×10 -3 (RMS) and 1.9×10 -2 (max) compared to the original noisy representation. The error in our energy minimized function obtained from the smoothed heart model is tabulated in Table 5. We observe that the on-surface RMS and max error of our representation is an order of magnitude less than the error introduced by the smoothing step. This implies that our representation of the smooth model is no worse than the smoothing step itself. We see that our method generates a function with a reliable zero level set of the smooth heart geometry, with an on-surface RMS error of 2.1 × 10 -4 . This error is lower compared to all the other examples in Table 3, and we attribute this to the smoothness of the heart geometry. We also note that the max on-and off-surface absolute errors of our representation are of the same order as the diameter of the CTR itself, typically 0.5-2.0 mm.

We now solve the two optimization subproblems using our energy minimized LSF of the smoothed heart model to enforce the geometric non-interference constraint. In the model from Lin et al. [START_REF] Lin | A Generalized Framework for Concentric Tube Robot Design Using Gradient-Based Optimization[END_REF], the non-interference constraint was imposed using a penalization function g(x), where it was defined as negative for the feasible region, and positive for the infeasible region. In our imple-Table 5: Accuracy of our method in representing a smoothed heart model. Minimum bounding box diagonal is 244.75 mm. The anatomy is sampled at N Γ = 100, 000, and our method uses a control point discretization of (28 × 23 × 37), λ n = 10 -2 , and λ r = 10 The results from this experiment are shown in Table 6, where the number of function evaluations and optimization time are tabulated for each subproblem and non-interference constraint method. The time to solve the energy minimization problem for our method is denoted as the setup time.

Between the two subproblems, the number of function evaluations and optimization time is significantly more for the design subproblem due to the inclusion of the kinematics models. Between the two non-interference constraint methods, we observe a significant decrease in optimization time by using our new method for both subproblems. Even when accounting for the setup time, our method provides a significant speedup for the design subproblem. However, we note that the speedup provided by our method for computationally inexpensive optimization problems, such as the path planning subproblem, may be negated by the setup time to solve for the energy minimized LSF.

For geometries with larger N Γ and more complex optimization problems requiring more function evaluations, we expect the speedup in optimization time to be more pronounced. In this paper, we investigated constraint modeling to prevent interference between geometric shapes in gradient-based optimization. In Sec. 1, we consolidated the terminology used in prior literature and call this category of constraints 'geometric non-interference constraints'. Additionally, we framed the set of optimization problems with geometric non-interference constraints into three groups: layout optimization, shape optimization, and optimal path planning problems. Section 2 reviewed the existing geometric non-interference constraint formulations in gradient-based optimization and contextualized our formulation within the field of surface reconstruction. In Sec.

3, we drew upon ideas from surface reconstruction techniques to construct our constraint formulation. Our formulation is based on an approximation of the signed distance function generated by solving an energy minimization problem for the values of the B-spline control points. Section 4 presented accuracy and scaling studies with our formulation. We also solved a path planning and shape optimization problem using our new formulation.

The contribution of this paper is a new formulation for representing geometric non-interference constraints in gradient-based optimization. This formulation involves a scalable, smooth, and fastto-evaluate constraint function that approximates the local signed distance to a geometric shape. The use of B-spline functions is key to our formulation being scalable, smooth, and fast-to-evaluate. We showed that our formulation achieves a level of accuracy on the same order of magnitude as surface reconstruction methods used in computer graphics. Additionally, our formulation scales better in accuracy, up to a certain limit, and computational time with respect to the number of points sampled on the geometric shape N Γ compared to previous non-interference constraint formulations used by the optimization community. Our resultant computational speed is on the order of 10 -6 seconds per point as measured on a modern desktop workstation, entirely independent of the number of sample points N Γ . The method results in a 78% and 56% speedup in optimization time for a path planning and design subproblem, respectively, for an existing concentric tube robot (CTR) gradient-based design optimization problem.

We identify multiple directions for future work. Adaptive octrees with B-splines can represent small-scale features such as edges and sharp corners more accurately. Using octrees for discretization instead of using a uniform grid can clearly yield faster and more accurate solutions in problems where any of the modeled geometries remain constant during optimization iterations, e.g., the CTR or wind farm layout optimization problems. However, it is worth restating that when geometries evolve during optimization, rediscretizing surfaces using octrees in each optimization iteration becomes unreasonably expensive, and we only recommend a uniform discretization in such cases. Acceleration with multi-threading or graphics processing units (GPUs) is another possible direction for future research.

Figure 3 :

 3 Figure 3: The -1, 0, 1, and 2 contours of the initialized (left) and the energy minimized (right) level set function ϕ.

φ

  Our function φExact SDF dΓ Non-differentiable points

Figure 4 :

 4 Figure 4: The exact SDF (top) and the energy minimized LSF ϕ along a 1D slice (bottom) for multiple geometric shapes.

Figure 5 :

 5 Figure 5: The Stanford Bunny model.

- 4 . 32 ×

 432 31 × 25 40 × 39 × 32 48 × 45 × 33 N cp = 24, 800 N cp = 49, 920 N cp = 71, 280

10 - 4

 4 10 -3 10 -2 10 -1 10 0

Figure 6 :

 6 Figure 6: The RMS errors for on-and off-surface points while varying λ n and λ r about 1. This result uses the Stanford Bunny sampled at N Γ = 25, 000, and a control point grid of 31 × 31 × 26.

Figure 7 :

 7 Figure 7: Isocontours of the energy minimized LSF for the Stanford Bunny. The colors represent the error of the isocontour to the true signed distance value. The isocontours and error are normalized by the minimum bounding box diagonal.

Figure 8 :

 8 Figure 8: Evaluation time per point (left), on-surface root-mean-square error (center), and offsurface root-mean-square error (right) varying the sampling (N Γ ) of the Stanford Bunny model. Our method was applied using a control point grid of (31 × 31 × 26) and λ n = 10 -2 , λ r = 5 × 10 -4 .

3 × 10 - 3

 3103 Max 9.0 × 10 -4 1.9 × 10 -3 2.0 × 10 -3 8.0 × 10 -4 7.6 × 10 -3 Dragon RMS 3.5 × 10 -4 8.0 × 10 -4 1.4 × 10 -3 5.1 × 10 -4 1.4 × 10 -3 Max ... 4.8 × 10 -3 ... 5.1 × 10 -3 1.0 × 10 -2

Figure 9 :

 9 Figure 9: Components necessary for spatial integration in aircraft design optimization. A cross section view is shown for an aircraft fuselage, wing, battery pack, human avatar, and luggage.

Figure 10 :

 10 Figure 10: Preprocessing the raw scanned data (left) to a smooth approximate model (right).

  represent this function with our energy minimized LSF in the form g(x) = -ϕ(x).

Table 4 :

 4 Relative on-surface error of our method on various components of engineering systems. All geometries were sampled at N Γ = 25, 000. The optimization weights λ n = 10 -2 and λ r = 5 × 10 -4 were used.

			Fuselage	Luggage	Wing	Battery pack Human avatar
	Relative error	RMS	7.6 × 10 -5	1.8 × 10 -4	2.5 × 10 -4	4.4 × 10 -4	7.8 × 10 -4
		Max	5.7 × 10 -4	1.8 × 10 -3	1.0 × 10 -3	1.8 × 10 -3	3.8 × 10 -3
	Absolute error	RMS	0.24 cm	0.01 cm	1.28 cm	0.66 cm	0.14 cm
		Max	1.78 cm	0.15 cm	5.29 cm	2.61 cm	0.67 cm
	Discretization		47 × 29 × 29 33 × 44 × 28 29 × 47 × 29 29 × 47 × 29	37 × 32 × 34
	Setup time (min:sec)		17:49	12:41	21:36	12:34	13:09

Table 6 :

 6 Constraint function evaluations and optimization time for the concentric tube robot's path planning and design optimization. The anatomy is represented using N Γ = 1, 842.

		Lin et al. [2]	Our method
	Subproblem	Path planning	Design	Path planning	Design
	Function evaluations	37	35,142	62	20,793
	Optimization time	4.2 sec	3 hr 11 min	0.9 sec	1 hr 24 min
	Setup time	N/A		3 min 51 sec
		5 Conclusion

Acknowledgments

The models from the computer graphics benchmarks are courtesy of the Stanford 3D Scanning Repository (Stanford Bunny, Armadillo, Dragon).

Declarations

The second author was supported by the National Science Foundation under grant no. 1917142.

Replication of results

The energy minimization problem was formulated in Python using OpenMDAO, an open-source optimization framework. The code to generate the Bsplines and solve these models are publicly available on a GitHub repository ( https://github.com/LSDOlab/lsdo noninterference ).