Automatic Generation of Interpolants for Lattice Samplings: Part I — Theory and Analysis

Given a set of discrete data points that reside on a lattice, one of the key fundamental operations that one can perform on those data is interpolation. This appears in many different contexts in science and engineering: one may use interpolation to fill in data between points in time-series data, to fill in missing data between pixels in an image, or missing data between voxels in a volumetric image. In fact, many scientific visualization algorithms require a continuous representation of a signal as an input.

As simple as it may seem, interpolation of lattice sampled data can quickly become complicated. To many, interpolation is synonymous with linear interpolation, and if we wish to interpolate in higher dimensions, we can simply linearly interpolate in each dimension separately (i.e. a tensor product extension). However, there is much more freedom available that we may take advantage of when we move past one dimension. For example, there are many alternative interpolants (or basis functions) we may choose instead of a simple linear interpolant, many with higher accuracy and/or higher smoothness than a simple linear interpolant. Not only this, but we may have more freedom in terms of where samples are placed when the dimension of the space is $\geq 2$. We say ‘may’ because for some problems, one is forced to use a specific lattice structure - typically the Cartesian lattice, whereas other problems allow lattice choice as a degree of freedom.

In this work, we build a framework for interpolation on lattices, and show how to do it fast. Our work originated in scientific visualization, in particular working with volumetric data, but the observations we make are generally applicable to other domains as well as higher dimensions. We will mainly keep to the 3-dimensional case, but we will note how to generalize to higher dimensions when it is necessary.

Returning to the context of scientific visualization, interpolation speed is important in practice because it facilitates interactivity which, in turn, allows for fast iteration over data-sets — the ability to quickly (visually) explore data is one of the key strengths of visualization. But more often than not, in practice, interpolation boils down to “use trilinear interpolation”, and it is clear why: trilinear interpolation is fast, it is implemented in hardware for all modern graphics processing units; it is straightforward to use, it inlvolves a single texture fetch instruction for modern GPUs; and it looks “good enough” in practice. However, linear interpolation has a number of shortcomings. Perhaps the most obvious is smoothness. Trilinear interpolation is based on a piece-wise linear univariate function extended via tensor product along all three dimensions; this ensures a continuous approximation but introduces discontinuities in the first and higher order derivatives. While this could be mitigated by using filtering methods, there is additional memory overhead associated with such methods (Alim et al., 2010). Linear interpolation also introduces visual artifacts when interpolating at lower resolutions. This is again related to the smoothness of the trilinear interpolant, but also has deeper roots in sampling and approximation theory. In short, the geometry and support of the chosen interpolant dictate how a reconstruction attenuates high-frequency content in the Fourier domain. This gives rise to different types of artifacts in data reconstruction (Csébfalvi and Rácz, 2016). Thus, there is value in exploring the design space of interpolants with the goal of minimizing error; in practice, one also needs to account for the computational costs associated with different interpolants.

Before moving further, we need to clarify what we mean by non-Cartesian computing. We use the term “non-Cartesian computing” to indicate a move away from tensor-product extensions. This can take the form of non-Cartesian lattices — take Figures 1 and 2 as examples that show the typical cubic Cartesian (CC) and square lattices in three and two dimensions, as well as more efficient non-Cartesian variants, namely the Body Centered Cubic (BCC) lattice, and the Face Centered Cubic (FCC) lattice, both of which have received a good amount of attention in the field of scientific visualization. We may also use the term “non-Cartesian computing” to denote the move away from tensor-product interpolants to non-separable interpolants. In the Cartesian world, one typically uses tensor product B-splines as basis functions; however, non tensor-product splines exist, and have been investigated in the literature (Csébfalvi and Rácz, 2016). One example is the ZP-element which requires fewer memory accesses per reconstruction than the cubic tensor-product B-spline but has similar approximation properties and provides smoother, more isotropic reconstructions (Entezari and Möller, 2006). When comparing different sampling lattices, one may argue that a true comparison for these lattices would be their Voronoi splines (Mirzargar and Entezari, 2010), which, until this work, have yet to have a dedicated GPU implementation.

