图像实现曲面屏效果

双线性插值

双线性插值是一种常用的图像插值方法，用于在图像中两个相邻像素之间进行插值，以获取介于它们之间某个位置的像素值。在透视变换等情况下，由于原始图像的像素点与目标图像的像素点位置不完全重合，因此需要对目标图像中的像素值进行插值。

以下是双线性插值的一般步骤：

1. 确定目标图像中的坐标 $(x^{\prime },y^{\prime })$，其中 $x^{\prime }$ 和 $y^{\prime }$ 可能是浮点数，不一定是整数像素坐标。
2. 将 $(x^{\prime },y^{\prime })$ 分解为四个最近的整数坐标 $(x_1,y_1)$、$(x_2,y_1)$、$(x_1,y_2)$ 和 $(x_2,y_2)$。
3. 计算目标坐标在原始图像中的相对位置 $dx$ 和 $dy$：$dx=x^{\prime }-x_1$，$dy=y^{\prime }-y_1$。
4. 分别获取这四个最近整数坐标的像素值 $f(x_1,y_1)$、$f(x_2,y_1)$、$f(x_1,y_2)$ 和 $f(x_2,y_2)$。
5. 使用双线性插值公式计算 $(x^{\prime },y^{\prime })$ 处的像素值：
   $f(x^{\prime },y^{\prime })=(1-dx)(1-dy)\cdot f(x_1,y_1)+dx(1-dy)\cdot f(x_2,y_1)+(1-dx)dy\cdot f(x_1,y_2)+dxdy\cdot f(x_2,y_2)$

通过这种方法，您可以根据目标坐标在原始图像中的位置，使用双线性插值得到相应的像素值。

代码实现：

// 定义双线性插值函数 - 通过(x,y)坐标获取像素值
    float bilinear_interpolation(const Eigen::MatrixXf& image, float x, float y)
    {
        // 获取图像的宽度和高度
        int width = image.cols();
        int height = image.rows();

        // 计算四个最近的像素的坐标
        int x0 = static_cast<int>(x);
        int y0 = static_cast<int>(y);
        //不能超出图像边界
        int x1 = std::min(x0 + 1, width - 1);
        int y1 = std::min(y0 + 1, height - 1);

        // 计算双线性插值系数
        float alpha = x - x0;
        float beta = y - y0;

        //1 3
        //2 4
        // 计算四个最近的像素的灰度值
        float f00 = image(y0, x0); //1
        float f10 = image(y1, x0); //2
        float f01 = image(y0, x1); //3
        float f11 = image(y1, x1); //4

        // 执行双线性插值
        return ((1 - alpha) * (1 - beta) * f00 +
                (1 - alpha) * beta * f10 +
                alpha * (1 - beta) * f01 +
                alpha * beta * f11);
    }

二维透视变换

二维透视变换是一种将二维图像中的点映射到另一个二维平面上的变换。通常，这种变换可以用一个3x3的矩阵表示，称为透视变换矩阵。透视变换矩阵可以根据具体的变换需求进行构造，例如平移、旋转、缩放和倾斜等。

假设我们有一个二维点 $(x,y)$，它经过透视变换后得到对应的点 $(x^{\prime },y^{\prime })$，则透视变换可以表示为以下形式：

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}\right)=\left(\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right)\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)$

其中，$(x,y)$ 是原始点的坐标，$(x^{\prime },y^{\prime })$ 是变换后的点的坐标，$(a,b,c,d,e,f,g,h,i)$ 是透视变换矩阵的参数。

透视变换矩阵的参数决定了变换的效果，例如：

• $a$ 和 $e$ 控制了水平和垂直方向的缩放；
• $b$ 和 $d$ 控制了水平和垂直方向的旋转；
• $c$ 和 $f$ 控制了水平和垂直方向的平移；
• $g$ 和 $h$ 控制了透视效果；
• $i$ 通常为1，用于保持齐次坐标的性质。

通过调整透视变换矩阵的参数，可以实现各种不同的二维透视变换，例如仿射变换、透视变换和投影变换等。

反推x

