线性插值

已知两点

A(x_1, y_1)

和

B(x_2, y_2)

，以及插值参数

t (0 \leq t \leq 1）

，可以通过以下公式计算任意中间点

C(x_t, y_t)

的坐标。

x 坐标的插值：
中间点的 x 坐标 $x_t$ 是 $A$ 和 $B$ 的 $x$ 坐标按照比例 $t$ 插值的结果，公式为：
$x_t = (1 - t) \cdot x_1 + t \cdot x_2$
- 当 $t = 0$ 时， $x_t = x_1$ 。
- 当 $t = 1$ 时， $x_t = x_2$ 。
- 在 $0 \leq t \leq 1$ 的范围内， $x_t$ 线性分布在 $x_1$ 和 $x_2$ 之间。
y 坐标的插值：
中间点的 $y$ 坐标 $y_t$ 是 $A$ 和 $B$ 的 $y$ 坐标按照比例 $t$ 插值的结果，公式为：
$y_t = (1 - t) \cdot y_1 + t \cdot y_2$
- 当 $t = 0$ 时， $y_t = y_1$ 。
- 当 $t = 1$ 时， $y_t = y_2$ 。
- 在 $0 \leq t \leq 1$ 的范围内， $y_t$ 线性分布在 $y_1$ 和 $y_2$ 之间。
y 坐标的插值：
中间点的 y 坐标 $y_t$ 是 $A$ 和 $B$ 的 y 坐标按照比例 $t$ 插值的结果，公式为：
$y_t = (1 - t) \cdot y_1 + t \cdot y_2$
综合表示：
中间点 $C(x_t, y_t)$ 的坐标公式为：
$C(x_t, y_t) = \big((1 - t) \cdot x_1 + t \cdot x_2, \, (1 - t) \cdot y_1 + t \cdot y_2 \big)$

双线性插值 bilinear

双线性插值

双线性插值公式推导

双线性插值是一种常见的插值方法，主要用于二维网格中根据已知点的值来估算某个未知点的值。

1. 问题描述

给定四个已知网格点 $Q_0, Q_1, Q_2, Q_3$ 的值 $z_{00}, z_{10}, z_{01}, z_{11}$ ，以及目标插值点 $P(tx, ty)$ 的相对位置，要求插值点 $P(tx, ty)$ 的值 $z$ 。在这里刚好 $m=tx, n=ty$ ，一般来说 $tx,ty$ 是 $m,n$ 的小数部分

四个网格点的坐标分别为：
$Q_0: (0, 0), \quad Q_1: (1, 0), \quad Q_2: (0, 1), \quad Q_3: (1, 1)$
插值点 $P(tx, ty)$ 的坐标为：
$P(tx, ty), \quad 其中 tx, ty \in [0, 1]$

目标是根据 $tx, ty$ 对 $z$ 进行估算。

2. 计算 $P_1$ 和 $P_2$

步骤 1：对 $x$ 方向插值

在 $Q_0$ 和 $Q_1$ 之间插值，计算 $P_1(tx, 0)$ 的值：
$P_1(tx, 0) = z_{00} \cdot (1 - tx) + z_{10} \cdot tx$

同理，在 $Q_2$ 和 $Q_3$ 之间插值，计算 $P_2(tx, 1)$ 的值：
$P_2(tx, 1) = z_{01} \cdot (1 - tx) + z_{11} \cdot tx$

3. 根据 $P_1$ 和 $P_2$ 计算 $P(tx, ty)$

步骤 2：对 $y$ 方向插值

在 $P_1(tx, 0)$ 和 $P_2(tx, 1)$ 之间插值，计算 $P(tx, ty)$ 的值：
$P(tx, ty) = P_1(tx, 0) \cdot (1 - ty) + P_2(tx, 1) \cdot ty$

将 $P_1(tx, 0)$ 和 $P_2(tx, 1)$ 的表达式代入上式，得到：
$P(tx, ty) = \big[ z_{00} \cdot (1 - tx) + z_{10} \cdot tx \big] \cdot (1 - ty) + \big[ z_{01} \cdot (1 - tx) + z_{11} \cdot tx \big] \cdot ty$

展开化简：
$P(tx, ty) = z_{00} \cdot (1 - tx) \cdot (1 - ty) + z_{10} \cdot tx \cdot (1 - ty) + z_{01} \cdot (1 - tx) \cdot ty + z_{11} \cdot tx \cdot ty$

4. 双线性插值公式

最终公式为：
$P(tx, ty) = z_{00} \cdot (1 - tx) \cdot (1 - ty) + z_{10} \cdot tx \cdot (1 - ty) + z_{01} \cdot (1 - tx) \cdot ty + z_{11} \cdot tx \cdot ty$

5. 总结

双线性插值公式分两步：

$x$ 方向插值，计算 $P_1$ 和 $P_2$ 的值。
$y$ 方向插值，根据 $P_1$ 和 $P_2$ 的值计算目标点 $P(tx, ty)$ 。

这是一种对二维网格点的插值方法，假设 $x$ 和 $y$ 方向是线性变化的，因此公式基于线性插值的原则推导而来。

用面积理解双线性插值示意图

面积理解双线性插值

双三次插值 bicubic

双三次插值

步骤 1：计算 $x$ 方向的中间插值点 $P_0, P_1, P_2, P_3$

确定插值点和网格点
- 插值点为 $P(n, m)$ ，位于 $(n, m)$ 。
- 最近的左上角网格点为 $Q_{11} = (\lfloor n \rfloor, \lfloor m \rfloor)$ 。
- 选择 $P$ 点周围的 $4 \times 4$ 网格点：
  $q_{i,j} = (\lfloor n \rfloor + j, \lfloor m \rfloor + i), \quad i, j \in \{-1, 0, 1, 2\}$
计算 $x$ 方向的偏移
- 计算 $t_x = n - \lfloor n \rfloor$ ，即插值点在 $x$ 方向上相对于网格点的偏移。
计算 $x$ 方向的中间值 $P_0, P_1, P_2, P_3$
- 对于每一行 $i \in \{-1, 0, 1, 2\}$ ，使用 $x$ 方向的插值权重 $w_x(j, t_x)$ 计算中间插值点：
  $P_i = \sum_{j=-1}^{2} w_x(j, t_x) \cdot q_{i,j}$
- $w_x(j, t_x)$ 的权重函数 $w(x)$ 定义如下（见权重公式部分）。

步骤 2：计算 $y$ 方向的最终插值点 $P$

计算 $y$ 方向的偏移
- 计算 $t_y = m - \lfloor m \rfloor$ ，即插值点在 $y$ 方向上相对于网格点的偏移。
计算 $y$ 方向的插值值 $P$
- 使用 $y$ 方向的插值权重 $w_y(i, t_y)$ 和中间插值点 $P_0, P_1, P_2, P_3$ 计算最终插值点的值：
  $P(n, m) = \sum_{i=-1}^{2} w_y(i, t_y) \cdot P_i$

双三次插值的权重函数

权重函数 $w(x)$ 的定义为：
$w(x) = \begin{cases} 1.5 |x|^3 - 2.5 |x|^2 + 1, & \text{if } 0 \leq |x| < 1 \\ -0.5 |x|^3 + 2.5 |x|^2 - 4 |x| + 2, & \text{if } 1 \leq |x| < 2 \\ 0, & \text{if } |x| \geq 2 \end{cases}$

权重函数中的变量说明

$x$ 表示插值点与网格点之间的相对距离，定义为：
$x = j - t_x \quad \text{或} \quad x = i - t_y$
根据 $x$ 的绝对值范围，选择对应的公式计算权重。

总结步骤

计算 $x$ 方向的中间插值点
- 使用 $t_x$ 和 $w(x)$ 对每一行进行插值，得到 $P_0, P_1, P_2, P_3$ 。
计算 $y$ 方向的最终插值点
- 使用 $t_y$ 和 $w(x)$ 对 $P_0, P_1, P_2, P_3$ 进行插值，得到插值点 $P$ 。
核心公式
- $x$ 方向中间值：
  $P_i = \sum_{j=-1}^{2} w_x(j, t_x) \cdot q_{i,j}$
- $y$ 方向最终值：
  $P(n, m) = \sum_{i=-1}^{2} w_y(i, t_y) \cdot P_i$
权重函数：
$w(x) = \begin{cases} 1.5 |x|^3 - 2.5 |x|^2 + 1, & \text{if } 0 \leq |x| < 1 \\ -0.5 |x|^3 + 2.5 |x|^2 - 4 |x| + 2, & \text{if } 1 \leq |x| < 2 \\ 0, & \text{if } |x| \geq 2 \end{cases}$

cuda 仿射变换实现

双线性插值

static .,__device__ float bicubic_weight(float x) {
    const float a = -0.5f;  // 常见的值，通常取-0.5f
    if (x < 0) x = -x;
    if (x < 1) return (a + 2.0f) * powf(x, 3) - (a + 3.0f) * powf(x, 2) + 1.0f;
    if (x < 2) return a * powf(x, 3) - 5.0f * a * powf(x, 2) + 8.0f * a * x - 4.0f * a;
    return 0.0f;
}

static __global__ void warp_affine_bicubic_plane_kernel(uint8_t *src,
                                                        int src_line_size,
                                                        int src_width,
                                                        int src_height,
                                                        uint8_t *dst,
                                                        int dst_width,
                                                        int dst_height,
                                                        uint8_t const_value_st,
                                                        float *warp_affine_matrix_2_3)
{
    int dx = blockDim.x * blockIdx.x + threadIdx.x;
    int dy = blockDim.y * blockIdx.y + threadIdx.y;
    if (dx >= dst_width || dy >= dst_height)
        return;

    float m_x1 = warp_affine_matrix_2_3[0];
    float m_y1 = warp_affine_matrix_2_3[1];
    float m_z1 = warp_affine_matrix_2_3[2];
    float m_x2 = warp_affine_matrix_2_3[3];
    float m_y2 = warp_affine_matrix_2_3[4];
    float m_z2 = warp_affine_matrix_2_3[5];

    float src_x = m_x1 * dx + m_y1 * dy + m_z1;
    float src_y = m_x2 * dx + m_y2 * dy + m_z2;

    if (src_x <= 1 || src_x >= src_width - 2 || src_y <= 1 || src_y >= src_height - 2)
    {
        dst[(dy * dst_width + dx) * 3] = const_value_st;
        dst[(dy * dst_width + dx) * 3 + 1] = const_value_st;
        dst[(dy * dst_width + dx) * 3 + 2] = const_value_st;
    }
    else
    {
        // 直接在插值时获取像素值
        float wts_x[4], wts_y[4];
        for (int i = 0; i < 4; i++)
        {
            wts_x[i] = bicubic_weight(src_x - ((int)floorf(src_x) - 1 + i));
            wts_y[i] = bicubic_weight(src_y - ((int)floorf(src_y) - 1 + i));
        }

        // x 方向插值
        float p_x[4][3] = {{0}};
        for (int i = 0; i < 4; i++) 
        {
            float px_r = 0, px_g = 0, px_b = 0;
            for (int j = 0; j < 4; j++) 
            {
                int y = (int)floorf(src_y) - 1 + i;
                int x = (int)floorf(src_x) - 1 + j;
                // 边界检查并计算像素
                if (y >= 0 && y < src_height && x >= 0 && x < src_width)
                {
                    uint8_t* pixel = &src[y * src_line_size + x * 3];
                    px_r += pixel[0] * bicubic_weight(src_x - x);
                    px_g += pixel[1] * bicubic_weight(src_x - x);
                    px_b += pixel[2] * bicubic_weight(src_x - x);
                }
                else
                {
                    px_r += const_value_st * bicubic_weight(src_x - x);
                    px_g += const_value_st * bicubic_weight(src_x - x);
                    px_b += const_value_st * bicubic_weight(src_x - x);
                }
            }
            p_x[i][0] = px_r;
            p_x[i][1] = px_g;
            p_x[i][2] = px_b;
        }

        // y 方向插值并输出到目标图像
        dst[(dy * dst_width + dx) * 3]     = (uint8_t)min(max(wts_y[0] * p_x[0][0] + wts_y[1] * p_x[1][0] + wts_y[2] * p_x[2][0] + wts_y[3] * p_x[3][0] + 0.5f, 0.0f), 255.0f);
        dst[(dy * dst_width + dx) * 3 + 1] = (uint8_t)min(max(wts_y[0] * p_x[0][1] + wts_y[1] * p_x[1][1] + wts_y[2] * p_x[2][1] + wts_y[3] * p_x[3][1] + 0.5f, 0.0f), 255.0f);
        dst[(dy * dst_width + dx) * 3 + 2] = (uint8_t)min(max(wts_y[0] * p_x[0][2] + wts_y[1] * p_x[1][2] + wts_y[2] * p_x[2][2] + wts_y[3] * p_x[3][2] + 0.5f, 0.0f), 255.0f);
    }
}

双三次插值

static __device__ float bicubic_weight(float x) {
    const float a = -0.5f;  // 常见的值，通常取-0.5f
    if (x < 0) x = -x;
    if (x < 1) return (a + 2.0f) * powf(x, 3) - (a + 3.0f) * powf(x, 2) + 1.0f;
    if (x < 2) return a * powf(x, 3) - 5.0f * a * powf(x, 2) + 8.0f * a * x - 4.0f * a;
    return 0.0f;
}

static __global__ void warp_affine_bicubic_plane_kernel(uint8_t *src,
                                                        int src_line_size,
                                                        int src_width,
                                                        int src_height,
                                                        uint8_t *dst,
                                                        int dst_width,
                                                        int dst_height,
                                                        uint8_t const_value_st,
                                                        float *warp_affine_matrix_2_3)
{
    int dx = blockDim.x * blockIdx.x + threadIdx.x;
    int dy = blockDim.y * blockIdx.y + threadIdx.y;
    if (dx >= dst_width || dy >= dst_height)
        return;

    float m_x1 = warp_affine_matrix_2_3[0];
    float m_y1 = warp_affine_matrix_2_3[1];
    float m_z1 = warp_affine_matrix_2_3[2];
    float m_x2 = warp_affine_matrix_2_3[3];
    float m_y2 = warp_affine_matrix_2_3[4];
    float m_z2 = warp_affine_matrix_2_3[5];

    float src_x = m_x1 * dx + m_y1 * dy + m_z1;
    float src_y = m_x2 * dx + m_y2 * dy + m_z2;

    if (src_x <= 1 || src_x >= src_width - 2 || src_y <= 1 || src_y >= src_height - 2)
    {
        dst[(dy * dst_width + dx) * 3] = const_value_st;
        dst[(dy * dst_width + dx) * 3 + 1] = const_value_st;
        dst[(dy * dst_width + dx) * 3 + 2] = const_value_st;
    }
    else
    {
        // 直接在插值时获取像素值
        float wts_y[4];
        for (int i = 0; i < 4; i++)
        {
            wts_y[i] = bicubic_weight(src_y - ((int)floorf(src_y) - 1 + i));
        }

        // x 方向插值
        float p_x[4][3] = {{0}};
        for (int i = 0; i < 4; i++) 
        {
            float px_r = 0, px_g = 0, px_b = 0;
            for (int j = 0; j < 4; j++) 
            {
                int y = (int)floorf(src_y) - 1 + i;
                int x = (int)floorf(src_x) - 1 + j;
                // 边界检查并计算像素
                if (y >= 0 && y < src_height && x >= 0 && x < src_width)
                {
                    uint8_t* pixel = &src[y * src_line_size + x * 3];
                    px_r += pixel[0] * bicubic_weight(src_x - x);
                    px_g += pixel[1] * bicubic_weight(src_x - x);
                    px_b += pixel[2] * bicubic_weight(src_x - x);
                }
                else
                {
                    px_r += const_value_st * bicubic_weight(src_x - x);
                    px_g += const_value_st * bicubic_weight(src_x - x);
                    px_b += const_value_st * bicubic_weight(src_x - x);
                }
            }
            p_x[i][0] = px_r;
            p_x[i][1] = px_g;
            p_x[i][2] = px_b;
        }

        // y 方向插值并输出到目标图像
        dst[(dy * dst_width + dx) * 3]     = (uint8_t)min(max(wts_y[0] * p_x[0][0] + wts_y[1] * p_x[1][0] + wts_y[2] * p_x[2][0] + wts_y[3] * p_x[3][0] + 0.5f, 0.0f), 255.0f);
        dst[(dy * dst_width + dx) * 3 + 1] = (uint8_t)min(max(wts_y[0] * p_x[0][1] + wts_y[1] * p_x[1][1] + wts_y[2] * p_x[2][1] + wts_y[3] * p_x[3][1] + 0.5f, 0.0f), 255.0f);
        dst[(dy * dst_width + dx) * 3 + 2] = (uint8_t)min(max(wts_y[0] * p_x[0][2] + wts_y[1] * p_x[1][2] + wts_y[2] * p_x[2][2] + wts_y[3] * p_x[3][2] + 0.5f, 0.0f), 255.0f);
    }
}

插值算法

插值算法

线性插值

双线性插值 bilinear

双线性插值公式推导

1. 问题描述

2. 计算 $P_1$ 和 $P_2$

3. 根据 $P_1$ 和 $P_2$ 计算 $P(tx, ty)$

4. 双线性插值公式

5. 总结

用面积理解双线性插值示意图

双三次插值 bicubic

步骤 1：计算 $x$ 方向的中间插值点 $P_0, P_1, P_2, P_3$

步骤 2：计算 $y$ 方向的最终插值点 $P$

双三次插值的权重函数

权重函数中的变量说明

总结步骤

cuda 仿射变换实现

双线性插值

双三次插值

推荐阅读更多精彩内容

插值算法

线性插值

双线性插值 bilinear

双线性插值公式推导

1. 问题描述

2. 计算 和

3. 根据 和 计算

4. 双线性插值公式

5. 总结

用面积理解双线性插值示意图

双三次插值 bicubic

步骤 1：计算 方向的中间插值点

步骤 2：计算 方向的最终插值点

双三次插值的权重函数

权重函数中的变量说明

总结步骤

cuda 仿射变换实现

双线性插值

双三次插值

推荐阅读更多精彩内容

2. 计算 $P_1$ 和 $P_2$

3. 根据 $P_1$ 和 $P_2$ 计算 $P(tx, ty)$

步骤 1：计算 $x$ 方向的中间插值点 $P_0, P_1, P_2, P_3$

步骤 2：计算 $y$ 方向的最终插值点 $P$