Colin B. Macdonald: Publications

We demonstrate a method for filtering images defined on curved surfaces embedded in 3D. Applications are noise removal and the creation of artistic effects. Our approach relies on in-surface diffusion: we formulate Weickert's edge/coherence enhancing diffusion models in a surface-intrinsic way. These diffusion processes are anisotropic and the equations depend non-linearly on the data. The surface-intrinsic equations are dealt with the closest point method, a technique for solving partial differential equations (PDEs) on general surfaces. The resulting algorithm has a very simple structure: we merely alternate a time step of a 3D analog of the in-surface PDE in a narrow 3D band containing the surface with a reconstruction of the surface function. Surfaces are represented by a closest point function. This representation is flexible and the method can treat very general surfaces. Experimental results include image filtering on smooth surfaces, open surfaces, and general triangulated surfaces.

We introduce a method-of-lines formulation of the closest point method, a numerical technique for solving partial differential equations (PDEs) defined on surfaces. This is an embedding method, which uses an implicit representation of the surface in a band containing the surface. We define a modified equation in the band, obtained in a straightforward way from the original evolution PDE, and show that the solutions of this equation are consistent with those of the surface equation. The resulting system can then be solved with standard implicit or explicit time-stepping schemes, and the solutions in the band can be restricted to the surface. Our derivation generalizes existing formulations of the closest point method and is amenable to standard convergence analysis.

Adaptive stepsize control is a critical feature for the robust and efficient numerical solution of initial-value problems in ordinary differential equations. In this paper, we show that adaptive stepsize control can be incorporated within a family of parallel time integrators known as Revisionist Integral Deferred Correction (RIDC) methods. The RIDC framework allows for various strategies to implement stepsize control, and we report results from exploring a few of them.

Elliptic partial differential equations are important both from application and analysis points of views. In this paper we apply the Closest Point Method to solving elliptic equations on general curved surfaces. Based on the closest point representation of the underlying surface, we formulate an embedding equation for the surface elliptic problem, then discretize it using standard finite differences and interpolation schemes on banded, but uniform Cartesian grids. We prove the convergence of the difference scheme for the Poisson's equation on a smooth closed curve. In order to solve the resulting large sparse linear systems, we propose a specific geometric multigrid method in the setting of the Closest Point Method. Convergence studies both in the accuracy of the difference scheme and the speed of the multigrid algorithm show that our approaches are effective.

We demonstrate a method for removing noise from images or other data on curved surfaces. Our approach relies on in-surface diffusion: we formulate both the Gaussian diffusion and Perona-Malik edge-preserving diffusion equations in a surface-intrinsic way. Using the Closest Point Method, a recent technique for solving partial differential equations (PDEs) on general surfaces, we obtain a very simple algorithm where we merely alternate a time step of the usual Gaussian diffusion (and similarly Perona-Malik) in a small 3D volume containing the surface with an interpolation step. The method uses a closest point function to represent the underlying surface and can treat very general surfaces. Experimental results include image filtering on smooth surfaces, open surfaces, and general triangulated surfaces.

We study spatially partitioned embedded Runge-Kutta (SPERK) schemes for partial differential equations (PDEs), in which each of the component schemes is applied over a different part of the spatial domain. Such methods may be convenient for problems in which the smoothness of the solution or the magnitudes of the PDE coefficients vary strongly in space. We focus on embedded partitioned methods as they offer greater efficiency and avoid the order reduction that may occur in non-embedded schemes. We demonstrate that the lack of conservation in partitioned schemes can lead to non-physical effects and propose conservative additive schemes based on partitioning the fluxes rather than the ordinary differential equations. A variety of SPERK schemes are presented, including an embedded pair suitable for the time evolution of fifth-order weighted non-oscillatory (WENO) spatial discretizations. Numerical experiments are provided to support the theory.

The study of reaction-diffusion processes is much more complicated on general curved surfaces than on standard Cartesian coordinate spaces. Here we show how to formulate and solve systems of reaction-diffusion equations on surfaces in an extremely simple way, using only the standard Cartesian form of differential operators, and a discrete unorganized point set to represent the surface. Our method decouples surface geometry from the underlying differential operators. As a consequence, it becomes possible to formulate and solve rather general reaction-diffusion equations on general surfaces without having to consider the complexities of differential geometry or sophisticated numerical analysis. To illustrate the generality of the method, computations for surface diffusion, pattern formation, excitable media, and bulk-surface coupling are provided for a variety of complex point cloud surfaces.

We apply the concept of effective order to strong stability preserving (SSP) explicit Runge-Kutta methods. Relative to classical Runge-Kutta methods, methods with an effective order of accuracy are designed to satisfy a relaxed set of order conditions, but yield higher order accuracy when composed with special starting and stopping methods. We show that this allows the construction of four-stage SSP methods with effective order four (such methods cannot have classical order four). However, we also prove that effective order five methods-like classical order five methods-require the use of non-positive weights and so cannot be SSP. By numerical optimization, we construct explicit SSP Runge-Kutta methods up to effective order four and establish the optimality of many of them. Numerical experiments demonstrate the validity of these methods in practice.

The Closest Point Method for solving partial differential equations (PDEs) posed on surfaces was recently introduced by Ruuth and Merriman [J. Comput. Phys. 2008] and successfully applied to a variety of surface PDEs. In this paper we study the theoretical foundations of this method. The main idea is that surface differentials of a surface function can be replaced with Cartesian differentials of its closest point extension, i.e., its composition with a closest point function. We introduce a general class of these closest point functions (a subset of differentiable retractions), show that these are exactly the functions necessary to satisfy the above idea, and give a geometric characterization this class. Finally, we construct some closest point functions and demonstrate their effectiveness numerically on surface PDEs

The Closest Point Method (CPM) is a method for numerically solving partial differential equations (PDEs) on arbitrary surfaces, independent of the existence of a surface parametrization. The CPM uses a closest point representation of the surface, to solve the unmodified Cartesian version of a surface PDE in a 3D volume embedding, using simple and well-understood techniques. In this paper we present the numerical solution of the wave equation and the incompressible Navier-Stokes equations on surfaces via the CPM, and we demonstrate surface appearance and shape variations in real-time using this method. To fully exploit the potential of the CPM, we present a novel GPU realization of the entire CPM pipeline. We propose a surface-embedding adaptive 3D spatial grid for efficient representation of the surface, and present a high-performance approach using CUDA for converting surfaces given by triangulations into this representation. For real-time performance, CUDA is also used for the numerical procedures of the CPM. For rendering the surface (and the PDE solution) directly from the closest point representation without the need to reconstruct a triangulated surface, we present a GPU ray-casting method that works on the adaptive 3D grid.

Eigenvalue problems are fundamental to mathematics and science. We present a simple algorithm for determining eigenvalues and eigenfunctions of the Laplace-Beltrami operator on rather general curved surfaces. Our algorithm, which is based on the Closest Point Method, relies on an embedding of the surface in a higher-dimensional space, where standard Cartesian finite difference and interpolation schemes can be easily applied. We show that there is a one-to-one correspondence between a problem defined in the embedding space and the original surface problem. For open surfaces, we present a simple way to impose Dirichlet and Neumann boundary conditions while maintaining second-order accuracy. Convergence studies and a series of examples demonstrate the effectiveness and generality of our approach.

We investigate the strong stability preserving (SSP) property of two-step Runge-Kutta (TSRK) methods. We prove that SSP TSRK methods belong to a simple subclass of TSRK methods, in which stages from the previous step are not used. We derive simple order conditions for this subclass. Whereas explicit SSP Runge-Kutta methods have order at most four, we prove that explicit SSP TSRK methods have order at most eight. We present TSRK methods of up to eighth order that were found by numerical search. These methods have larger SSP coefficients than any known methods of the same order of accuracy, and may be implemented in a form with relatively modest storage requirements. The usefulness of the TSRK methods is demonstrated through numerical examples, including integration of very high-order WENO discretizations

We study the linear stability of the fifth-order Weighted Essentially Non-Oscillatory spatial discretization (WENO5) combined with explicit time stepping applied to the one-dimensional advection equation. We show that it is not necessary for the stability domain of the time integrator to include a part of the imaginary axis. In particular, we show that the combination of WENO5 with either the forward Euler method or a two-stage, second-order Runge-Kutta method is linearly stable provided very small time step-sizes are taken. We also consider fifth-order multistep time discretizations whose stability domains do not include the imaginary axis. These are found to be linearly stable with moderate time steps when combined with WENO5. In particular, the fifth-order extrapolated BDF scheme gave superior results in practice to high-order Runge-Kutta methods whose stability domain includes the imaginary axis. Numerical tests are presented which confirm the analysis.

In this work we discuss a class of defect correction methods which is easily adapted to create parallel time integrators for multicore architectures and is ideally suited for developing methods which can be order adaptive in time. The method is based on integral deferred correction (IDC), which was itself motivated by spectral deferred correction by Dutt, Greengard, and Rokhlin [BIT, 40 (2000), pp. 241-266]. The method presented here is a revised formulation of explicit IDC, dubbed revisionist IDC (RIDC), which can achieve pth-order accuracy in "wall-clock time" comparable to a single forward Euler simulation on problems where the time to evaluate the right-hand side of a system of differential equations is greater than latency costs of interprocessor communication, such as in the case of the N-body problem. The key idea is to rewrite the defect correction framework so that, after initial start-up costs, each correction loop can be lagged behind the previous correction loop in a manner that facilitates running the predictor and M=p-1 correctors in parallel on an interval which has K steps, where K much greater than p. We prove that given an rth-order Runge-Kutta method in both the prediction and M correction loops of RIDC, then the method is order r×(M+1). The parallelization in RIDC uses a small number of cores (the number of processors used is limited by the order one wants to achieve). Using a four-core CPU, it is natural to think about fourth-order RIDC built with forward Euler, or eighth-order RIDC constructed with second-order Runge-Kutta. Numerical tests on an N-body simulation show that RIDC methods can be significantly faster than popular Runge-Kutta methods such as the classical fourth-order Runge-Kutta scheme. In a PDE setting, one can imagine coupling RIDC time integrators with parallel spatial evaluators, thereby increasing the parallelization. The ideas behind RIDC extend to implicit and semi-implicit IDC and have high potential in this area.