方程可以表示为 v =Mx。其中，M 是你的转换矩阵，x 是你想要求解的向量，v 是结果向量。

在你的例子中，你已经有了矩阵 M 和结果向量 (x1, y1, z1)。要求解原始的向量 (x, y, z)，你可以将方程重写为 x = M^(-1) * v，其中 M^(-1) 是 M 的逆矩阵，* 表示矩阵乘法。

以下是一个示例代码，演示了如何使用 Eigen 库求解线性方程组：

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

int main() {
    // 定义转换矩阵 M
    Matrix3f M;
    M << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;

    // 定义结果向量
    Vector3f v_result;
    v_result << 10, 11, 12;

    // 求解原始向量 x
    Vector3f x = M.inverse() * v_result;

    // 打印结果
    std::cout << "Original vector (x, y, z) = (" << x(0) << ", " << x(1) << ", " << x(2) << ")" << std::endl;

    return 0;
}

在这个示例中，我们首先定义了转换矩阵 M 和结果向量 v_result。然后，我们使用 M.inverse() * v_result 来求解原始向量 x。最后，我们打印了求解出的原始向量的值。

c++ eigen 实现

#include <Eigen/Dense>
#include <iostream>
#include <cmath>
#include <cassert>
#include <QImage>
struct Point
{
    float x = -1;
    float y = -1;
};

using MatrixP = Eigen::Matrix<Point,Eigen::Dynamic,Eigen::Dynamic>;
class HyperboloidMapping
{
public:
    //1 2
    //4 3
    HyperboloidMapping(const Eigen::Vector2f src_points[4], const Eigen::Vector2f dst_points[4],int src_height,int src_width,int dst_height = -1,int dst_width = -1)
        :_height(src_height),_width(src_width)
    {
        _transform_matrix = get_perspective_transform(src_points,dst_points);

        //找出目标边界
        _min_width = std::min(std::min(dst_points[0](0),dst_points[1](0)),std::min(dst_points[2](0),dst_points[3](0)));
        _min_height = std::min(std::min(dst_points[0](1),dst_points[1](1)),std::min(dst_points[2](1),dst_points[3](1)));

        _max_width = std::max(std::max(dst_points[0](0),dst_points[1](0)),std::max(dst_points[2](0),dst_points[3](0)));
        _max_height = std::max(std::max(dst_points[0](1),dst_points[1](1)),std::max(dst_points[2](1),dst_points[3](1)));

        if(dst_height <= 0)
        {
            dst_height = _height;
        }

        if(dst_width <= 0)
        {
            dst_width = _width;
        }
        _perspective_mapping_table = get_perspective_mapping_table(_transform_matrix,dst_height,dst_width);
    }