What do we gain by adding additional complexity to our interpolation and data processing schemes? The seemingly benign difference in lattice and basis function changes the properties of the resulting approximation scheme. In fact, in 2D, an optimal hexagonal sampling scheme (Fig. 2) yields $14\%$ better approximations when combined with an appropriate interpolant, whereas on the BCC lattice (Fig. 1)— the optimal sampling lattice in 3D — the improvement is around $30\%$ (Entezari, 2006). Despite these advantages, the benefits of non-Cartesian computing have remained elusive. This is largely due to the complex algebro-geometric interplay between the lattices and the basis function they are paired with. These basis functions are typically composed of box splines which have intricate piece-wise polynomial structures which are challenging to evaluate in an efficient manner. Most treatments have dealt with evaluation and implementation issues in a somewhat adhoc manner.

Our goal in this work is to make the advantage of non-Cartesian computing more tangible — explicitly, we want to make the use of alternative interpolants/lattices more practical and usable. Our main contribution is a holistic analysis of approximation spaces with a careful focus on fast implementation. We outline a framework for translating interpolation bases into fast interpolation schemes for piece-wise polynomial interpolants. We then show how to use this pipeline to create unified fast implementations of many of the interpolants in the non-Cartesian volumetric visualization literature. It is important to mention that we obtain a form of evaluation that is agnostic to platform. We could generate CPU code, GPU code, or Verilog from our intermediate representation. Our GPU implementation takes advantage of the in-built tensor product linear filtering hardware (i.e. linear texture fetches) on contemporary GPUs to reduce the amount of memory accesses needed for reconstruction, however this does not always lead to optimal results.

We provide the theoretical analysis in this paper, and in the subsequent work we provide extensive implementation details. In our supplementary material we share the worksheet we use to automate the analysis as well as the code generation module, and pre-generated CUDA PTX code. To summarise, in this part, our contributions are as follows:

• •

  We provide a unified framework of analyses for combinations of splines and lattices from the context of scientific visualization.

• •

  We show how to use trilinear interpolation to reduce memory fetches and increase performance of our implementations.

• •

  We provide implementations for many interpolants in the literature that have not yet had robust GPU implementations.

• •

  We show the relative performance between GPU implementations, and show generality by moving towards a 4 dimensional example.

The remainder of this paper is organized as follows. In Section 2, we review the literature in our native context, scientific visualisation, but we will also pay note to some of the work done in one dimension and in image processing. In Section 3, we provide the background necessary for our framework. In Section 4, we attack the problem at its root, the convolution sum. We manipulate the convolution sum until it emits a form suitable for fast evaluation — this can be thought of as unrolling the convolution sum loop. In Section 5 we generate code for various interpolants (box splines and Voronoi splines) some of which have not yet appeared in the literature. We then evaluate their performance on the CPU and GPU with respect to the application of volume rendering.

Interpolation is a fundamental operation in scientific visualization. Many visualization tasks (volumetric rendering, iso-surface extraction, flow visualization etc.) rely on the interpolation of a continuous function from a discrete set of values. In one dimension, the family of B-splines are the prototypical interpolant — the linear B-spline is the least computationally expensive in this family, and has received the most attention in practice (Unser, 1999). However, it is well understood what B-splines are members of the maximal order minimal support (MOMS) family, which are the optimal interpolants for a given support (Blu et al., 2001). This is not the only option though, there are other constraints one may wish to design around. For example: arithmetical complexity; total number of memory accesses (i.e. support) needed for a single point reconstruction; and smoothness are all parameters that can be tweaked. If one requires infinite differentiability, then CINAPACT splines are an appropriate interpolant candidate (Akram et al., 2015).

The design space for uni-variate interpolants is well understood. Appropriate choice of interpolant is not so clear in higher dimensions. For the domain of scientific visualization, the move to three dimensions allows practitioners the freedom to consider non-seperable splines. In this case, it is clear that the tensor product B-splines are no longer MOMS splines; the generalized ZP-element shows this, since it has lower support than (but the same approximation order as) the cubic B-spline (Entezari and Möller, 2006).

There are multiple works that investigate non-separable basis functions on non-Cartesian lattices; box splines and weighted linear combinations of box splines often appear as popular candidates (de Boor et al., 1993; Domonkos and Csébfalvi, 2010; Csebfalvi, 2013). However, it is fairly well known that evaluating a box spline numerically is quite difficult; the recursive form for box splines is unstable if naively implemented (de Boor, 1993). One can rectify this by ensuring that any decision made while evaluating a point on the spline’s separating hyperplanes is consistent (Kobbelt, 1997). Even then, recursive box spline evaluation is unsuitable for use on the GPU — the conditional nature and large branching behaviour of the recursive form will lead to branch divergence and stall execution units on the GPU, leaving its resources underutilized. It is more convenient to work with either the Bernstein Bézier (BB) form, or the explicit piece-wise polynomial (PP) form of a box spline which has been characterized in closed form (Horacsek and Alim, 2016). Mathematically these two forms are equivalent, the BB form is simply one specific factorization of the PP form that lends itself to evaluation with De Castlejau’s algorithm. Generally, a polynomial can be factorized in many different ways.

