1. Introduction
Existing monotone cubic interpolation methods [1] [2] [3] [4] are successful in solving many practical problems. However, they are C1 and not fit for certain applications that require a higher degree of smoothness. For a C2 continuous monotone interpolation, a quintic polynomial is required [5] [6] [7] , or some subdivision of intervals needs to be performed [8] [9] [10] . An issue with these methods is that the derivative of the interpolation curve can have global oscillations and this limits their usage. There are methods existing for either monotone [11] [12] or with nonoscillation derivatives [13] [14] . In general, a monotone polynomial spline method requires certain constraints on their slope estimate to be satisfied.
In this article, we propose an unconditionally monotone interpolation method that is only quartic, with the C2 continuity over the entire domain. The new method has no nonphysical oscillation with its derivative and is 3rd order accurate in space.
2. The Solution Procedure
2.1. Description of the Problem
A set of points
are specified for
with
ordered increasingly as well as
, i.e.
and
hold for each integer i in range, as shown in Figure 1. Here N is the number of intervals defined by the
given points.
The problem is to construct an explicit polynomial curve
that passes all the given data points, i.e.,
, for
, and has a continuous second derivative. In addition, the curve must be strictly increasing (being monotonic
). Furthermore, the derivative of the interpolation function,
should have no unnecessary oscillations. Figure 2 demonstrates the oscillation in the derivative space of a monotone interpolation obtained with a simple Hermite-interpolation method applied to a radioactive particle energy distribution problem.
Figure 1. A monotone interpolation of a strictly increasing data set
.
Figure 2. Oscillation in the slope-space of a monotone interpolation. Although the slope of the spline function (red curve) is positive and matches the given areas in each interval defined by the black polylines, thus, exactly passes each data point, such a monotone interpolation is unfit in certain applications like rebinning a data set of radiation energy counts.
2.2. Reduction of the Problem
We consider a reduced problem in the slope space of the original problem. Let the slope function
, then if
satisfy that
(1)
(2)
and
has a continuous first derivative (being C1). Then, we can see that the integral of
(3)
not only passes through all the data points, but also has a well-defined second derivative. Therefore, we pursue the solution of the reduced problem by finding certain
that satisfy the above given constraints.
The domain
is divided to N intervals with interval i defined by
. The partial area
is bounded by two lines:
, and the boundaries of interval i, here
(4)
where
. A solution of the problem must satisfy Equations (1), (2), and (3), and without unnecessary oscillation.
The above description corresponds to certain statistics problems such as re-binning of a radioactive particle energy distribution (see Figure 3 on the next page), with the x-axis being the energy and y-axis being the denisty of particles in each energy bin. A requirement for energy re-binning is that the curve
has to be smooth for differentiation, i.e.
is C2. Not only that, the solution needs to have no oscillation for a minimal slope variation. Therefore, as the solution is mapped to a different bin-structure, it still makes physical sense.
There can be an infinite number of candidates of the solution. To limit the choices we consider a Hermite cubic-spline between an arbitrarily given pair of data points
and
(here L and R stand for the “left” and “right” boundaries of an interval) such that
(5)
where
is a non-dimensional parameter.
are the estimates of
,
are the estimates of
on the left and the right ends of a given interval.
are the Hermit cubic spline base functions that
(6)
Next, we will show that how
can be constructed with the above Hermite spline to satisfy Equations (1), (2), and (3), by properly choosing a set of control points.
Figure 3. A good monotone interpolation
should have no nonphysical oscillations in its slope space
besides matching the given areas under the black polylines bounded inside each interval exactly. The green curve is such a positive function
that satisfies the said constraints, obtained with the proposed interpolation method.
2.3. Area-Matching for Selecting Control Points on Interval-Walls
The slope of
,
at each inner data point can be estimated numerically. For example, using a quadratic interpolation on the three points
and
, one is able to obtain an estimate of
, except at the left and right boundaries. We consider the reduced problem as an interface reconstruction problem for volume conservation. The approach is to construct the geometry of the interface contained in interval i and to match the partial volume (area)
for each i. The interface piece constructed in interval i in general does not match with the pieces constructed in its neighbor intervals on boundaries. We will apply a Hermit spline later to eliminate the gaps and ensure a global slope continuation of
.
To start with, we consider a given internal
as shown in Figure 4. The average of
in this interval, hi is defined in Equation (4). We replace the subscripts “i” with “L” (left) and “i + 1” with “R” (right) for generality and let the width of the interval be
. We build a local Cartesian coordinate system
with its origin at
, and let a quadratic curve represent the interface, such that
For slope match at the left end
and the right end
, we have
The area matching means the integral of
over the interval is 0, or
Figure 4. A quadratic fitting in the local coordinate system
in
. The two end slopes and the area under
are exactly matched. The blue dots on the ends are to be repositioned by finding the intersection between the interval wall
and a Hermite spline passing the current yellow control points. The yellow control points in the middle are to be lifted/lowered later with another round of area-matching.
Solving for a, b and c, one arrives at
(7)
The constant term c carries the position of the interface at the middle of the interval to the 3rd order accuracy (for our fitting here is quadratic).
We will use these mid interval interface positions obtained from the above quadratic fitting as a set of control points to construct our first approximation of the solution. Each interval wall
intersects the Hermit spline curve passing the set of the mid-interval control points and the intersection is taken as a fixed control point. We have
control points
. The
of them with even subscripts are on the walls of the intervals and are fixed, the rest N of them above the middle points of intervals are to be shifted vertically by matching volumes again.
2.4. Area-Matching for Selecting Mid Interval Control Points
We are to construct a Hermit spline that passes all the control points. For an exact area match, we break each original interval into two subintervals about the control point in the middle of the interval, see Figure 5. With the two neighbor mid interval points, there are three mid interval control points involved in the area-matching. We compute the heights of the mid interval control points by solving a tri-linear linear system.
In Figure 6, the shadowed area
under a Hermite cubic spline in the interval
can be calculated with
(8)
Let the control point immediately left to
be
, and the
Figure 5. The yellow dots at the middle of each interval on their quadratic area fitting curves (green curves) are to be lifted or lowered. The blue control points are the intersection between the interval walls and a cubic spline (red curve) passing the yellow control points. The interval
is broken to two subintervals for another round of area-matching.
Figure 6. The area under a Hermit spline curve defined in a general interior interval
. To compute the area bounded for the original interval, two areas from the two subintervals are to be added together. In another word, each of the intervals
,
and
in this figure should be considered as a subinterval.
one immediately right to
be
. Evaluation of the above area integral provides that
(9)
where
(10)
Our choices of the x-slope terms are
and the y-slope terms are
(11)
They have been explicitly substituted in Equation (9) the expression of the physical area. The
terms depend only on the x-coordinates of the boundaries of an interval. The terms
are area integrals of the Hermit spline defined as
(12)
Each of
takes the values 0 or 1. These integrals are evaluated as the following
Now let us consider the ith interval, with which 3 control points are involved, see in Figure 7. They are
and
with
, and
to be determined. For area matching we
Figure 7. The area under the Hermit spline curves on interval i is the sum of two half-interval areas defined with Equation (13). Each area matching involves three mid interval control points (yellow). The solution of a tri-diagonal linear system determines their heights.
need to enforce that the sum of the area contribution from the left subinterval and the right subinterval to equal an known value that
(13)
Utilizing Equation (9) one arrives at
Because each of the
terms involves only the x-coordinates of the control points, it is a constant. The sum of areas under each Hermit cubit splines defined on a subinterval is simply a linear combination of
,
, and
. Therefore, we have a tri-diagonal linear system involving all intervals to solve in order to match the area exactly in each interval, with certain boundary conditions provided for the left-most and the right-most intervals.
2.5. Boundary Conditions
For the interval on the left boundary, we must specify the slope terms
. For the interval on the right boundary, we must provide the slope terms
as well. Currently we assume two kinds of boundary conditions. The first is a symmetrical boundary condition with which one sets a ghost control points at the reflection point of the nearest inner control point. The other one is a counter-symmetric condition by setting a ghost control point out of the boundary by extending the line-segment defined by the two nearest known control points involved (see Figure 8). These ghost points provide closure of the solution of the tri-diagonal linear system mentioned in the last subsection.
Figure 8. A demonstration of the symmetric (the right side) and the counter-symmetric (the left side) boundary conditions. A ghost control point (pink) is employed in each case.
2.6. Positivity in the Slope-Space
We have obtained
with a Hermite cubic spline on the control points computed in the previous sections. This solution may not necessarily be positive when
is close or equal to zero. Thus, we need to have a treatment in the possible case this negativity occurs. Fortunately, we have a nearly trivial treatment that is rather easy to perform, at least for isolated cases.
Considering the worst case that
, in the
space this means
must be a constant in the corresponding interval. Equivalently, the slope
must be zero everywhere in interval i, otherwise any variation would create some negative slopes then the monotone condition is violated. Specifically, the slopes at the two ends must be zeros.
Thus, we choose to enforce the slope terms
to zeros in the Hermite cubic spline in case an interval contains a point with a negative
. This means the two control points on interval boundaries are at the same height. Since we also assume a troubled interval is isolated, we lift the neighbor control
points at
and
to match the areas in the two neighbor
intervals. Because the area match condition is linear for a single variable
or
in either neighbor interval so the solution is trivial and does not affect rest of the intervals.
After the above treatment there will be no occurrence of
anymore. Therefore, the monotonicity of
is satisfied unconditionally, see Figure 9. Not to mention that no unnecessary oscillations are introduced with this local treatment.
Figure 9. With setting the end slopes to zeros (the green curve) and matching the area under the polylines again, the isolated slope-space negativity in a single interval (the red curve) is fixed. The locally modified spline still has a continuous derivative crossing the end points of the interval and does not affect the solution elsewhere.
3. An Unconditional Monotone C2 Spline
3.1. Constructing the Proposed Method
We have obtained a Hermite cubic spline in previous sections. The spline is non-negative, differentiable, and the area under it bounded by the walls of intervals and the x-axis exactly matches a specified area for each interval. Therefore, the integration of
provides a monotone interpolation of the data set
with an excellent quality. We have picked a general parametric form of
and
in Equation (5). In practice one can simply pick that
for interval “i”. All the previous discussions would still hold. This means we have an unconditionally monotone interpolation that is only quartic. It is an integral of the Hermit-cubic splines. Its order is lower than some of the existing monotone interpolation methods. Although we have split each interval to two, thus increased the number of knots. Nevertheless, since our analysis is done in the slope space for a cubic spline, the proposed method may be easier to handle than other monotone interpolation methods. The new monotone interpolation can be explicitly expressed as the follows
(14)
for the left subinterval of the original interval “i” with
.
The quartic functions
are integrations of the Hermite base functions
and
Let us define that
Then for the right subinterval
, we have
(15)
with
this time. The slope terms are defined as
if interval “i” contains no negativeness. In the case of a possible isolated negativeness the corresponding q terms are taken to zeros. Because the y (thus q) terms in the above equations are obtained with area matching, Equations (1), (2), and (3) are all satisfied. Because
is differentiable,
is an unconditionally monotone C2 quartic spline. Besides all the above, the proposed method does not have nonphysical oscillations on its derivative for we have explicit slope control with a Hermite cubic spline. This is a desirable feature. It is possible to further refine the solution by minimize the curvature (say) of
to adjust the level of the control points for minimal variation in slope. However, we are confident that the proposed solution is of third order accuracy with a quadratic interface reconstruction. Therefore, the control points are almost right on the ideal solution assuming which exists. Then, a further refinement is hardly necessary for a practical purpose.
3.2. In Two-Dimensions
Consider a set of data triplets
given for
and
with the properties
, and
,
for an arbitrary pair of integers i, j in the range. Can one construct a polynomial surface
that passes every data point (i.e.
for all feasible integer pair
), and has positive spatial derivatives
with no oscillation for both
and
?
If the given dataset is also consistent with a sufficient condition for data monotonicity
(16)
for all valid i and j, then, a method for interface-reconstruction in three-dimensions using volume conservation may be applied to constructing a two-dimensional monotone spline that is at least twice differentiable, with no oscillation on the first derivatives of
.
Without loss of generality, we assume
and modify a given data set satisfying Equation (16), by adding a layer of phantom data points by simply choosing
,
, and assigning
for all
,
. The monotone feature that
can be derived from Equation (16) using the boundary values.
The interface reconstruction is done in the
space. The “fluid volume” confined in each rectangular cell
, is defined by Equation (16), which is a difference expression for the cross-derivative at the cell center. Because slope estimates can be computed with certain difference schemes, there are enough data to support the definition of a local quadratic polynomial (or some other algebraic expression) that bounds
exactly in this cell.
The average of the intersections between the straight line
and the reconstructed local interface expressions in the surrounding cells can be used as a fixed control point. The control points above the middle point of each edge of a given cell can be set similarly (see Figure 10). A bicubic spline polynomial can be defined over all the control points and an exact volume match for all the cells can be performed to determine the heights of the control points at the center of cells. Solution of a sparse linear system is required because a volume integral under such a bicubic spline is a linear combination of the heights of neighbor control points above cell centers.
Finally, the monotone spline can be obtained with an integral of the bicubic spline
described above
Because
holds everywhere by construction and on the boundaries
and
the first derivatives of S are nonnegative, as the spatial integrals of
,
and
must hold everywhere inside the computational domain. Thus the spline is monotone. Because explicit slope control is provided for
with setting the control points by volume conservation, there is no oscillation on the derivatives of
for
has no oscillation. In addition, the spatial slopes is continuous for a bicubic polynomial
,
is then at least twice differentiable.
Figure 10. The top view of an internal cell (solid line) and its neighbors (dashed lines) are shown above. The blue control points are computed from local interface reconstruction and are fixed. The yellow control point at the center of each cell is going to be shifted in the direction perpendicular to the page to match volumes defined by Equation (16) by solving a linear system in order to construct a bicubic spline
.
4. An Application: Arbitrary Rebinning of Statistical Data
A set of statistical data pair
is given. Which can be considered as density of radioactive particle counts vs. energy, as shown in Figure 11. It is desired to bin the data in such a way that the number of counts is even in each bin. We are going to demonstrate how to solve this problem with constructing the proposed monotone interpolation.
4.1. Slope Estimate
At each internal knot
, we find its neighbor knots
and
and fit a quadratic polynomial
that passes the three data points and the solution is explicit that
For an estimate of slopes we can differentiate
and take
4.2. A Quadratic Polynomial Area-Match to Estimate Mid-Interval Control Points
For each interval
, we fit a quadratic polynomial
with
the slope estimates described above. One
Figure 11. The statistical counts of certain radio-active events are binned. The x-axis is energy density and the y-axis stands for the average event count per energy unit. This plot has the same interpretation as shown in Figure 2.
finds that
The value of
provides the estimate of the height of the mid interval control point
. The quadratic area fitting function
described
above has the 3rd order spatial accuracy and provides the initial heights of the mid interval control points.
4.3. A Hermit Spline to Determine Control Points above Knots
A Hermit cubic spline passing the set of mid interval area matching points obtained above intersects the line
and the intersection is taken as the position of control point
. In the unlikely case that this control point is below the x-axis, for positivity, we would lift it to
thus move this control-point closer to the true solution which is by definition above the x-axis. These control points are not to be moved.
4.4. Adjusting Mid-Interval Control Points
At this step, each original interval is broken to two subintervals of equal lengths. Each interval
corresponds to a mid-interval control point
.
is computed by solving a tri-diagonal linear system
defined by the area-matching Equation (13) for all i’s.
4.5. Constructing the Proposed Monotone Spline
All the control points are now determined. They define a unique Hermit spline curve
that passes all the control points, has a continuous derivative crossing the bin-walls, and matches the counts in each bin. We are able to integrate
to construct the monotone interpolation function
as described with Equation (14), Equation (15).
4.6. Solving for Evenly Dividing Counts in Each New Bins
Because
is monotone, one is able to evenly divide the vertical axis in the
space between
to M intervals. Let the new accumulated bin counts be
Then one solves for each j that
for the new bin-system. Because the derivative
is available at any point in the domain, a Newton- Raphson method
Figure 12. The statistical counts in Figure 11 are rebinned with the proposed monotone C2 non-oscillation quartic interpolation. The black curve is the slope function
. The red bin structure is the original and the green bin structure is for even bin counts in each bin.
will generate the solution quickly for
is monotone. The first guess can be picked with a search to find the ith bin that contains the solution, i.e.
then, simply taken to
.
A solution in the case of
and
is shown in Figure 12. One should easily understand the feature of no oscillation with
is necessary for a sensible solution of rebinning the statistical data.
5. Conclusion
A one dimensional monotone C2 quartic interpolation method is proposed. The original problem is reduced to a volume matching problem in the slope space to locate a set of control points on the solution curve. A quadratic fitting in each interval is employed to first approximately locate the control points with area matching. The solution of a tri-diagonal linear system is employed to relocate the control points above the middle of each interval to exactly matching the areas under a Hermit-cubic spline curve while fixing the control points above the original knots. The integration of the solution in the slope-space provides the desired unconditional monotone C2 interpolation. The proposed method has the feature of being of lower order (quartic) and with no undesired oscillation in the slope space. An application to the practical problem of rebinning a set of nonnegative statistical data shows the proposed method is effective in one-dimension. The proposed method may be extended to two-dimensions in a similar fashion with a higher smoothness for data sets that satisfy an extra constraint.
Acknowledgements
This work is performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
The authors would like to express appreciation to Simon Labov and Chad Noble of Lawrence Livermore National Laboratory for their encouragement and support.