devbox
2023面试高手
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

双目立体成像（三）视差图空洞填充_双目视差图

作者：2023面试高手 | 2024-02-16 03:28:47

赞

踩

双目视差图

在双目立体成像（二）中，对图像进行畸变矫正和立体校正之后生成的点云图效果有所改进，但是可以看出恢复出的点云图有很多洞。视差图中视差不可靠值大多数是由于遮挡或光照不均匀引起的，可以用附近可靠视差值进行填充。步骤如下：

① 计算视差图的积分图integral，并保存对应积分图中每个积分值处所有累加的像素点个数n（空洞处的像素点不计入n中，因为空洞处像素值为0，对积分值没有任何作用，反而会平滑图像）。

② 采用多层次均值滤波。首先以一个较大的初始窗口去做均值滤波，将大区域的空洞赋值。然后下次滤波时，将窗口尺寸缩小为原来的一半，利用原来的积分图再次滤波，给较小的空洞赋值（覆盖原来的值）；依次类推，直至窗口大小变为3x3，此时停止滤波，得到最终结果。

③ 多层次滤波考虑的是对于初始较大的空洞区域，需要参考更多的邻域值，如果采用较小的滤波窗口，不能够完全填充，而如果全部采用较大的窗口，则图像会被严重平滑。因此根据空洞的大小，不断调整滤波窗口。先用大窗口给所有空洞赋值，然后利用逐渐变成小窗口滤波覆盖原来的值，这样既能保证空洞能被填充上，也能保证图像不会被过度平滑。

空洞填充的函数代码如下：

  1. void insertDepth32f(cv::Mat& depth)
  2. {
  3.     const int width = depth.cols;
  4.     const int height = depth.rows;
  5.     float* data = (float*)depth.data;
  6.     cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
  7.     cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
  8.     double* integral = (double*)integralMap.data;
  9.     int* ptsIntegral = (int*)ptsMap.data;
  10.     memset(integral, 0, sizeof(double) * width * height);
  11.     memset(ptsIntegral, 0, sizeof(int) * width * height);
  12.     for (int i = 0; i < height; ++i)
  13.     {
  14.         int id1 = i * width;
  15.         for (int j = 0; j < width; ++j)
  16.         {
  17.             int id2 = id1 + j;
  18.             if (data[id2] > 1e-3)
  19.             {
  20.                 integral[id2] = data[id2];
  21.                 ptsIntegral[id2] = 1;
  22.             }
  23.         }
  24.     }
  25.     // 积分区间
  26.     for (int i = 0; i < height; ++i)
  27.     {
  28.         int id1 = i * width;
  29.         for (int j = 1; j < width; ++j)
  30.         {
  31.             int id2 = id1 + j;
  32.             integral[id2] += integral[id2 - 1];
  33.             ptsIntegral[id2] += ptsIntegral[id2 - 1];
  34.         }
  35.     }
  36.     for (int i = 1; i < height; ++i)
  37.     {
  38.         int id1 = i * width;
  39.         for (int j = 0; j < width; ++j)
  40.         {
  41.             int id2 = id1 + j;
  42.             integral[id2] += integral[id2 - width];
  43.             ptsIntegral[id2] += ptsIntegral[id2 - width];
  44.         }
  45.     }
  46.     int wnd;
  47.     double dWnd = 2;
  48.     while (dWnd > 1)
  49.     {
  50.         wnd = int(dWnd);
  51.         dWnd /= 2;
  52.         for (int i = 0; i < height; ++i)
  53.         {
  54.             int id1 = i * width;
  55.             for (int j = 0; j < width; ++j)
  56.             {
  57.                 int id2 = id1 + j;
  58.                 int left = j - wnd - 1;
  59.                 int right = j + wnd;
  60.                 int top = i - wnd - 1;
  61.                 int bot = i + wnd;
  62.                 left = max(0, left);
  63.                 right = min(right, width - 1);
  64.                 top = max(0, top);
  65.                 bot = min(bot, height - 1);
  66.                 int dx = right - left;
  67.                 int dy = (bot - top) * width;
  68.                 int idLeftTop = top * width + left;
  69.                 int idRightTop = idLeftTop + dx;
  70.                 int idLeftBot = idLeftTop + dy;
  71.                 int idRightBot = idLeftBot + dx;
  72.                 int ptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] - (ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
  73.                 double sumGray = integral[idRightBot] + integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
  74.                 if (ptsCnt <= 0)
  75.                 {
  76.                     continue;
  77.                 }
  78.                 data[id2] = float(sumGray / ptsCnt);
  79.             }
  80.         }
  81.         int s = wnd / 2 * 2 + 1;
  82.         if (s > 201)
  83.         {
  84.             s = 201;
  85.         }
  86.         cv::GaussianBlur(depth, depth, cv::Size(s, s), s, s);
  87.     }
  88. }