The BB form has been used in fast GPU evaluation schemes for box splines on non-Cartesian lattices. In the work of Kim et al. (Kim, 2017), the symmetry of a box spline is used to create look-up tables of the BB polynomial coefficients. Since these coefficients are rational, they are stored as pairs of integers and any division occurs at run time. This is particularly convenient on the GPU, as it allows one to write one function for De Castlejau’s algorithm, and each separate region of the spline to be evaluated is a set of integers that can be looked up based on an analysis of the regions of the spline’s mesh at runtime. However, this is somewhat wasteful since most regions within a spline’s mesh are related by a rotation and/or a shift (i.e. an affine change of variables). As we will see in Section 4, this change of variables allows us to re-write the polynomial within the region of evaluation as a transformation followed by a polynomial evaluation of some chosen reference polynomial (provided the spline has an appropriate structure). This reference polynomial can be evaluated in BB form, however, we choose to use a Horner factorization of the PP form — this allows us to reduce the amount of operations needed to calculate a given polynomial (Ceberio and Kreinovich, 2004). This is a generalization of the approach used in the work of Finkbeiner et al. (Finkbeiner et al., 2010).

On the BCC lattice multiple box splines have been investigated for volumetric visualization. The linear and quintic box splines were among the first used by Entezari et al. (Entezari et al., 2008, 2009) for volumetric data visualization. These splines are particularly interesting as they mimic the geometry of the BCC lattice and they attain the same order of approximation as the trilinear and tricubic splines on the CC lattice, but with smaller support. Then there is the quartic box spline of Kim et al. which has even smaller support than the quintic box spline, but the same order of approximation (Kim, 2013c). Additionally, in that work, Kim et al. proposed a 12-direction box spline with tensor product structure; this box spline has a large support size, but is reasonably fast compared to the other proposed box splines due to its tensor product flavour — in our implementation, we call the 12-direction box spline the tricubic box spline on the BCC lattice.

The FCC lattice has not received as much attention from researchers as the BCC lattice; we are only aware of few works that investigate box splines on the FCC lattice. In particular, Kim et al. have investigated a six direction box spline that respects the geometry of the FCC lattice (Kim et al., 2008a; Kim, 2013a). However, it is also true that the 12-direction box-spline (Kim, 2013c) is usable on the FCC lattice. The proposed splines of Csebfalvi et al. may also be generalized to the FCC lattice (Csebfalvi, 2013; Domonkos and Csébfalvi, 2010).

We are only aware few works that attempt the use of a non tensor-product box spline on the CC lattice for visualization. The work of Entezari et al. generalizes the Zwart-Powell (ZP) element of two dimensions to three dimensions. This also maintains the same order of approximation as the cubic tensor product spline, but with smaller support (Entezari and Möller, 2006). Even though this spline tends toward higher fidelity approximations, we are not currently aware of any GPU implementation of this spline, and we provide one in the supplementary material.

There is also the work of Csebfalvi et al., where a shifted and re-weighted box spline is designed so as to map easily to linear interpolation among the cosets of the BCC lattice (Csebfalvi, 2013). This method has the advantage of being extremely easy to implement on the GPU, has respectable reconstruction quality and runs at decent speeds compared to trilinear interpolation. There has also been some work investigating the use of splines designed for one lattice on another (Csébfalvi and Rácz, 2016). The use of direction vectors that do not correspond to principle lattice directions helps distribute frequency content more evenly in the Fourier domain.

Another intuitive idea is to look at the Voronoi cell of a lattice and convolve that with itself to obtain an interpolant, Figure 1 shows the CC, BCC and FCC lattices with their Voronoi cells. This produces a valid approximation scheme, but is quite expensive to compute (Mirzargar and Entezari, 2010; Van De Ville et al., 2004; Mirzargar and Entezari, 2011). Until this work, there was no GPU implementation available for these splines. These fit within our pipeline, and we provide an implementation for the BCC and FCC lattices — Voronoi splines on the CC lattice correspond to the tensor product B-spline.

