Skip to content

Instantly share code, notes, and snippets.

@onsummer

onsummer/bilinear-interpolation.js

Last active March 7, 2022 02:48
Show Gist options
  • Save onsummer/7c01d7d624fc9270e7454df267996755 to your computer and use it in GitHub Desktop.
Save onsummer/7c01d7d624fc9270e7454df267996755 to your computer and use it in GitHub Desktop.
Download ZIP
bilinear-interpolation
Raw
bilinear-interpolation.js
/**
* 计算双线性内插点的值
* @param {Point} pt 目标点位
* @param {Point} p1 左上
* @param {Point} p2 右上
* @param {Point} p3 左下
* @param {Point} p4 右下
*
* @returns {Number} pt的插值结果
*/
const bilinearInterpolation2D = (pt, p1, p2, p3, p4) => {
const xt1 = linearInterpolation2D(pt, p1, p2)
const xt2 = linearInterpolation2D(pt, p3, p4)
const val = ((p1.y - pt.y) * xt1 + (pt.y - p3.y) * xt2) / (p1.y - p3.y)
pt.v = val
return val
}
/**
* 计算双线性内插点的值
*
* pt[待插值变量] = a·(b·p1[待插值的变量] + c·p2[待插值的变量] + d·p3[待插值的变量] + e·p4[待插值的变量])
*
* 其中
* a = 1/((p2.x - p4.x)·(p2.y - p4.y))
* b = (p2.x - pt.x)(p2.y - pt.y)
* c = (pt.x - p1.x)(p1.y - pt.y)
* d = (p3.x - pt.x)(pt.y - p3.y)
* e = (pt.x - p4.x)(pt.y - p4.y)
*
* @param {Point} pt 目标点位
* @param {Point} p1 左上
* @param {Point} p2 右上
* @param {Point} p3 左下
* @param {Point} p4 右下
*/
const bilinearInterpolation2D_2 = (pt, p1, p2, p3, p4) => {
const a = (p2.x - p4.x) * (p2.y - p4.y)
const b = (p2.x - pt.x) * (p2.y - pt.y)
const c = (pt.x - p1.x) * (p1.y - pt.y)
const d = (p3.x - pt.x) * (pt.y - p3.y)
const e = (pt.x - p4.x) * (pt.y - p4.y)
return (b * p4.v + c * p3.v + d * p2.v + e * p1.v) / a
}
/**
* 线性插值
* @param {Point} pt 待插值点
* @param {Point} p1 左点
* @param {Point} p2 右点
* @param {String} order 顺序，可以是 'x' 或 'y'，默认是 'x'
*/
const linearInterpolation2D = (pt, p1, p2, order = 'x') => {
if (order === 'x' || order === 'y') {
const k1 = (p2[order] - pt[order]) / (p2[order] - p1[order])
const k2 = (pt[order] - p1[order]) / (p2[order] - p1[order])
return k1 * p1.v + k2 * p2.v
}
else {
return undefined
}
}
/**
* 逆双线性插值算法，计算凸四边形内任意一点的插值。
*
* 插值公式
*
* value_pt =
* value_p1 +
* (value_p2 - value_p1) * u +
* (value_p4 - value_p1) * v +
* (value_p1 - value_p2 + value_p3 - value_p4) * uv
*
* 而 u、v 则满足如下方程组
*
* · av^2 + bv + c = 0
* · u = [ (pt.x - p1.x) - (p4.x - p1.x) * v ] / [ (p2.x - p1.x) + (p1.x - p2.x + p3.x - p4.x) * v ]
*
* 其中，a、b、c 参数可由坐标计算而来
*
* · a = (p1.x - p2.x + p3.x - p4.x) * (p4.y - p1.y) - (p1.y - p2.y + p3.y - p4.y) * (p4.x - p1.x)
* · b = (p2.x - p1.x) * (p4.y - p1.y) - (p2.y - p1.y) * (p4.x - p1.x)
* + (pt.x - p1.x) * (p1.y - p2.y + p3.y - p4.y) - (pt.y - p1.y) * (p1.x - p2.x + p3.x - p4.x)
* · c = (pt.x - p1.x) * (p2.y - p1.y) - (pt.y - p1.y) * (p2.x - p1.x)
*
* 求根公式可得到 v 的解，取 [0,1] 区间内的值即可。
*
* @todo
* 注意，若 p1p2 和 p3p4 共线，则 a = 0，即解一元一次方程 bv + c = 0 => v = -c/b.
*
* @see https://blog.csdn.net/lk040384/article/details/104939742
*
* @param {Point} pt 待插值 2D 点
* @param {Point} p1 点1
* @param {Point} p2 点2
* @param {Point} p3 点3
* @param {Point} p4 点4
*/
const inverseBilinearInterpolation2D = (pt, p1, p2, p3, p4) => {
if (isPointWithinQuadrilateral(pt, p1, p2, p3, p4) === false) {
return
}
let u = 0, v = 0
const a = (p1.x - p2.x + p3.x - p4.x) * (p4.y - p1.y) - (p1.y - p2.y + p3.y - p4.y) * (p4.x - p1.x)
const b = (p2.x - p1.x) * (p4.y - p1.y) - (p2.y - p1.y) * (p4.x - p1.x) + (pt.x - p1.x) * (p1.y - p2.y + p3.y - p4.y) - (pt.y - p1.y) * (p1.x - p2.x + p3.x - p4.x)
const c = (pt.x - p1.x) * (p2.y - p1.y) - (pt.y - p1.y) * (p2.x - p1.x)
const delta = b ** 2 - 4 * a * c
if (delta < 0) {
throw new RangeError('delta is smaller than zero.')
}
const sqrtDelta = Math.sqrt(delta)
/** calc v */
v = (-b + sqrtDelta) / (2 * a)
if (v < 0 || v > 1) {
v = (-b - sqrtDelta) / (2 * a)
}
/** calc u */
u = ((pt.x - p1.x) - (p4.x - p1.x) * v) / ((p2.x - p1.x) + (p1.x - p2.x + p3.x - p4.x) * v)
/** calc value of pt */
const val = p1.v + (p2.v - p1.v) * u + (p4.v - p1.v) * v + (p1.v - p2.v + p3.v - p4.v) * u * v
pt.v = val
return val
}
/**
* 判断点是否在四边形内，使用顺时针逆时针算法：若均为顺时针，则点在四边形内
* 其中，p1、p2、p3、p4 为四边形顺时针的四个顶点
*
* p1 -———> p2
* ↑ ↘↖ ↙↗ |
* | pt |
* | ↙↗ ↘↖ ↓
* p4 <-——— p3
*
* @param {Point} pt 待判断的 2D 点
* @param {Point} p1 四边形的点1
* @param {Point} p2 四边形的点2
* @param {Point} p3 四边形的点3
* @param {Point} p4 四边形的点4
* @returns {Boolean}
*/
const isPointWithinQuadrilateral = (pt, p1, p2, p3, p4) => {
const booleanArr = []
booleanArr.push(isClockWise(p1, p2, pt))
booleanArr.push(isClockWise(p2, p3, pt))
booleanArr.push(isClockWise(p3, p4, pt))
booleanArr.push(isClockWise(p4, p1, pt))
return booleanArr.every(v => v === true)
}
/**
* 判断三个 2D 点的顺序是否为顺时针，是则返回 true，否则返回 false。
*
* p1
* ↗ ↘
* p3 ← p2
*
* @param {Point} p1 点1
* @param {Point} p2 点2
* @param {Point} p3 点3
* @returns {Boolean}
*/
const isClockWise = (p1, p2, p3) => {
return (p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y) < 0
}
/**
*
* @param {Point} pt 待测试点
* @param {Point} p1 点1
* @param {Point} p2 点2
* @param {Point} p3 点3
* @param {Point} p4 点4
*
* @returns {Boolean[][]} 若返回 [[1,2]] 则代表pt在1、2点外，若返回[[1,2],[2,3]]则代表在1、2、3点外
*/
const onWhichSide = (pt, p1, p2, p3, p4) => {
const flags = [
isClockWise(p1, p2, pt),
isClockWise(p2, p3, pt),
isClockWise(p3, p4, pt),
isClockWise(p4, p1, pt),
]
let sidePointIds = []
flags.forEach((v, i, _) => {
if (v == false) {
// 若i是3，i+2必定是5，那么得让它i-2才行
sidePointIds.push([i + 1, i + 2 > flags.length ? i - 2 : i + 2])
}
})
return sidePointIds
}
export {
bilinearInterpolation2D,
bilinearInterpolation2D_2,
inverseBilinearInterpolation2D,
isPointWithinQuadrilateral,
isClockWise,
onWhichSide
}
Raw
bilinear-interpolation.ts
interface IPoint {
x: Number,
y: Number,
v: Number
}
/**
* 计算双线性内插点的值
* @param {IPoint} pt 目标点位
* @param {IPoint} p1 左上
* @param {IPoint} p2 右上
* @param {IPoint} p3 左下
* @param {IPoint} p4 右下
*
* @returns {Number} pt的插值结果
*/
const bilinearInterpolation2D = (
pt: IPoint,
p1: IPoint, p2: IPoint, p3: IPoint, p4: IPoint): number => {
const xt1 = linearInterpolation2D(pt, p1, p2)
const xt2 = linearInterpolation2D(pt, p3, p4)
const val = ((p1.y - pt.y) * xt1 + (pt.y - p3.y) * xt2) / (p1.y - p3.y)
pt.v = val
return val
}
/**
* 计算双线性内插点的值
*
* pt[待插值变量] = a·(b·p1[待插值的变量] + c·p2[待插值的变量] + d·p3[待插值的变量] + e·p4[待插值的变量])
*
* 其中
* a = 1/((p2.x - p4.x)·(p2.y - p4.y))
* b = (p2.x - pt.x)(p2.y - pt.y)
* c = (pt.x - p1.x)(p1.y - pt.y)
* d = (p3.x - pt.x)(pt.y - p3.y)
* e = (pt.x - p4.x)(pt.y - p4.y)
*
* @param {Point} pt 目标点位
* @param {Point} p1 左上
* @param {Point} p2 右上
* @param {Point} p3 左下
* @param {Point} p4 右下
*/
const bilinearInterpolation2D_2 = (
pt: IPoint,
p1: IPoint, p2: IPoint, p3: IPoint, p4: IPoint): number => {
const a = (p2.x - p4.x) * (p2.y - p4.y)
const b = (p2.x - pt.x) * (p2.y - pt.y)
const c = (pt.x - p1.x) * (p1.y - pt.y)
const d = (p3.x - pt.x) * (pt.y - p3.y)
const e = (pt.x - p4.x) * (pt.y - p4.y)
return (b * p4.v + c * p3.v + d * p2.v + e * p1.v) / a
}
/**
* 线性插值
* @param {IPoint} pt 待插值点
* @param {IPoint} p1 左点
* @param {IPoint} p2 右点
* @param {String} order 顺序，可以是 'x' 或 'y'，默认是 'x'
*/
const linearInterpolation2D = (pt: IPoint, p1: IPoint, p2: IPoint, order: string = 'x'): number | undefined => {
if (order === 'x' || order === 'y') {
const k1 = (p2[order] - pt[order]) / (p2[order] - p1[order])
const k2 = (pt[order] - p1[order]) / (p2[order] - p1[order])
return k1 * p1.v + k2 * p2.v
}
else {
return undefined
}
}
/**
* 逆双线性插值算法，计算凸四边形内任意一点的插值。
*
* 插值公式
*
* value_pt =
* value_p1 +
* (value_p2 - value_p1) * u +
* (value_p4 - value_p1) * v +
* (value_p1 - value_p2 + value_p3 - value_p4) * uv
*
* 而 u、v 则满足如下方程组
*
* · av^2 + bv + c = 0
* · u = [ (pt.x - p1.x) - (p4.x - p1.x) * v ] / [ (p2.x - p1.x) + (p1.x - p2.x + p3.x - p4.x) * v ]
*
* 其中，a、b、c 参数可由坐标计算而来
*
* · a = (p1.x - p2.x + p3.x - p4.x) * (p4.y - p1.y) - (p1.y - p2.y + p3.y - p4.y) * (p4.x - p1.x)
* · b = (p2.x - p1.x) * (p4.y - p1.y) - (p2.y - p1.y) * (p4.x - p1.x)
* + (pt.x - p1.x) * (p1.y - p2.y + p3.y - p4.y) - (pt.y - p1.y) * (p1.x - p2.x + p3.x - p4.x)
* · c = (pt.x - p1.x) * (p2.y - p1.y) - (pt.y - p1.y) * (p2.x - p1.x)
*
* 求根公式可得到 v 的解，取 [0,1] 区间内的值即可。
*
* @todo
* 注意，若 p1p2 和 p3p4 共线，则 a = 0，即解一元一次方程 bv + c = 0 => v = -c/b.
*
* @see https://blog.csdn.net/lk040384/article/details/104939742
*
* @param {IPoint} pt 待插值 2D 点
* @param {IPoint} p1 点1
* @param {IPoint} p2 点2
* @param {IPoint} p3 点3
* @param {IPoint} p4 点4
*/
const inverseBilinearInterpolation2D = (
pt: IPoint,
p1: IPoint, p2: IPoint, p3: IPoint, p4: IPoint
): number | undefined => {
if (isPointWithinQuadrilateral(pt, p1, p2, p3, p4) === false) {
return undefined
}
let u = 0, v = 0
const a = (p1.x - p2.x + p3.x - p4.x) * (p4.y - p1.y) - (p1.y - p2.y + p3.y - p4.y) * (p4.x - p1.x)
const b = (p2.x - p1.x) * (p4.y - p1.y) - (p2.y - p1.y) * (p4.x - p1.x) + (pt.x - p1.x) * (p1.y - p2.y + p3.y - p4.y) - (pt.y - p1.y) * (p1.x - p2.x + p3.x - p4.x)
const c = (pt.x - p1.x) * (p2.y - p1.y) - (pt.y - p1.y) * (p2.x - p1.x)
const delta = b ** 2 - 4 * a * c
if (delta < 0) {
throw new RangeError('delta is smaller than zero.')
}
const sqrtDelta = Math.sqrt(delta)
/** calc v */
v = (-b + sqrtDelta) / (2 * a)
if (v < 0 || v > 1) {
v = (-b - sqrtDelta) / (2 * a)
}
/** calc u */
u = ((pt.x - p1.x) - (p4.x - p1.x) * v) / ((p2.x - p1.x) + (p1.x - p2.x + p3.x - p4.x) * v)
/** calc value of pt */
const val = p1.v + (p2.v - p1.v) * u + (p4.v - p1.v) * v + (p1.v - p2.v + p3.v - p4.v) * u * v
pt.v = val
return val
}
/**
* 判断点是否在四边形内，使用顺时针逆时针算法：若均为顺时针，则点在四边形内
* 其中，p1、p2、p3、p4 为四边形顺时针的四个顶点
*
* p1 ———————> p2
* ↑ ↘↖ ↙↗ |
* | pt |
* | ↙↗ ↘↖ ↓
* p4 <——————— p3
*
* @param {IPoint} pt 待判断的 2D 点
* @param {IPoint} p1 四边形的点1
* @param {IPoint} p2 四边形的点2
* @param {IPoint} p3 四边形的点3
* @param {IPoint} p4 四边形的点4
* @returns {Boolean}
*/
const isPointWithinQuadrilateral = (
pt: IPoint,
p1: IPoint, p2: IPoint, p3: IPoint, p4: IPoint
): boolean => {
const booleanArr = []
booleanArr.push(isClockWise(p1, p2, pt))
booleanArr.push(isClockWise(p2, p3, pt))
booleanArr.push(isClockWise(p3, p4, pt))
booleanArr.push(isClockWise(p4, p1, pt))
return booleanArr.every(v => v === true)
}
/**
* 判断三个 2D 点的顺序是否为顺时针，是则返回 true，否则返回 false。
*
* p1
* ↗ ↘
* p3 ← p2
*
* @param {IPoint} p1 点1
* @param {IPoint} p2 点2
* @param {IPoint} p3 点3
* @returns {Boolean}
*/
const isClockWise = (
p1: IPoint, p2: IPoint, p3: IPoint
): boolean => {
return (p2.x - p1.x) * (p3.y - p1.y) - (p3.x - p1.x) * (p2.y - p1.y) < 0
}
/**
*
* @param {IPoint} pt 待测试点
* @param {IPoint} p1 点1
* @param {IPoint} p2 点2
* @param {IPoint} p3 点3
* @param {IPoint} p4 点4
*
* @returns {Boolean[][]} 若返回 [[1,2]] 则代表pt在1、2点外，若返回[[1,2],[2,3]]则代表在1、2、3点外
*/
const onWhichSide = (
pt: IPoint,
p1: IPoint, p2: IPoint, p3: IPoint, p4: IPoint
): bool[][] => {
const flags = [
isClockWise(p1, p2, pt),
isClockWise(p2, p3, pt),
isClockWise(p3, p4, pt),
isClockWise(p4, p1, pt),
]
let sidePointIds = []
flags.forEach((v, i, _) => {
if (v == false) {
// 若i是3，i+2必定是5，那么得让它i-2才行
sidePointIds.push([i + 1, i + 2 > flags.length ? i - 2 : i + 2])
}
})
return sidePointIds
}
export {
bilinearInterpolation2D,
bilinearInterpolation2D_2,
inverseBilinearInterpolation2D,
isPointWithinQuadrilateral,
isClockWise,
onWhichSide
}
Raw
random-point-generator.js
function randomPoint2D(count, zoom = [-100, 100], numberFixed = 4, handler = result => result) {
let result = []
if (count < 1) {
return result
}
const scale = zoom[1] - zoom[0]
for (let i = 0; i < count; i++) {
result.push([
round(zoom[0] + scale * Math.random(), numberFixed),
round(zoom[0] + scale * Math.random(), numberFixed)
])
}
/** 加入自定义处理机制 */
result = handler(result)
return result
}
function round(value, n) {
if (isNaN(value)) {
return 0
}
const p1 = Math.pow(10, n + 1)
const p2 = Math.pow(10, n)
return Math.round(value * p1 / 10) / p2
}
export {
round,
randomPoint2D
}
Raw
test-data
{x: -6.29, y: 2.45} {x: 5.68, y: 3.42}
{x: -2.33, y: 1.15}
{x: 4.12, y: -2.77}
{x: -5.94, y: -4.89}
params:
pt, p1, p2, p3, p4
{x: -2.33, y: 1.15, v: 0}, {x: -6.29, y: 2.45, v: 1}, {x: 5.68, y: 3.42, v: 2}, {x: 4.12, y: -2.77, v: 3}, {x: -5.94, y: -4.89, v: 2}
without value:
pt, p1, p2, p3, p4
{x: -2.33, y: 1.15}, {x: -6.29, y: 2.45}, {x: 5.68, y: 3.42}, {x: 4.12, y: -2.77}, {x: -5.94, y: -4.89}
一个点在1、4号点的外围
外围点：{x: -10, y: 0}
{x: -8.0336, y: 6.1536}
{x: 4.0793, y: 0.5787}
{x: -9.0849, y: -2.0879}
{x: -4.3149, y: -4.833}
双线性插值点
{x: 17, y: 30, v: 2} {x: 24, y: 30, v: 2}
{x: 19, y: 25, v: 0}
{x: 17, y: 21, v: 2} {x: 24, y: 21, v: 2}
{x: 19, y: 25, v: 0}, {x: 17, y: 30, v: 2}, {x: 24, y: 30, v: 2}, {x: 24, y: 21, v: 2}, {x: 17, y: 21, v: 2}
Raw
test-pointside-lonlatdemo.js
import { onWhichSide } from './src/fn.js'
import { randomPoint2D } from './src/utils/random-point-generator.js'
function Point(x, y, v) {
this.x = x
this.y = y
this.v = v
}
// const randomPoints = randomPoint2D(4, [-10, 10], 4, result => {
// result.forEach((v, i, arr) => {
// arr[i] = new Point(v[0], v[1], 0)
// })
// return result
// })
/* test group 1
const p1 = new Point(-8.0336, 6.1536, 1)
const p2 = new Point(4.0793, 0.5787, 2)
const p3 = new Point(-4.3149, -4.833, 3)
const p4 = new Point(-9.0849, -2.0879, 4)
const pout = new Point(-10, 0, 0)
*/
const p1 = new Point(114.12931061, 22.43935585, 1)
const p2 = new Point(114.13014221, 22.43938446, 2)
const p3 = new Point(114.13018036, 22.43873405, 3)
const p4 = new Point(114.12931824, 22.43870163, 4)
const pout = new Point(114.1289873, 22.4397613, 0)
const sidePtIds = onWhichSide(pout, p1, p2, p3, p4)
//#region
// 若 sidePtIds.length == 0，那么进行插值
/* 若不为0，则：
[1,2] -> 2：整体行列号上移1格
[2,3] -> 6: 整体行列号右移1格
[3,4] -> 12: 整体行列号下移1格
[4,1] -> 4: 整体行列号左移1格
*/
//#endregion
sidePtIds.forEach(group => {
const tag = group[0] * group[1]
let orient = undefined
if (tag === 2)
orient = '上'
else if (tag === 6)
orient = '右'
else if (tag === 12)
orient = '下'
else if (tag === 4)
orient = '左'
else {
inverseBilinearInterpolation2D(pout, p1, p2, p3, p4)
console.log(`插值结果：${JSON.stringify(pout)}`)
}
console.log(`点位于第 ${group[0]} 个点和第 ${group[1]} 个点的 ${orient} 侧。`)
})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment