插值算法的Overshoot

高阶插值很多都有overshoot,我们期望插值kernel是单调的,这篇博客介绍几种无overshoot的插值算法

一、附加约束的Cubic Spline

https://github.com/conveyal/gridualizer/pull/24/commits/0953f195e17cf2890b6a883e082bd6ccf0217b2e

copy的js代码

Plain text
Copy to clipboard
Open code in new window
EnlighterJS 3 Syntax Highlighter
import {getGridValue} from '../util'
/**
* A 2D constrained spline interpolator.
*
* Polynomial (e.g. bicubic) interpolation is prone to oscillation and overshoot.
* Splines are better behaved but can still overshoot. We can sacrifice smoothness and add a no-overshoot constraint. * See: http://www.korf.co.uk/spline.pdf
* Kruger, CJC. Constrained Cubic Spline Interpolation for Chemical Engineering Applications.
*
* A pre-calculated bicubic interpolation patch.
* For a 4x4 grid of samples, this allows us to calculate interpolated values between the central four samples.
* By pre-fitting the curves in one dimension (y) and proceeding with the interpolation row by row,
* we re-use most of the computation from one output pixel to the next.
*/
const SplineInterpolator = (grid, gridX, gridY) => {
// First, find grid coordinates for the sixteen cells.
// Produces a bicubic interpolation patch for a one-cell square.
// The patch extends one cell east and south of the specified grid position, but uses 16 cells in the grid.
// Deal with the edges of the input grid by duplicating adjacent values.
// It's tempting to do this with typed arrays and slice(), but we need special handling for the grid edges.
const x0 = (gridX === 0) ? gridX : gridX - 1 // Handle left edge
const x1 = gridX
const x2 = (gridX + 1 >= grid.width) ? gridX : gridX + 1 // Handle right edge
const x3 = (gridX + 2 >= grid.width) ? gridX : gridX + 2 // Handle right edge
const y0 = (gridY === 0) ? gridY : gridY - 1 // Handle top edge
const y1 = gridY
const y2 = (gridY + 1 >= grid.height) ? gridY : gridY + 1 // Handle bottom edge
const y3 = (gridY + 2 >= grid.height) ? gridY : gridY + 2 // Handle bottom edge
const p00 = getGridValue(grid, x0, y0)
const p01 = getGridValue(grid, x0, y1)
const p02 = getGridValue(grid, x0, y2)
const p03 = getGridValue(grid, x0, y3)
const p10 = getGridValue(grid, x1, y0)
const p11 = getGridValue(grid, x1, y1)
const p12 = getGridValue(grid, x1, y2)
const p13 = getGridValue(grid, x1, y3)
const p20 = getGridValue(grid, x2, y0)
const p21 = getGridValue(grid, x2, y1)
const p22 = getGridValue(grid, x2, y2)
const p23 = getGridValue(grid, x2, y3)
const p30 = getGridValue(grid, x3, y0)
const p31 = getGridValue(grid, x3, y1)
const p32 = getGridValue(grid, x3, y2)
const p33 = getGridValue(grid, x3, y3)
// Create interpolations through each of the four columns
const columnInterpolator0 = ConstraintedSplineInterpolator(p00, p01, p02, p03)
const columnInterpolator1 = ConstraintedSplineInterpolator(p10, p11, p12, p13)
const columnInterpolator2 = ConstraintedSplineInterpolator(p20, p21, p22, p23)
const columnInterpolator3 = ConstraintedSplineInterpolator(p30, p31, p32, p33)
return function (yFraction) {
// Perform curve fitting in the second (x) dimension based on the pre-fit curves in the y dimension.
const p0 = columnInterpolator0(yFraction)
const p1 = columnInterpolator1(yFraction)
const p2 = columnInterpolator2(yFraction)
const p3 = columnInterpolator3(yFraction)
// Return the one-dimensional interpolator for this row.
return ConstraintedSplineInterpolator(p0, p1, p2, p3)
}
}
SplineInterpolator.tileOffset = 0.5
export default SplineInterpolator
/**
* Given four adjacent values a, b, c, d, fit a constrained cubic spline to them,
* sacrificing smoothness to prevent overshoot.
* The returned function provides interpolated values between b and c using a and d to
* determine the slope going into and out of the b-c interval.
* The original paper handles the general case where data points are (x,y) pairs.
* In our case, the four points are always evenly spaced, so we assign X coordinates of
* -1, 0, 1, 2 knowing we will perform interpolation between the second and third points.
* This greatly simplifies the equations, because it gives many differences, multipliers,
* and denominators have a value of 1.
*/
const ConstraintedSplineInterpolator = (a, b, c, d) => {
const test = function (a, b, c, d) {
let out = []
let interpolator = new ConstraintedSplineInterpolator (a, b, c, d)
for (let i = 0; i < 100; i++) {
out[i] = interpolator(i / 100)
}
console.log(out.join(","))
}
// Optimization: if b and c are equal, interpolate a straight line
if (b === c) return (x) => b
const bSlope = slope(a, b, c) // equation 7a
const cSlope = slope(b, c, d) // equation 7a
const bSlope2 = (2 * cSlope + 4 * bSlope) + (6 * (c - b)) // equation 8
const cSlope2 = (4 * cSlope + 2 * bSlope) - (6 * (c - b)) // equation 9
const kd = (cSlope2 - bSlope2) / 6 // equation 10
const kc = bSlope2 / 2 // equation 11
const kb = (c - b) - kc - kd // equation 12
const ka = b // equation 13
// return (x) => b + (x * (c - b)) // TEST - should be exactly the same as bilinear
// The returned function takes an x value in [0, 1] expressing the position between b and c,
// and returns the interpolated value.
return (x) => ka + kb * x + kc * x ** 2 + kd * x ** 3
}
/**
* Find the target slope at a particular data point.
* y here is the same as y sub i in the equations in the paper.
* yPrev here is y sub i-1 and yNext here is y sub i+1.
* The equations are simplified significantly by the fact that we know all of our points are exactly one unit apart.
*/
const slope = function (yPrev, y, yNext) {
const prevSlope = y - yPrev
const postSlope = yNext - y
if (prevSlope === 0 || postSlope === 0) return 0 // necessary condition for no overshoot
if (Math.sign(prevSlope) !== Math.sign(postSlope)) return 0 // includes case where only one slope is zero
return 2 / (1 / postSlope + 1 / prevSlope)
}
// Finding the slope at the endpoints is messy because it requries recursively computing slopes.
import {getGridValue} from '../util' /** * A 2D constrained spline interpolator. * * Polynomial (e.g. bicubic) interpolation is prone to oscillation and overshoot. * Splines are better behaved but can still overshoot. We can sacrifice smoothness and add a no-overshoot constraint. * See: http://www.korf.co.uk/spline.pdf * Kruger, CJC. Constrained Cubic Spline Interpolation for Chemical Engineering Applications. * * A pre-calculated bicubic interpolation patch. * For a 4x4 grid of samples, this allows us to calculate interpolated values between the central four samples. * By pre-fitting the curves in one dimension (y) and proceeding with the interpolation row by row, * we re-use most of the computation from one output pixel to the next. */ const SplineInterpolator = (grid, gridX, gridY) => { // First, find grid coordinates for the sixteen cells. // Produces a bicubic interpolation patch for a one-cell square. // The patch extends one cell east and south of the specified grid position, but uses 16 cells in the grid. // Deal with the edges of the input grid by duplicating adjacent values. // It's tempting to do this with typed arrays and slice(), but we need special handling for the grid edges. const x0 = (gridX === 0) ? gridX : gridX - 1 // Handle left edge const x1 = gridX const x2 = (gridX + 1 >= grid.width) ? gridX : gridX + 1 // Handle right edge const x3 = (gridX + 2 >= grid.width) ? gridX : gridX + 2 // Handle right edge const y0 = (gridY === 0) ? gridY : gridY - 1 // Handle top edge const y1 = gridY const y2 = (gridY + 1 >= grid.height) ? gridY : gridY + 1 // Handle bottom edge const y3 = (gridY + 2 >= grid.height) ? gridY : gridY + 2 // Handle bottom edge const p00 = getGridValue(grid, x0, y0) const p01 = getGridValue(grid, x0, y1) const p02 = getGridValue(grid, x0, y2) const p03 = getGridValue(grid, x0, y3) const p10 = getGridValue(grid, x1, y0) const p11 = getGridValue(grid, x1, y1) const p12 = getGridValue(grid, x1, y2) const p13 = getGridValue(grid, x1, y3) const p20 = getGridValue(grid, x2, y0) const p21 = getGridValue(grid, x2, y1) const p22 = getGridValue(grid, x2, y2) const p23 = getGridValue(grid, x2, y3) const p30 = getGridValue(grid, x3, y0) const p31 = getGridValue(grid, x3, y1) const p32 = getGridValue(grid, x3, y2) const p33 = getGridValue(grid, x3, y3) // Create interpolations through each of the four columns const columnInterpolator0 = ConstraintedSplineInterpolator(p00, p01, p02, p03) const columnInterpolator1 = ConstraintedSplineInterpolator(p10, p11, p12, p13) const columnInterpolator2 = ConstraintedSplineInterpolator(p20, p21, p22, p23) const columnInterpolator3 = ConstraintedSplineInterpolator(p30, p31, p32, p33) return function (yFraction) { // Perform curve fitting in the second (x) dimension based on the pre-fit curves in the y dimension. const p0 = columnInterpolator0(yFraction) const p1 = columnInterpolator1(yFraction) const p2 = columnInterpolator2(yFraction) const p3 = columnInterpolator3(yFraction) // Return the one-dimensional interpolator for this row. return ConstraintedSplineInterpolator(p0, p1, p2, p3) } } SplineInterpolator.tileOffset = 0.5 export default SplineInterpolator /** * Given four adjacent values a, b, c, d, fit a constrained cubic spline to them, * sacrificing smoothness to prevent overshoot. * The returned function provides interpolated values between b and c using a and d to * determine the slope going into and out of the b-c interval. * The original paper handles the general case where data points are (x,y) pairs. * In our case, the four points are always evenly spaced, so we assign X coordinates of * -1, 0, 1, 2 knowing we will perform interpolation between the second and third points. * This greatly simplifies the equations, because it gives many differences, multipliers, * and denominators have a value of 1. */ const ConstraintedSplineInterpolator = (a, b, c, d) => { const test = function (a, b, c, d) { let out = [] let interpolator = new ConstraintedSplineInterpolator (a, b, c, d) for (let i = 0; i < 100; i++) { out[i] = interpolator(i / 100) } console.log(out.join(",")) } // Optimization: if b and c are equal, interpolate a straight line if (b === c) return (x) => b const bSlope = slope(a, b, c) // equation 7a const cSlope = slope(b, c, d) // equation 7a const bSlope2 = (2 * cSlope + 4 * bSlope) + (6 * (c - b)) // equation 8 const cSlope2 = (4 * cSlope + 2 * bSlope) - (6 * (c - b)) // equation 9 const kd = (cSlope2 - bSlope2) / 6 // equation 10 const kc = bSlope2 / 2 // equation 11 const kb = (c - b) - kc - kd // equation 12 const ka = b // equation 13 // return (x) => b + (x * (c - b)) // TEST - should be exactly the same as bilinear // The returned function takes an x value in [0, 1] expressing the position between b and c, // and returns the interpolated value. return (x) => ka + kb * x + kc * x ** 2 + kd * x ** 3 } /** * Find the target slope at a particular data point. * y here is the same as y sub i in the equations in the paper. * yPrev here is y sub i-1 and yNext here is y sub i+1. * The equations are simplified significantly by the fact that we know all of our points are exactly one unit apart. */ const slope = function (yPrev, y, yNext) { const prevSlope = y - yPrev const postSlope = yNext - y if (prevSlope === 0 || postSlope === 0) return 0 // necessary condition for no overshoot if (Math.sign(prevSlope) !== Math.sign(postSlope)) return 0 // includes case where only one slope is zero return 2 / (1 / postSlope + 1 / prevSlope) } // Finding the slope at the endpoints is messy because it requries recursively computing slopes.
import {getGridValue} from '../util'

/**
 * A 2D constrained spline interpolator.
 *
 * Polynomial (e.g. bicubic) interpolation is prone to oscillation and overshoot.
 * Splines are better behaved but can still overshoot. We can sacrifice smoothness and add a no-overshoot constraint.    * See: http://www.korf.co.uk/spline.pdf
 * Kruger, CJC. Constrained Cubic Spline Interpolation for Chemical Engineering Applications.
 *
 * A pre-calculated bicubic interpolation patch.
 * For a 4x4 grid of samples, this allows us to calculate interpolated values between the central four samples.
 * By pre-fitting the curves in one dimension (y) and proceeding with the interpolation row by row,
 * we re-use most of the computation from one output pixel to the next.
 */
const SplineInterpolator = (grid, gridX, gridY) => {
  // First, find grid coordinates for the sixteen cells.
  // Produces a bicubic interpolation patch for a one-cell square.
  // The patch extends one cell east and south of the specified grid position, but uses 16 cells in the grid.

  // Deal with the edges of the input grid by duplicating adjacent values.
  // It's tempting to do this with typed arrays and slice(), but we need special handling for the grid edges. 
  const x0 = (gridX === 0) ? gridX : gridX - 1 // Handle left edge
  const x1 = gridX
  const x2 = (gridX + 1 >= grid.width) ? gridX : gridX + 1 // Handle right edge
  const x3 = (gridX + 2 >= grid.width) ? gridX : gridX + 2 // Handle right edge

  const y0 = (gridY === 0) ? gridY : gridY - 1 // Handle top edge
  const y1 = gridY
  const y2 = (gridY + 1 >= grid.height) ? gridY : gridY + 1 // Handle bottom edge
  const y3 = (gridY + 2 >= grid.height) ? gridY : gridY + 2 // Handle bottom edge

  const p00 = getGridValue(grid, x0, y0)
  const p01 = getGridValue(grid, x0, y1)
  const p02 = getGridValue(grid, x0, y2)
  const p03 = getGridValue(grid, x0, y3)

  const p10 = getGridValue(grid, x1, y0)
  const p11 = getGridValue(grid, x1, y1)
  const p12 = getGridValue(grid, x1, y2)
  const p13 = getGridValue(grid, x1, y3)

  const p20 = getGridValue(grid, x2, y0)
  const p21 = getGridValue(grid, x2, y1)
  const p22 = getGridValue(grid, x2, y2)
  const p23 = getGridValue(grid, x2, y3)

  const p30 = getGridValue(grid, x3, y0)
  const p31 = getGridValue(grid, x3, y1)
  const p32 = getGridValue(grid, x3, y2)
  const p33 = getGridValue(grid, x3, y3)

  // Create interpolations through each of the four columns
  const columnInterpolator0 = ConstraintedSplineInterpolator(p00, p01, p02, p03)
  const columnInterpolator1 = ConstraintedSplineInterpolator(p10, p11, p12, p13)
  const columnInterpolator2 = ConstraintedSplineInterpolator(p20, p21, p22, p23)
  const columnInterpolator3 = ConstraintedSplineInterpolator(p30, p31, p32, p33)

  return function (yFraction) {
    // Perform curve fitting in the second (x) dimension based on the pre-fit curves in the y dimension.
    const p0 = columnInterpolator0(yFraction)
    const p1 = columnInterpolator1(yFraction)
    const p2 = columnInterpolator2(yFraction)
    const p3 = columnInterpolator3(yFraction)
    // Return the one-dimensional interpolator for this row.
    return ConstraintedSplineInterpolator(p0, p1, p2, p3)
  }
}
SplineInterpolator.tileOffset = 0.5
export default SplineInterpolator

/**
 * Given four adjacent values a, b, c, d, fit a constrained cubic spline to them, 
 * sacrificing smoothness to prevent overshoot.
 * The returned function provides interpolated values between b and c using a and d to
 * determine the slope going into and out of the b-c interval.
 * The original paper handles the general case where data points are (x,y) pairs.
 * In our case, the four points are always evenly spaced, so we assign X coordinates of 
 * -1, 0, 1, 2 knowing we will perform interpolation between the second and third points.
 * This greatly simplifies the equations, because it gives many differences, multipliers, 
 * and denominators have a value of 1.
 */
const ConstraintedSplineInterpolator = (a, b, c, d) => {
  const test = function (a, b, c, d) {
    let out = []
    let interpolator = new ConstraintedSplineInterpolator (a, b, c, d)
    for (let i = 0; i < 100; i++) {
      out[i] = interpolator(i / 100)
    }
    console.log(out.join(","))
  }

  // Optimization: if b and c are equal, interpolate a straight line
  if (b === c) return (x) => b 
  const bSlope = slope(a, b, c) // equation 7a
  const cSlope = slope(b, c, d) // equation 7a
  const bSlope2 = (2 * cSlope + 4 * bSlope) + (6 * (c - b)) // equation 8
  const cSlope2 = (4 * cSlope + 2 * bSlope) - (6 * (c - b)) // equation 9
  const kd = (cSlope2 - bSlope2) / 6                    // equation 10
  const kc = bSlope2 / 2                                // equation 11
  const kb = (c - b) - kc - kd                            // equation 12
  const ka = b                                          // equation 13
  // return (x) => b + (x * (c - b)) // TEST - should be exactly the same as bilinear
  // The returned function takes an x value in [0, 1] expressing the position between b and c,
  // and returns the interpolated value.
  return (x) => ka + kb * x + kc * x ** 2 + kd * x ** 3
}

/**
 * Find the target slope at a particular data point.
 * y here is the same as y sub i in the equations in the paper.
 * yPrev here is y sub i-1 and yNext here is y sub i+1.
 * The equations are simplified significantly by the fact that we know all of our points are exactly one unit apart.
 */
const slope = function (yPrev, y, yNext) {
  const prevSlope = y - yPrev
  const postSlope = yNext - y
  if (prevSlope === 0 || postSlope === 0) return 0 // necessary condition for no overshoot
  if (Math.sign(prevSlope) !== Math.sign(postSlope)) return 0 // includes case where only one slope is zero
  return 2 / (1 / postSlope + 1 / prevSlope)
}
// Finding the slope at the endpoints is messy because it requries recursively computing slopes.

二、单调Cubic 插值

下面这个demo可以保证绝对单调

MonotonicCubicPy/monotonic_cubic.py at master · bssrdf/MonotonicCubicPy (github.com)

另一种常用的方法是pchip

multivariate-c2/pchip.cpp at master · vasily-kartashov/multivariate-c2 · GitHub

发表评论