Finally, while not strictly related to spline evaluation, it is important to discuss quasi-interpolation, since an approximation space may not harness the full approximation power of a given basis function without proper pre-filtering. Most of the basis functions discussed so far are not interpolating. If a basis is stable (i.e it forms a Riesz basis) then it is possible to process the input data so that the resulting reconstruction interpolates data (Alim, 2012). However, this does not always yield the best reconstruction. Moreover, some bases cannot be made interpolating, so what is to be done in those cases? Quasi-interpolants ensure that a reconstruction harnesses the full approximation order of a space — in general it is a good idea to prefilter data with a quasi-interpolant if a basis is not interpolating. This is related to an error kernel analysis, where the data convolved with a filter is guaranteed to decay with the order of the basis (Blu and Unser, 1999; Alim, 2012). While this is an important ingredient to a good reconstruction scheme, in this work we are not concerned with the approximation properties of a basis function, only on the speed at which it can be evaluated.

When we work with data that are uniformly sampled in one dimension there is only one possible sampling scheme — take equally spaced samples along the indepenent axis. The straightforward extension of this to multiple dimensions is to sample evenly in all cardinal directions, i.e. a Cartesian sampling. However, the situation is generally more complex in higher dimensions.

To move away from Cartesian lattices and account for this complexity, a function can be sampled according to an invertible matrix. This matrix — known as the generating matrix — represents a lattice structure. That is, if we have a real valued function $f(\mathbf{x})$ in some dimension $s$, we define the sampled version of that function with respect to a given lattice as $f_{\mathbf{n}}:=f(hL\mathbf{n})$. Here, $h$ is a real valued scale parameter that controls the granularity or coarseness of the sampling pattern (i.e. the sampling rate), $\mathbf{n}\in\mathbb{Z}^{s}$ and $L$ is the generating matrix of the lattice. In any dimension, the $s\times s$ identity matrix generates a Cartesian sampling. In two dimensions the hexagonal and quincunx lattices are generated by the matrices

respectively (Figure 2). In three dimensions, the BCC and FCC lattices are generated by

respectively (Figure 1). It is worth noting that a lattice may be generated by any number of different generating matrices — one may construct an equivalent lattice by carefully choosing vectors in $L\mathbb{Z}^{s}$. In general, finding a minimal basis for a given generating matrix is an NP-hard problem (Conway and Sloane, 2013), as such it is difficult to find a relationship between matrices that generate the same lattice.

While lattice theory is quite deep and general, we only adhere to the case where $L\in\mathbb{Z}^{s\times s}$ and $L$ has full rank. This is because, in practice, we use lattices as array-like data structures and it is not possible to index via real valued numbers. As a consequence, it is not possible for us to directly represent functions on the hexagonal lattice within our framework. However, it is possible to stretch an integer lattice into another non-integer lattice; for example, in Figure 2, the hexagonal lattice can be seen as a distorted quincunx lattice.

Another fundamental object related to a lattice is its Voronoi region, which is the set of points closest to the lattice site situated at the origin (both Figure 1 and Figure 2 show each lattice’s Voronoi region). The Voronoi region of a lattice embodies all the properties of a lattice — from it one may derive the shortest vector of a lattice (i.e. the shortest vector problem which is NP-complete) as well as set of shortest spanning vectors, the symmetry group of the lattice, packing efficiency etc.  (Conway and Sloane, 2013; Viterbo and Biglieri, 1996). Unsurprisingly, computing the Voronoi region of a lattice is NP-complete as well, but for reasonably specified bases in lower dimensions, the problem is tractable via the Diamond Cutting algorithm (Viterbo and Biglieri, 1996).

To reconstruct a function on a given lattice, we choose a basis function $\varphi(\mathbf{x})$, shift it to each lattice site, and modulate it by the coefficient stored at that lattice site. The space of all functions spanned by the lattice shifts of $\varphi$ is defined as the set

where $\mathbf{c}$ is a real-valued coefficient sequence associated with the lattice. Obtaining a good approximation $\tilde{f}(\mathbf{x})\in\mathbb{V}(\varphi,L,h)$ to $f(\mathbf{x})$ equates to finding an appropriate set $\{c_{\mathbf{n}}\}$; this is highly dependent on the choice of $\varphi(\mathbf{x})$ and $L$. If we are so lucky that $\varphi(\mathbf{x})$ is interpolating on $L$ (i.e. $\varphi(\mathbf{0})=1$ and $\varphi(L\mathbf{n})=0$ when $\mathbf{n}\in\mathbb{Z}^{s}/\{\mathbf{0}\}$), then we may choose $c_{\mathbf{n}}=f(hL\mathbf{n})$. However, in general, some preprocessing must be performed in order to achieve optimal error decay between the original function and its approximation. This is beyond the scope of this article, but the interested reader may refer to some of the related work (Unser, 2000; Blu and Unser, 1999; Alim, 2012). Since any asymptotic analysis is outside of the scope of this paper, without loss of generalization we set $h=1$.

As for specific choices of $\varphi$, box splines have been investigated extensively in scientific visualization. Box splines have many equivalent definitions, but perhaps the most intuitive is the convolutional definition where one successively convolves an indicator function along the $\boldsymbol{\xi}_{i}$ direction vectors (de Boor et al., 1993). We collect $n$ of these vectors in an $s\times n$ matrix

In our case, we are interested only in matrices where $n\geq s$ and the dimension of the range of $\Xi$ is equal to $s$; we require these for a valid approximation scheme (de Boor et al., 1993). With this, a box spline can be defined recursively as

where ${\Xi/\{\boldsymbol{\xi}\}}$ is interpreted as removing a direction vector $\boldsymbol{\xi}$ from the matrix $\Xi$. When $n=s$, we define $M_{\Xi}(\mathbf{x})$ as

which provides a base case for the recursion.

The box spline $M_{\Xi}(\mathbf{x})$ is compactly supported within a polytope defined by the Minkowski sum of the direction vectors in $\Xi$. We denote the support as $\text{supp}\,M_{\Xi}$; it is given by

$M_{\Xi}$ is piecewise polynomial within its support. The polynomial pieces are delineated by all $(s-1)$-dimensional hyperplanes spanned by the direction vectors in $\Xi$. Denoting $\mathbb{H}$ as the set of all such hyperplanes, we can define a fine mesh that is formed by the lattice shifts of $M_{\Xi}$:

Each region of the mesh is a convex polytope that contains a polynomial, typically these are different polynomials, however, it is possible that some regions may share the same polynomial. In general, it is possible to derive the explicit piece-wise polynomial form for each of these regions (Horacsek and Alim, 2016)— we will use this explicit form in Section 4 to determine a simpler representation of the convolution sum in Equation 3.

Our framework also accounts for Voronoi splines, which can be intuitively defined as successive convolutions of a lattice’s indicator region with itself. More precisely, Voronoi splines are defined as

where $\ast$ denotes convolution, and $\chi_{L}(\mathbf{x})=1/|\det{L}|$ if $\mathbf{x}$ is in the Voronoi region of the lattice $L$, and $0$ otherwise. Voronoi splines are also piece-wise polynomial, and their piece-wise polynomial form may be obtained by first redefining $V_{0}(\mathbf{x})$ as the sum of a collection of constant valued box splines; higher order $V_{n}(\mathbf{x})$ may be obtained by the convolution definition (Mirzargar and Entezari, 2010). Unrolling the recursion, we obtain $V_{n}(\mathbf{x})=V_{0}(\mathbf{x})^{\ast n}$ where $\ast n$ denotes successive convolutions. If $V_{0}(\mathbf{x})$ is the sum of box-splines, then one may use the multinomial theorem to expand $V_{0}(\mathbf{x})^{\ast n}$ into the sum of many box-spline convolutions; convolving a box-spline with another box-spline simply produces a (possibly shifted) box-spline. Thus, by summing all these box-splines, we arrive at the final piece-wise polynomial form. However, these calculations are not trivial. For the BCC lattice, $V_{0}(\mathbf{x})$ is the sum of 16 box-splines, $V_{1}(\mathbf{x})$ (analogous to the tri-linear B-spline on the CC lattice) is the sum of 136 box-splines, $V_{2}(\mathbf{x})$ is the sum of 816 and $V_{3}(\mathbf{x})$, the sum of 3876 (Mirzargar and Entezari, 2010). While these numbers are not intractably large, keep in mind that one must compute a box-spline’s polynomial representation at every iteration, which can take a non-trivial amount of time; on the order of a few dozen minutes, for the splines that sum to $V_{3}(\mathbf{x})$. Spread across 4 CPU cores, this took us about a month to compute.