    //应用双曲
    /
    /// \brief warp_hyperbola
    /// \param src_image    //输入数据
    /// \param dst_image    //输出数据
    /// \param radian       //曲面弧度 限制在 [0, 1] 范围内
    /// \param crop         //裁剪有效数据到dst_image
    ///
    void warp_hyperbola(const Eigen::MatrixXf (&src_image)[3],Eigen::MatrixXf (&dst_image)[3],float radian = 0.1,bool crop = false)
    {

        assert(src_image[0].rows() == _height);
        assert(src_image[1].rows() == _height);
        assert(src_image[2].rows() == _height);

        assert(src_image[0].cols() == _width);
        assert(src_image[1].cols() == _width);
        assert(src_image[2].cols() == _width);

        Eigen::MatrixXf perspective_dst_image[3];
        //先进行透视变换
        warp_perspective(src_image,perspective_dst_image,_perspective_mapping_table);

        int src_rows = perspective_dst_image[0].rows();
        int src_cols = perspective_dst_image[0].cols();
        for(int i = 0;i < 3; ++i)
        {
            dst_image[i] = Eigen::MatrixXf::Zero(src_rows,src_cols);
        }

        //用于记录当前位置的累计像素个数 ，最后求平局像素值
        Eigen::MatrixXf accumulative_total = Eigen::MatrixXf::Zero(src_rows,src_cols);

        //坐标原点
        int center_x = src_cols >> 1;
        int center_y = src_rows >> 1;

        //radian 的值限制在 [0, 1] 范围内
        radian = std::min(1.0f, radian);
        radian = std::max(0.0f, radian);

        //内缩
        int retract = center_y * radian;

        //系数
        float coefficient = (float)retract / (center_x * center_x);
        //系数步长
        float coefficient_step = coefficient / center_y;

        // 计算映射表 y位置
        for (int y = 0; y < center_y; ++y)
        {
            //变量x步长
            float x_step = y * coefficient_step;
            //变量常数大小
            float y_step = (center_y - y) - ((center_x) * (center_x) * (coefficient - x_step));


            //获取映射y坐标
            for (int x = center_x; x < src_cols; ++x)
            {
                //获取第一象限y位置
                int y1 = center_y - static_cast<float>((x - center_x) * (x - center_x) * (coefficient - x_step) + (y_step));
                accumulative_total(y1, x)++;

                //获取第二象限y位置
                int x2 = center_x - (x - center_x);
                int y2 = y1;
                accumulative_total(y2, x2)++;

                //获取第三象限y位置
                int x3 = x2;
                int y3 = src_rows - y1 - 1;
                accumulative_total(y3, x3)++;

                //获取第四象限y位置
                int x4 = x;
                int y4 = y3;
                accumulative_total(y4, x4)++;

                for(int i = 0;i < 3; ++i)
                {
                    //获取第一象限y位置
                    dst_image[i](y1, x) += perspective_dst_image[i](y,x);

                    //获取第二象限y位置
                    dst_image[i](y2, x2) += perspective_dst_image[i](y,x2);

                    //获取第三象限y位置
                    dst_image[i](y3, x3) += perspective_dst_image[i](src_rows - y - 1,x3);

                    //获取第四象限y位置
                    dst_image[i](y4, x4) += perspective_dst_image[i](src_rows - y - 1,x4);
                }
            }
        }

        for(int i = 0;i < 3; ++i)
        {
            // 对位相除
            dst_image[i] = dst_image[i].cwiseQuotient(accumulative_total);
        }

        //裁剪
        if(crop)
        {
            int min_width = _min_width > src_cols ? 0 : _min_width;
            int max_width = _max_width > src_cols ? src_cols : _max_width;
            int min_height = _min_height > src_rows ? 0 : _min_height;
            int max_height = _max_height > src_rows ? src_rows : _max_height;
            for(int i = 0;i < 3; ++i)
            {
                Eigen::MatrixXf temp = dst_image[i].block(min_height, min_width, max_height-min_height, max_width-min_width);
                dst_image[i] = temp;
            }
        }
    }

    //qimage 转 Matrix
    static void image_2_matrix(const QImage& image,Eigen::MatrixXf (&matrix)[3])
    {
        // 将图像转换为 RGB888 格式
        QImage convertedImage = image.convertToFormat(QImage::Format_RGB888);

        for(int i=0;i<3;++i)
        {
            matrix[i] = Eigen::MatrixXf(convertedImage.height(), convertedImage.width());
        }

        // 从 QImage 中提取每个通道的数据
        for (int i = 0; i < convertedImage.height(); ++i)
        {
            for (int j = 0; j < convertedImage.width(); ++j)
            {
                // 获取像素的 RGB 值
                QRgb pixel = convertedImage.pixel(j, i);

                // 将 RGB 值拆分为每个通道的值，并保存到对应的 Eigen Matrix 中
                matrix[0](i, j) = qRed(pixel);
                matrix[1](i, j) = qGreen(pixel);
                matrix[2](i, j) = qBlue(pixel);
            }
        }

    }