We propose a method to detect objects and patterns in textures on general surfaces. Our approach applies the Chan-Vese variational model for active contours without edges to the problem of segmentation of scalar surface data. This leads to gradient descent equations which are level set equations on surfaces. These equations are evolved using the Closest Point Method, which is a recent technique for solving partial differential equations (PDEs) on surfaces. The final algorithm has a particularly simple form: it merely alternates a time step of the usual Chan-Vese model in a small 3D neighborhood of the surface with an interpolation step. We remark that the method can treat very general surfaces since it uses a closest point function to represent the underlying surface. Various experimental results are presented, including segmentation on smooth surfaces, non-smooth surfaces, open surfaces, and general triangulated surfaces.

Many applications in the natural and applied sciences require the solutions of partial differential equations (PDEs) on surfaces or more general manifolds. The Closest Point Method is a simple and accurate embedding method for numerically approximating PDEs on rather general smooth surfaces. However, the original formulation is designed to use explicit time stepping. This may lead to a strict time-step restriction for some important PDEs such as those involving the Laplace-Beltrami operator or higher-order derivative operators. To achieve improved stability and efficiency, we introduce a new implicit Closest Point Method for surface PDEs. The method allows for large, stable time steps while retaining the principal benefits of the original method. In particular, it maintains the order of accuracy of the discretization of the underlying embedding PDE, it works on sharply defined bands without degrading the accuracy of the method, and it applies to general smooth surfaces. It also is very simple and may be applied to a rather general class of surface PDEs. Convergence studies for the in-surface heat equation and a fourth-order biharmonic problem are given to illustrate the accuracy of the method. We demonstrate the flexibility and generality of the method by treating flows involving diffusion, reaction-diffusion and fourth-order spatial derivatives on a variety of interesting surfaces including surfaces of mixed codimension.

Strong stability preserving (SSP) time discretizations were developed for use with spatial discretizations of partial differential equations that are strongly stable under forward Euler time integration. SSP methods preserve convex boundedness and contractivity properties satisfied by forward Euler, under a modified timestep restriction. We turn to implicit Runge-Kutta methods to alleviate this timestep restriction, and present implicit SSP Runge-Kutta methods which are optimal in the sense that they preserve convex boundedness properties under the largest timestep possible among all methods with a given number of stages and order of accuracy. We consider methods up to order six (the maximal order of SSP Runge-Kutta methods) and up to eleven stages. The numerically optimal methods found are all diagonally implicit, leading us to conjecture that optimal implicit SSP Runge-Kutta methods are diagonally implicit. These methods allow a significant increase in the SSP timestep limit, compared to explicit methods of the same order and number of stages. Numerical tests verify the order and the SSP property of the methods.

Level set methods have been used in a great number of applications in R² and R³ and it is natural to consider extending some of these methods to problems defined on surfaces embedded in R³ or higher dimensions. In this paper we consider the treatment of level set equations on surfaces via a recent technique for solving partial differential equations (PDEs) on surfaces, the Closest Point Method [?]. Our main modification is to introduce a Weighted Essentially Non-Oscillatory (WENO) interpolation step into the Closest Point Method. This, in combination with standard WENO for Hamilton-Jacobi equations, gives high-order results (up to fifth-order) on a variety of smooth test problems including passive transport, normal flow and redistancing. The algorithms we propose are straightforward modifications of standard codes, are carried out in the embedding space in a well-defined band around the surface and retain the robustness of the level set method with respect to the self-intersection of interfaces. Numerous examples are provided to illustrate the flexibility of the method with respect to geometry.

Diagonally split Runge-Kutta (DSRK) time discretization methods are a class of implicit time-stepping schemes which offer both high-order convergence and a form of nonlinear stability known as unconditional contractivity. This combination is not possible within the classes of Runge-Kutta or linear multistep methods and therefore appears promising for the strong stability preserving (SSP) time-stepping community which is generally concerned with computing oscillation-free numerical solutions of PDEs. Using a variety of numerical test problems, we show that although second- and third-order unconditionally contractive DSRK methods do preserve the strong stability property for all time step-sizes, they suffer from order reduction at large step-sizes. Indeed, for time-steps larger than those typically chosen for explicit methods, these DSRK methods behave like first-order implicit methods. This is unfortunate, because it is precisely to allow a large time-step that we choose to use implicit methods. These results suggest that unconditionally contractive DSRK methods are limited in usefulness as they are unable to compete with either the first-order backward Euler method for large step-sizes or with Crank-Nicolson or high-order explicit SSP Runge-Kutta methods for smaller step-sizes. We also present stage order conditions for DSRK methods and show that the observed order reduction is associated with the necessarily low stage order of the unconditionally contractive DSRK methods.

Keywords: diagonally split Runge-Kutta methods; Runge-Kutta methods; unconditional contractivity; strong stability preserving; time discretization; order reduction; stage order; hyperbolic PDEs