The convolution sum in Equation 3 is the focus of this section. There are two important parts to evaluating a function of this form; the first is the representation of the coefficient set $c_{\mathbf{n}}$ (i.e. the memory layout), the second is the evaluation of the basis function $\varphi$. Concerning $c_{\mathbf{n}}$, we choose a specific representation in which each $c_{\mathbf{n}}$ corresponds to a fetch from an $s$-dimensional array. Thus, for any lattice, we choose to treat it as a union of shifted Cartesian lattices (also known as a Cartesian coset decomposition). This allows us to use texture memory on GPUs to store each coset, as shown in Figure 3. Provided we combine all reads from a given coset, this also allows us to take advantage of data locality (i.e. cache). It also allows us to use trilinear fetches as building blocks for our interpolation schemes — that is, we attempt to write our interpolation schemes as weighted sums of trilinear texture fetches on the different cosets of the lattice. This does not necessarily reduce computation per evaluation, but rather reduces the number of texture fetches needed to evaluate the convolution sum, and therefore reduces bandwidth.

With the memory layout solidified, we turn to the basis function — we apply a series of manipulations to “unroll” the convolution sum. We build some insight into how to transform $\varphi$ into a more convenient form for evaluation with some running examples. We consider the linear tensor product spline and the Zwart-Powell (ZP) element on the two-dimensional Cartesian lattice. These have direction matrices

respectively. On the quincunx lattice, we use a tensor product style spline defined by

These box splines are shown in Figure 4. While this is nowhere near an exhaustive list of interpolants, the examples we provide are simple enough to understand in 2D and show the intricacies of the procedure.

In this work, we validate our methodology by demonstrating that our framework works for many different test cases. While performance is important, and we do report limited performance metrics in this work, we defer in-depth performance analysis to Part II (Horacsek and Alim, 2021). To this end, we choose two simple applications that demonstrate the generality of our framework. The first is volumetric rendering, in which we generate implementations for various different spline lattice combinations, many of which have not had fast GPU implementations. The second application is simply function approximation, however, we do this on the $\mathcal{D}_{4}$ lattice on the CPU. To our knowledge, there relatively little works that investigate function approximation in this space. We are not aware of any that assess the convergence of such spaces, nor any that provide implementations.

We do not “implement” Algorithm 1 directly; we use it simply as a template for code generation. That is, the input to our pipeline is a list of (polyhedral region, polynomial) pairs, and the output is generated LLVM code (Lattner and Adve, 2004). For our volumetric rendering experiments, we use the LLVM code to generate PTX code, for our other tests we emit x86_64 assembly. We have constrained ourselves to the CUDA ecosystem as it is the most popular GPGPU platform, however, LLVM contains both a backend for AMD GPUs and a platform agnostic SPIR-V backend that is capable of targing OpenCL + OpenGL devices (Kessenich et al., 2018). The details of the code generation stage, and a more detailed breakdown of performance is detailed in part II (Horacsek and Alim, 2021).

The processing time for a given interpolant varies based on the spline, higher order splines with large support take much longer, but the process takes typically on the order of hours on an Intel i7 7700 with 16GB of DDR4 RAM at 2333 MHz. However, the processing time of the pipeline is a one time computation, and need not be repeated once an effecient form has been derived. We run our tests on the generated code.

We presented a generalized framework for analysing multidimensional splines on non-Cartesian (and Cartesian grids), with the target of producing fast evaluation schemes for said spline spaces. While this is the main contribution of the work, we also produced performant code for the notoriously complex Voronoi splines on the FCC and BCC lattice which have not yet had efficient implementations. We also demonstrated the computational behaviour of our approach as spline size increased, showing reasonable computational increase as the support of a spline increases—however we saw that it is difficult to predict performance based on support size alone. Finally, we investigated the performance in 4-dimensions, and provided an imlpementation (and quasi-interpolant filter) for a spline on the $\mathcal{D}_{4}$ lattice. The entire pipeline of our worksheet is implemented within a SageMath worksheet, and is available on github (Developers, 2016; Horacsek, 2021). Further details related to the code generation step of our pipeline, as well as detailed performance results are presented in part II of this work (Horacsek and Alim, 2021). Future theoretical work is focused on extending this framework to more classes of splines (i.e. the exponential box splines), designing optimized interpolants and assessing the quality of interpolants in a systematic manner.