    // 将三个 Eigen Matrix 转换为 QImage
    static void matrix_2_image(const Eigen::MatrixXf (&matrix)[3],QImage& image)
    {
        // 创建一个空的 QImage
        image = QImage(matrix[0].cols(), matrix[0].rows(), QImage::Format_RGB888);

        // 设置图像的像素值
        for (int i = 0; i < image.height(); ++i)
        {
            for (int j = 0; j < image.width(); ++j)
            {
                // 创建 QColor 对象并设置像素值
                QColor color(static_cast<int>(matrix[0](i, j)), static_cast<int>(matrix[1](i, j)), static_cast<int>(matrix[2](i, j)));
                image.setPixel(j, i, color.rgb());
            }
        }
    }
protected:
    //src_points
    //1 2
    //4 3
    //获取透视变化矩阵
    Eigen::Matrix3f get_perspective_transform(const Eigen::Vector2f src_points[4], const Eigen::Vector2f dst_points[4])
    {
        Eigen::Matrix3f perspective_matrix;

        // 构造线性方程组
        Eigen::Matrix<float, 8, 8> A;
        Eigen::Matrix<float, 8, 1> b;
        for (int i = 0; i < 4; ++i) {
            A.row(i * 2) << src_points[i].x(), src_points[i].y(), 1, 0, 0, 0, -dst_points[i].x() * src_points[i].x(), -dst_points[i].x() * src_points[i].y();
            A.row(i * 2 + 1) << 0, 0, 0, src_points[i].x(), src_points[i].y(), 1, -dst_points[i].y() * src_points[i].x(), -dst_points[i].y() * src_points[i].y();
            b.row(i * 2) << dst_points[i].x();
            b.row(i * 2 + 1) << dst_points[i].y();
        }

        // 解线性方程组 
         /* 在 Eigen 库中，`colPivHouseholderQr()` 是用于执行列主元素高斯-约当消元法的方法，
         * 用于解线性方程组。它返回一个对象，该对象提供了一种求解线性方程组的方式。
         * 在你的代码中，`A.colPivHouseholderQr().solve(b)` 表示对矩阵 `A` 应用列主元素高斯-约当消元法，
         * 并解出线性方程组 `Ax = b`，其中 `b` 是右侧的常数向量，`x` 是未知向量。解出的向量 `x` 包含了方程组的解。
         */
        Eigen::Matrix<float, 8, 1> x = A.colPivHouseholderQr().solve(b);

        // 构造透视变换矩阵
        perspective_matrix << x[0], x[1], x[2],
            x[3], x[4], x[5],
            x[6], x[7], 1;

        return perspective_matrix;
    }

    // 定义双线性插值函数 - 通过(x,y)坐标获取像素值
    float bilinear_interpolation(const Eigen::MatrixXf& image, float x, float y)
    {
        // 获取图像的宽度和高度
        int width = image.cols();
        int height = image.rows();

        // 计算四个最近的像素的坐标
        int x0 = static_cast<int>(x);
        int y0 = static_cast<int>(y);
        //不能超出图像边界
        int x1 = std::min(x0 + 1, width - 1);
        int y1 = std::min(y0 + 1, height - 1);

        // 计算双线性插值系数
        float alpha = x - x0;
        float beta = y - y0;

        //1 3
        //2 4
        // 计算四个最近的像素的灰度值
        float f00 = image(y0, x0); //1
        float f10 = image(y1, x0); //2
        float f01 = image(y0, x1); //3
        float f11 = image(y1, x1); //4

        // 执行双线性插值
        return ((1 - alpha) * (1 - beta) * f00 +
                (1 - alpha) * beta * f10 +
                alpha * (1 - beta) * f01 +
                alpha * beta * f11);
    }


    //透视变换映射表
    MatrixP get_perspective_mapping_table(const Eigen::Matrix3f& transform,int dst_height,int dst_width)
    {
        //找出目标边界
        int dst_min_width = _min_width;
        int dst_min_height = _min_height;

        int dst_max_width = _max_width;
        int dst_max_height = _max_height;


        //求逆
        Eigen::Matrix3f transform_inv = transform.inverse();

        //映射表
        MatrixP matp(dst_height,dst_width);

        // 遍历目标图像的每个像素，并进行透视变换
        for (int y = dst_min_height; y < dst_max_height; ++y)
        {
            for (int x = dst_min_width; x < dst_max_width; ++x)
            {
                //在规定的范围呢
                if(y < dst_height && x < dst_width)
                {
                    // 构建齐次坐标向量
                    Eigen::Vector3f src_point(x, y, 1);

                    // 应用透视变换
                    Eigen::Vector3f dst_point = transform_inv * src_point;

                    // 归一化坐标
                    float u = dst_point[0] / dst_point[2];
                    float v = dst_point[1] / dst_point[2];


                    //构建映射表
                    matp(y,x) = {u,v};
                }

            }
        }
        return matp;
    }