将左右目图片去畸变和立体校正之后计算视差，再对视差图进行空洞填补。

原图：

去畸变&立体校正后：

生成的点云图：

空洞填充后生成的点云图：

可以看出虽然物体部分恢复的细节比未填充多，但是白墙这种无纹理得区域还是空洞的。一般默认算的是左视差图，如果需要计算右视差图，则设置参数是要注意最小视差值和视差范围。由于视差值计算出来为负值，disp类型为16SC1，因此需要取绝对值，然后保存，代码如下：

  1. #include <opencv2/opencv.hpp>
  2. #include <vector>
  3. #include <string>
  4. #include <Eigen/Core>
  5. #include <pangolin/pangolin.h>
  6. #include <unistd.h>
  7. using namespace std; 
  8. using namespace Eigen;
  9.  
  10. void showPointCloud(const vector<Vector4d, aligned_allocator<Vector4d>> &pointcloud);
  11.  
  12.  
  13. void insertDepth32f(cv::Mat& depth)
  14. {
  15.     const int width = depth.cols;
  16.     const int height = depth.rows;
  17.     float* data = (float*)depth.data;
  18.     cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
  19.     cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
  20.     double* integral = (double*)integralMap.data;
  21.     int* ptsIntegral = (int*)ptsMap.data;
  22.     memset(integral, 0, sizeof(double) * width * height);
  23.     memset(ptsIntegral, 0, sizeof(int) * width * height);
  24.     for (int i = 0; i < height; ++i)
  25.     {
  26.         int id1 = i * width;
  27.         for (int j = 0; j < width; ++j)
  28.         {
  29.             int id2 = id1 + j;
  30.             if (data[id2] > 1e-3)
  31.             {
  32.                 integral[id2] = data[id2];
  33.                 ptsIntegral[id2] = 1;
  34.             }
  35.         }
  36.     }
  37.     // 积分区间
  38.     for (int i = 0; i < height; ++i)
  39.     {
  40.         int id1 = i * width;
  41.         for (int j = 1; j < width; ++j)
  42.         {
  43.             int id2 = id1 + j;
  44.             integral[id2] += integral[id2 - 1];
  45.             ptsIntegral[id2] += ptsIntegral[id2 - 1];
  46.         }
  47.     }
  48.     for (int i = 1; i < height; ++i)
  49.     {
  50.         int id1 = i * width;
  51.         for (int j = 0; j < width; ++j)
  52.         {
  53.             int id2 = id1 + j;
  54.             integral[id2] += integral[id2 - width];
  55.             ptsIntegral[id2] += ptsIntegral[id2 - width];
  56.         }
  57.     }
  58.     int wnd;
  59.     double dWnd = 2;
  60.     while (dWnd > 1)
  61.     {
  62.         wnd = int(dWnd);
  63.         dWnd /= 2;
  64.         for (int i = 0; i < height; ++i)
  65.         {
  66.             int id1 = i * width;
  67.             for (int j = 0; j < width; ++j)
  68.             {
  69.                 int id2 = id1 + j;
  70.                 int left = j - wnd - 1;
  71.                 int right = j + wnd;
  72.                 int top = i - wnd - 1;
  73.                 int bot = i + wnd;
  74.                 left = max(0, left);
  75.                 right = min(right, width - 1);
  76.                 top = max(0, top);
  77.                 bot = min(bot, height - 1);
  78.                 int dx = right - left;
  79.                 int dy = (bot - top) * width;
  80.                 int idLeftTop = top * width + left;
  81.                 int idRightTop = idLeftTop + dx;
  82.                 int idLeftBot = idLeftTop + dy;
  83.                 int idRightBot = idLeftBot + dx;
  84.                 int ptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] - (ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
  85.                 double sumGray = integral[idRightBot] + integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
  86.                 if (ptsCnt <= 0)
  87.                 {
  88.                     continue;
  89.                 }
  90.                 data[id2] = float(sumGray / ptsCnt);
  91.             }
  92.         }
  93.         int s = wnd / 2 * 2 + 1;
  94.         if (s > 201)
  95.         {
  96.             s = 201;
  97.         }
  98.         cv::GaussianBlur(depth, depth, cv::Size(s, s), s, s);
  99.     }
  100. }
  101.  
  102.  
  103. int main(int argc,char **argv)
  104. {
  105.     if(argc!=3)
  106.     {
  107.         cout<<"ERROR:Please input pictures!"<<endl;
  108.     }
  109.     cv::Mat left_image=cv::imread(argv[1],0);
  110.     cv::Mat right_image=cv::imread(argv[2],0);
  111.     // 内参
  112.     double fx = 392.7540541781677, fy = 395.9067564822596, cx = 320.5167125964447, cy = 214.0481721219656;
  113.     // 基线
  114.     double b=60;//基线单位是mm
  115.     cv::Mat disparity_sgbm,disparity;
  116.     //参数设置
  117.         enum { STEREO_BM = 0, STEREO_SGBM = 1, STEREO_HH = 2, STEREO_VAR = 3, STEREO_3WAY = 4 };
  118.     int numberOfDisparities = ((left_image.size().width / 8) + 15) & -16;
  119.     cv::Ptr<cv::StereoSGBM> sgbm = cv::StereoSGBM::create(0, 16, 3);
  120.     sgbm->setPreFilterCap(63);
  121.     int SADWindowSize = 9;
  122.     int sgbmWinSize = SADWindowSize > 0 ? SADWindowSize : 3;
  123.     sgbm->setBlockSize(sgbmWinSize);
  124.     int cn = left_image.channels();
  125.     sgbm->setP1(8 * cn*sgbmWinSize*sgbmWinSize);
  126.     sgbm->setP2(32 * cn*sgbmWinSize*sgbmWinSize);
  127.     sgbm->setMinDisparity(0);
  128.     sgbm->setNumDisparities(numberOfDisparities);
  129.     sgbm->setUniquenessRatio(10);
  130.     sgbm->setSpeckleWindowSize(100);
  131.     sgbm->setSpeckleRange(32);
  132.     sgbm->setDisp12MaxDiff(1);
  133.     int alg = STEREO_SGBM;
  134.     if (alg == STEREO_HH)
  135.         sgbm->setMode(cv::StereoSGBM::MODE_HH);
  136.     else if (alg == STEREO_SGBM)
  137.         sgbm->setMode(cv::StereoSGBM::MODE_SGBM);
  138.     else if (alg == STEREO_3WAY)
  139.         sgbm->setMode(cv::StereoSGBM::MODE_SGBM_3WAY);
  140.     //参数设置    
  141.     sgbm->compute(left_image,right_image,disparity_sgbm);
  142.     disparity_sgbm.convertTo(disparity,CV_32F,1.0/16.0f);    
  143.     // 视差图后处理
  144.     insertDepth32f(disparity);
  145.     cv::imwrite("./disparity.jpg",disparity);
  146.     vector<Vector4d,Eigen::aligned_allocator<Vector4d>> pointcloud;
  147.     for(int v=0;v<left_image.rows;v++)
  148.     {
  149.         for(int u=0;u<left_image.cols;u++)
  150.         {
  151.             if(disparity.at<float>(v,u)<=0.0||disparity.at<float>(v,u)>=96.0)
  152.             {
  153.                 continue;
  154.             }
  155.             Vector4d point(0,0,0,left_image.at<uchar>(v,u)/255.0);
  156.             double x = (u-cx)/fx;
  157.             double y = (v-cy)/fy;
  158.             double depth=fx*b/(disparity.at<float>(v,u));
  159.             // double depth=depthImg_left.at<float>(v,u);
  160.             point[0]=x*depth;
  161.             point[1]=y*depth;
  162.             point[2]=depth;
  163.  
  164.             pointcloud.push_back(point);
  165.         }
  166.     }
  167.  
  168.     // 右视差图
  169.     cv::Mat disparity_sgbm_r,disparity_r;
  170.     int numberOfDisparities_r = ((right_image.size().width / 8) + 15) & -16;
  171.     cv::Ptr<cv::StereoSGBM> sgbm_r = cv::StereoSGBM::create(0, 16, 3);
  172.     sgbm_r->setPreFilterCap(63);
  173.     sgbm_r->setBlockSize(sgbmWinSize);
  174.     int cn_r = right_image.channels();
  175.     sgbm_r->setP1(8 * cn_r*sgbmWinSize*sgbmWinSize);
  176.     sgbm_r->setP2(32 * cn_r*sgbmWinSize*sgbmWinSize);
  177.     sgbm_r->setMinDisparity(-numberOfDisparities_r);
  178.     sgbm_r->setNumDisparities(numberOfDisparities_r);
  179.     sgbm_r->setUniquenessRatio(10);
  180.     sgbm_r->setSpeckleWindowSize(100);
  181.     sgbm_r->setSpeckleRange(32);
  182.     sgbm_r->setDisp12MaxDiff(1);
  183.     if (alg == STEREO_HH)
  184.         sgbm_r->setMode(cv::StereoSGBM::MODE_HH);
  185.     else if (alg == STEREO_SGBM)
  186.         sgbm_r->setMode(cv::StereoSGBM::MODE_SGBM);
  187.     else if (alg == STEREO_3WAY)
  188.         sgbm_r->setMode(cv::StereoSGBM::MODE_SGBM_3WAY);
  189.  
  190.     sgbm_r->compute(right_image,left_image,disparity_sgbm_r);
  191.     disparity_sgbm_r.convertTo(disparity_r,CV_32F,1.0/16.0f);
  192.     disparity_r=abs(disparity_r);
  193.     insertDepth32f(disparity_r);
  194.     cv::imwrite("./disparity_r.jpg",disparity_r); 
  195.     vector<Vector4d,Eigen::aligned_allocator<Vector4d>> pointcloud_r;
  196.     for(int v=0;v<right_image.rows;v++)
  197.     {
  198.         for(int u=0;u<right_image.cols;u++)
  199.         {
  200.             if(disparity_r.at<float>(v,u)<=10.0||disparity_r.at<float>(v,u)>=96.0)
  201.             {
  202.                 continue;
  203.             }
  204.             Vector4d point1(0,0,0,right_image.at<uchar>(v,u)/255.0);
  205.             double x1 = (u-cx)/fx;
  206.             double y1 = (v-cy)/fy;
  207.             double depth1=fx*b/(disparity_r.at<float>(v,u));
  208.             // double depth1=depthImg_right.at<float>(v,u);
  209.             point1[0]=x1*depth1;
  210.             point1[1]=y1*depth1;
  211.             point1[2]=depth1;
  212.  
  213.             pointcloud_r.push_back(point1);
  214.         }
  215.     }
  216.     //右视差图
  217.  
  218.         cv::imshow("disparity",disparity/96.0);
  219.         cv::imshow("disparity_r",disparity_r/96.0);
  220.  
  221.         // cv::imshow("depth_left",depthImg_left/96.0);
  222.         // cv::imshow("depth_right",depthImg_right/96.0);
  223.  
  224.         cv::waitKey(0);
  225.         cout<<argv[3]<<endl;
  226.         //showPointCloud(pointcloud);
  227.         showPointCloud(pointcloud_r);
  228.         return 0; 
  229.  
  230. }
  231.  
  232. void showPointCloud(const vector<Vector4d, Eigen::aligned_allocator<Vector4d>> &pointcloud) 
  233. {
  234.  
  235.     if (pointcloud.empty()) {
  236.         cerr << "Point cloud is empty!" << endl;
  237.         return;
  238.     }
  239.     //创建一个pangolin窗口，窗口名“Point Cloud Viewer”,窗口宽度1024,窗口高度768
  240.     pangolin::CreateWindowAndBind("Point Cloud Viewer", 1024, 768);
  241.     glEnable(GL_DEPTH_TEST);//根据物体远近，实现遮挡效果
  242.     glEnable(GL_BLEND);//使用颜色混合模型，让物体显示半透明效果
  243.     //GL_SRC_ALPHA表示使用源颜色的alpha值作为权重因子，GL_ONE_MINUS_SRC_ALPHA表示使用(1-源颜色的alpha值)作为权重因子
  244.     glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  245.     /*
  246.     //ProjectionMatrix()中各参数依次为
  247.     图像宽度=1024
  248.     图像高度=768
  249.     fx=500、fy=500
  250.     cx=512、cy=389
  251.     最近距离=0.1和最远距离=1000
  252.     *****************************************************************
  253.     ModelViewLookAt()中各参数为相机位置，被观察点位置和相机哪个轴朝上
  254.     ModelViewLookAt(0, -0.1, -1.8, 0, 0, 0, 0.0, -1.0, 0.0)
  255.     表示相机在(0, -0.1, -1.8)位置处观看视点(0, 0, 0)，
  256.     并设置相机XYZ轴正方向为（0，-1，0），即右上前
  257.     */
  258.     pangolin::OpenGlRenderState s_cam(
  259.         pangolin::ProjectionMatrix(1024, 768, 500, 500, 512, 389, 0.1, 1000),
  260.         pangolin::ModelViewLookAt(0, -0.1, -1.8, 0, 0, 0, 0.0, -1.0, 0.0)
  261.     );
  262.     //创建交互视图，显示上一帧图像内容
  263.     pangolin::View &d_cam = pangolin::CreateDisplay()
  264.         .SetBounds(0.0, 1.0, pangolin::Attach::Pix(175), 1.0, -1024.0f / 768.0f)
  265.         .SetHandler(new pangolin::Handler3D(s_cam));
  266.     /*
  267.     SetBounds()内的前4个参数分别表示交互视图的大小，均为相对值，范围在0.0至1.0之间
  268.     第1个参数表示bottom，即为视图最下面在整个窗口中的位置
  269.     第2个参数为top，即为视图最上面在整个窗口中的位置
  270.     第3个参数为left，即视图最左边在整个窗口中的位置
  271.     第4个参数为right，即为视图最右边在整个窗口中的位置
  272.     第5个参数为aspect，表示横纵比
  273.     */
  274.     //如果pangolin窗口没有关闭，则执行
  275.     while (pangolin::ShouldQuit() == false) 
  276.     {
  277.         //清空颜色和深度缓存，使得前后帧不会互相干扰
  278.         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  279.         //激活显示，并设置相机状态
  280.         d_cam.Activate(s_cam);
  281.         //设置背景颜色为白色
  282.         glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
  283.  
  284.         glPointSize(2);
  285.         //绘制点云
  286.         glBegin(GL_POINTS);
  287.         for (auto &p: pointcloud) 
  288.         {
  289.             //设置颜色信息
  290.             glColor3f(p[3], p[3], p[3]);
  291.             //设置位置信息
  292.             glVertex3d(p[0], p[1], p[2]);
  293.         }
  294.         glEnd();
  295.         //按照上面的设置执行渲染
  296.         pangolin::FinishFrame();
  297.         //停止执行5ms
  298.         usleep(5000);   // sleep 5 ms
  299.     }
  300.     return;
  301. }

右视差图

生成的点云图：

其他场景测试

原图校正后视差图点云图
左
右
原图校正后视差图点云图
左
右
原图校正后视差图点云图
左
右

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/2023面试高手/article/detail/90365
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号