图像实现曲面屏效果

双线性插值

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

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

确定目标图像中的坐标 $(x^{'}, y^{'})$ ，其中 $x^{'}$ 和 $y^{'}$ 可能是浮点数，不一定是整数像素坐标。
将 $(x^{'}, y^{'})$ 分解为四个最近的整数坐标 $x_1, y_1)$ 、 $x_2, y_1)$ 、 $x_1, y_2)$ 和 $x_2, y_2)$ 。
计算目标坐标在原始图像中的相对位置 $d x$ 和 $d y$ ： $dx = x’ – x_1$ ， $dy = y’ – y_1$ 。
分别获取这四个最近整数坐标的像素值 $f(x_1, y_1)$ 、 $f(x_2, y_1)$ 、 $f(x_1, y_2)$ 和 $f(x_2, y_2)$ 。
使用双线性插值公式计算 $(x^{'}, y^{'})$ 处的像素值：
$\cdot f(x_1, y_1) + dx(1 – dy) \cdot f(x_2, y_1) + (1 – dx)dy \cdot f(x_1, y_2) + dx dy \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); //1float f10 = image(y1, x0); //2float f01 = image(y0, x1); //3float f11 = image(y1, x1); //4// 执行双线性插值return ((1 - alpha) * (1 - beta) * f00 +(1 - alpha) * beta * f10 +alpha * (1 - beta) * f01 +alpha * beta * f11);}

二维透视变换

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

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

$\begin{pmatrix} x’ \\ y’ \\ w’ \end{pmatrix} = \begin{pmatrix} a & b & c \\ d & e & f \\ g & h & i \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix}$

其中， $(x, y)$ 是原始点的坐标， $(x^{'}, y^{'})$ 是变换后的点的坐标， $(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 #include using namespace Eigen;int main() {// 定义转换矩阵 MMatrix3f M;M << 1, 2, 3, 4, 5, 6, 7, 8, 9;// 定义结果向量Vector3f v_result;v_result << 10, 11, 12;// 求解原始向量 xVector3f 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 #include #include #include #include struct Point{float x = -1;float y = -1;};using MatrixP = Eigen::Matrix<Point,Eigen::Dynamic,Eigen::Dynamic>;class HyperboloidMapping{public://1 2//4 3HyperboloidMapping(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 转 Matrixstatic 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 转换为 QImagestatic void matrix_2_image(const Eigen::MatrixXf (&matrix)[3],QImage& image){// 创建一个空的 QImageimage = 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); //1float f10 = image(y1, x0); //2float f01 = image(y0, x1); //3float 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);// 对坐标进行边界检查sif (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;//转换imageQImage sss;HyperboloidMapping::matrix_2_image(dst_image,sss);

OpenCV实现

#include#include #include 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 = 图像高度的一半

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

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

图像实现曲面屏效果