    //应用透视
    void warp_perspective(const Eigen::MatrixXf (&src_image)[3],Eigen::MatrixXf (&dst_image)[3],const MatrixP& mapping_table)
    {
        int src_rows = src_image[0].rows();
        int src_cols = src_image[0].cols();
        for(int i = 0;i < 3; ++i)
        {
            dst_image[i] = Eigen::MatrixXf::Zero(mapping_table.rows(),mapping_table.cols());
        }
        // 遍历目标图像的每个像素，并进行透视变换
        for (int y = 0; y < mapping_table.rows(); ++y)
        {
            for (int x = 0; x < mapping_table.cols(); ++x)
            {
                Point point = mapping_table(y,x);
                // 对坐标进行边界检查s
                if (point.x >= 0 && point.x < src_cols && point.y >= 0 && point.y < src_rows)
                {
                    for(int i = 0;i < 3; ++i)
                    {
                        dst_image[i](y,x) = bilinear_interpolation(src_image[i],point.x,point.y);
                    }
                }
            }
        }

    }

private:
    //输入图大小
    int _height;
    int _width;

    //有效数据范围
    int _min_width;
    int _max_width;
    int _min_height;
    int _max_height;

    Eigen::Matrix3f _transform_matrix; //透视变换矩阵
    MatrixP _perspective_mapping_table; //透视映射表
};

 // 读取图像
    QImage image("G:/Snipaste_2024-02-28_20-06-36.jpg");
    Eigen::MatrixXf src_image[3];
    HyperboloidMapping::image_2_matrix(image,src_image);


    // 定义原始图像的四个角点和目标图像的四个角点
    Eigen::Vector2f src_points[] = {
        Eigen::Vector2f(0, 0),
        Eigen::Vector2f(src_image[0].cols() - 1, 0),
        Eigen::Vector2f(src_image[0].cols() - 1, src_image[0].rows() - 1),
        Eigen::Vector2f(0, src_image[0].rows() - 1)
    };
    Eigen::Vector2f dst_points[] = {
        Eigen::Vector2f(0, 0),
        Eigen::Vector2f(src_image[0].cols() * 0.2f,0),
        Eigen::Vector2f(src_image[0].cols() * 0.2f, src_image[0].rows()),
        Eigen::Vector2f(0, src_image[0].rows() - 1)
    };




    HyperboloidMapping hy(src_points, dst_points,src_image[0].rows(),src_image[0].cols());
    Eigen::MatrixXf dst_image[3];

    // 开始计时
    auto start = std::chrono::high_resolution_clock::now();
    hy.warp_hyperbola(src_image,dst_image,0.5,true);
    // 结束计时
    auto end = std::chrono::high_resolution_clock::now();
    // 计算持续时间
    auto duration = end - start;

    auto durationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(duration);
    std::cout << "::" << durationInMilliseconds.count() << "ms" << std::endl;

    //转换image
    QImage sss;
    HyperboloidMapping::matrix_2_image(dst_image,sss);
    

OpenCV实现

#include<iostream>
#include <opencv2/opencv.hpp>
#include <cmath>
int main() {
    // 读取图像
    cv::Mat image = cv::imread("G:/Snipaste_2024-02-28_20-06-36.jpg");

    cv::resize(image, image, cv::Size(2000, 1200));
    cv::Mat result;

    // 定义原始图像的四个角点和目标图像的四个角点
    std::vector<cv::Point2f> srcPoints = {
        cv::Point2f(0, 0),
        cv::Point2f(image.cols - 1, 0),
        cv::Point2f(image.cols - 1, image.rows - 1),
        cv::Point2f(0, image.rows - 1)
    };
    std::vector<cv::Point2f> dstPoints = {
        cv::Point2f(100, 0),
        cv::Point2f(image.cols * 0.5f - 1,0),
        cv::Point2f(image.cols * 0.5f - 1, image.rows - 1),
        cv::Point2f(100, image.rows - 1)
    };

    // 计算透视变换矩阵
    cv::Mat perspectiveMatrix = cv::getPerspectiveTransform(srcPoints, dstPoints);

    std::cout << perspectiveMatrix << std::endl;
    auto start = std::chrono::high_resolution_clock::now();
    // 应用变换
    cv::warpPerspective(image, result, perspectiveMatrix, cv::Size(image.cols, image.rows));
    // 结束计时
    auto end = std::chrono::high_resolution_clock::now();
    // 计算持续时间
    auto duration = end - start;

    auto durationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(duration);
    std::cout << "::" << durationInMilliseconds.count() << "ms" << std::endl;
    



    // 创建映射表
    cv::Mat map_x(image.size(), CV_32FC1);
    cv::Mat map_y(image.size(), CV_32FC1);
    //初始化，建立映射表

    int center_x = (image.cols) >> 1;
    int center_y = (image.rows) >> 1;


    //内缩
    int retract = 100;

    //系数
    float coefficient = (float)retract / (center_x * center_x);
    //系数步长
    float coefficient_step = coefficient / center_y;

    // 计算映射表 x位置保持不变
    for (int y = 0; y < image.rows; ++y)
    {
        for (int x = 0; x < image.cols; ++x)
        {
            map_x.at<float>(y, x) = x;
            map_y.at<float>(y, x) = -1;
        }
    }
    // 计算映射表 y位置
    for (int y = 0; y < center_y; ++y)
    {
        //变量x步长
        float x_step = y * coefficient_step;
        //变量常数大小
        float y_step = (center_y - y) - ((center_x) * (center_x) * (coefficient - x_step));


        //获取映射y坐标
        for (int x = center_x; x < image.cols; ++x)
        {

            //获取第一象限y位置
            int y1 = center_y - static_cast<double>((x - center_x) * (x - center_x) * (coefficient - x_step) + (y_step));
            map_y.at<float>(y1, x) = y;

            //获取第二象限y位置
            int x2 = center_x - (x - center_x);
            int y2 = y1;
            map_y.at<float>(y2, x2) = y;

            //获取第三象限y位置
            double x3 = x2;
            double y3 = image.rows - y1 - 1;
            map_y.at<float>(y3, x3) = image.rows - y - 1;

            //获取第四象限y位置
            double x4 = x;
            double y4 = y3;
            map_y.at<float>(y4, x4) = image.rows - y - 1;

        }
    }
    {
        // 开始计时
        auto start = std::chrono::high_resolution_clock::now();
        // 应用映射表
        cv::Mat mapped_image;
        cv::remap(result, mapped_image, map_x, map_y, cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0, 0, 0));
        // 结束计时
        auto end = std::chrono::high_resolution_clock::now();
        // 计算持续时间
        auto duration = end - start;

        auto durationInMilliseconds = std::chrono::duration_cast<std::chrono::milliseconds>(duration);
        std::cout << "::" << durationInMilliseconds.count() << "ms" << std::endl;
        cv::imwrite("G:/mapped_image.jpg", mapped_image);
        // 显示结果
        cv::imshow("Original Image", image);
        cv::imshow("Mapped Image", mapped_image);
        cv::waitKey(0);

    }
    return 0;
}

对比 OpenCV实现速度 10ms 而c++和eigen实现100ms，快了10倍。

实现效果

原理：使用一元二次函数

这种效果就如同圆柱体上贴纸形状，图像围绕圆柱体3D效果。
想到了双曲线函数，最后发现一元二次函数可以替代， 实现比较方便。

数学离散描述地址


a = 表示内缩
w = 图像宽度的一半
h = 图像高度的一半

按照上面的公式，把像素映射到函数的位置，注意：（中间过程会重叠像素，当前位置累加像素，最后除以累加次数，求平均值。）

有更好的实现方法，欢迎讨论。