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

双目视觉（五）立体匹配算法之动态规划全局匹配_全局匹配算法

作者：繁依Fanyi0 | 2024-03-23 23:03:09

赞

踩

全局匹配算法

系列文章：

  1. 双目视觉（一）双目视觉系统
  2. 双目视觉（二）双目匹配的困难和评判标准
  3. 双目视觉（三）立体匹配算法
  4. 双目视觉（四）匹配代价
  5. 双目视觉（五）立体匹配算法之动态规划全局匹配
  6. 双目视觉（六）U-V视差
  7. 【项目实战】利用Ｕ-V视差进行地面检测
  8. 【项目实践】U-V视差路面检测之动态规划　　　　　　

代码：GitHub stereo matchhttps://github.com/Xke1718He/stereo-matchicon-default.png?t=N3I4https://github.com/Xke1718He/stereo-match

注：分享的代码中，双目匹配的代价只有一种，如果觉得不错，可以给予支持，谢谢^0^

全局匹配算法：

全局匹配主要是利用图像的全局约束信息，通过构建全局能量函数，然后通过优化方法最小化全局能量函数以求得致密视差图。目前优化方法主要有：动态规划（ＤＰ）、置信传播（ＢＰ）、模拟退火、图割法（ＧＣ）等。

动态规划：

动态规划的思想就是把求解整个过程分解为若干子过程，逐个求解子过程。例如斐波那契数列，fib(6)分解为fib(5)和fib(4),而fib(5)分解为fib(4)和fib(3),fib(4)为一个重复子问题，而动态规划就是解决这类问题的方式之一。

这是我学习动态规划相关知识的链接：

动态规划 (第1讲)_哔哩哔哩_bilibili

【动态规划专题班】ACM总冠军、清华+斯坦福大神带你入门动态规划算法_哔哩哔哩_bilibili

在双目匹配中，动态规划立体匹配是基于极线约束，通过依次寻找每条极线上匹配点对的最小代价路径的动态寻优方法求解全局能量最小化，得到匹配视差图。

步骤大致如下：

  • step1.构建全局能量函数

其中，为数据项，一般表示为代价函数；为视差平滑约束项。

  • step2. 代价函数构建

代价函数构建这里采用的是块匹配，也就是以第一图像待匹配点为中心像素创建一个n*n的窗口，在第二幅图像中，沿着极线取出与基准点领域同样大小为n*n的领域，进行领域内的相似度函数的计算。

说明：代价函数采用的是Ｃ＋＋中的继承构建

  • imageType：判断左右图像的信息相似性。比如大小是否一样大呀。
  • aggregate：这也就是所谓的代价聚合函数，计算左右图像领域内的相似性度量的。
  1. class CostFunction {
  2. public:
  3.     float lambda;
  4.  
  5.     cv::Mat left;//左图像
  6.     cv::Mat right;//右图像
  7.  
  8.     int blocksize;//block大小
  9.     int margin;//块的边界
  10.  
  11.     float normCost;
  12.     float normWin;
  13.     //构造函数
  14.     CostFunction( cv::Mat left, cv::Mat right, float lambda)
  15.     {
  16. 		lambda = 1.0;
  17.         this->left = left;
  18.         this->right = right;
  19.         imageType(left, right);
  20.         this->lambda = lambda;
  21.     }
  22.     //虚构函数
  23.     virtual bool imageType(cv::Mat left, cv::Mat right) {
  24.         assert(left.size() == right.size());
  25.  
  26.         return true;
  27.     }
  28.     //代价聚合
  29.     virtual float aggregate(int x1, int x2, int y) = 0;
  30.  
  31.     float p(float cost) {
  32.         return 1 - exp(-cost / lambda);
  33.     }
  34.  
  35.     ~CostFunction() {}
  36. };
  37.  
  38. class RGBCost : public CostFunction {
  39. public:
  40.     RGBCost(cv::Mat left, cv::Mat right, float lambda) : CostFunction(left, right, lambda) {}
  41.  
  42.     bool imageType(cv::Mat left, cv::Mat right) {
  43.         assert(left.type() == right.type() && "imgL imgR types not equal");
  44.         assert(left.type() == CV_8UC3 && "img type not supported");
  45.  
  46.         return true;
  47.     }
  48.  
  49.     // aggregate over a ROI of input images
  50.     float aggregate(int x1, int x2, int y) {
  51.         float sum = 0;
  52.         for (int i = y - margin; i <= y + margin; ++i) {
  53.             cv::Vec3b* lptr = left.ptr<cv::Vec3b>(i);
  54.             cv::Vec3b* rptr = right.ptr<cv::Vec3b>(i);
  55.  
  56.             for ( int j = -margin; j <= margin; ++j) {
  57.                 sum += eukl(lptr[x1 + j], rptr[x2 + j]);      // cost function
  58.             }
  59.         }
  60.  
  61.         return sum / sqrt(255*255 + 255*255+255*255);      // normalize to winsize*1.0
  62.     }
  63.  
  64.     float eukl(cv::Vec3b l, cv::Vec3b r) {
  65.         float a = l[0] - r[0];
  66.         float b = l[1] - r[1];
  67.         float c = l[2] - r[2];
  68.         return std::sqrt(a*a + b*b + c*c);
  69.     }
  70.  
  71.     ~RGBCost() {}
  72. };
  73.  
  74. class FloatCost : public CostFunction {
  75. public:
  76.     FloatCost(cv::Mat left, cv::Mat right, float lambda) : CostFunction(left, right, lambda) {}
  77.  
  78.     bool imageType(cv::Mat left, cv::Mat right) {
  79.         assert(left.type() == right.type() && "imgL imgR types not equal");
  80.         assert(left.type() == CV_32F && "img type not supported");
  81.  
  82.         return true;
  83.     }
  84.  
  85.     // aggregate over a ROI of input images
  86.     float aggregate(int x1, int x2, int y) {
  87.         float sum = 0;
  88.         for (int i = y - margin; i <= y + margin; ++i) {
  89.             float* lptr = left.ptr<float>(i);
  90.             float* rptr = right.ptr<float>(i);
  91.  
  92.             for ( int j = -margin ; j <= margin; ++j) {
  93.                 sum += abs(lptr[x1 + j] - rptr[x2 + j]);      // cost function
  94.             }
  95.         }
  96.  
  97.         return sum / (blocksize*blocksize*lambda);
  98.     }
  99. };
  100.  
  101. class CondHistCost : public CostFunction {
  102. public:
  103.     cv::Mat nuLeft, nuRight;
  104.     CondHistCost(cv::Mat left, cv::Mat right, float lambda) : CostFunction(left, right, lambda) {
  105.         cv::Mat histl = condHist(left, 3);
  106.         nuLeft = switchColors(left, histl);
  107.         cv::Mat histr = condHist(right, 3);
  108.         nuRight = switchColors(right, histr);
  109.     }
  110.  
  111.     bool imageType(cv::Mat left, cv::Mat right) {
  112.         assert(left.type() == right.type() && "imgL imgR types not equal");
  113.         assert(left.type() == CV_32F && "img type not supported");
  114.  
  115.         return true;
  116.     }
  117.  
  118.     // aggregate over a ROI of input images
  119.     float aggregate(int x1, int x2, int y) {
  120.         float sum = 0;
  121.         for (int i = y - margin; i <= y + margin; ++i) {
  122.             float* lptr = nuLeft.ptr<float>(i);
  123.             float* rptr = nuRight.ptr<float>(i);
  124.  
  125.             for ( int j = -margin ; j <= margin; ++j) {
  126.                 sum += abs(lptr[x1 + j] - rptr[x2 + j]);      // cost function
  127.             }
  128.         }
  129.  
  130.         return sum / (blocksize*blocksize*lambda);
  131.     }
  132. };
  133.  
  134.  
  135. class GrayCost : public CostFunction {
  136. public:
  137.     GrayCost(cv::Mat left, cv::Mat right, float lambda) : CostFunction(left, right, lambda) {}
  138.  
  139.     bool imageType(cv::Mat left, cv::Mat right) {
  140.         assert(left.type() == right.type() && "imgL imgR types not equal");
  141.         assert(left.type() == CV_8UC1 && "img type not supported");
  142.  
  143.         return true;
  144.     }
  145.  
  146.     // aggregate over a ROI of input images
  147.     float aggregate(int x1, int x2, int y) {
  148.         float sum = 0;
  149.         for (int i = y - margin; i <= y + margin; ++i) {
  150.             uchar* lptr = left.ptr<uchar>(i);
  151.             uchar* rptr = right.ptr<uchar>(i);
  152.  
  153.             for ( int j = -margin; j <= margin; ++j) {
  154.                 sum += abs(lptr[x1 + j] - rptr[x2 + j]);      // cost function
  155.             }
  156.         }
  157.  
  158.         return sum / (blocksize*blocksize*255.0);
  159.     }
  160. };
  161.  
  162. class GradientCost : public CostFunction {
  163. public:
  164.     cv::Mat l_grad;   // 3 channel float
  165.     cv::Mat r_grad;   // 3 channel float
  166.  
  167.     GradientCost(const cv::Mat left, const cv::Mat right, float lambda) : CostFunction(left, right, lambda) {
  168.         l_grad = getRGBGradientAngle(left);
  169.         r_grad = getRGBGradientAngle(right);
  170.  
  171.         //displayGradientPic(l_grad);
  172.         //displayGradientPic(r_grad);
  173.     }
  174.  
  175.     bool imageType(cv::Mat left, cv::Mat right) {
  176.         assert(left.type() == right.type() && "imgL imgR types not equal");
  177.         assert(left.type() == CV_8UC3 && "img type not supported");
  178.  
  179.         return true;
  180.     }
  181.  
  182.     // aggregate over a ROI of input images
  183.     float aggregate(int x1, int x2, int y) {
  184.         float sum = 0;
  185.         for (int i = y - margin; i <= y + margin; ++i) {
  186.             cv::Vec3f* lptr = l_grad.ptr<cv::Vec3f>(i);
  187.             cv::Vec3f* rptr = r_grad.ptr<cv::Vec3f>(i);
  188.  
  189.             for ( int j = -margin; j <= margin; ++j) {
  190.                 sum += eukl(lptr[x1 + j], rptr[x2 + j]);      // cost function
  191.             }
  192.         }
  193.  
  194.         return sum / sqrt(255*255 + 255*255 + 255*255);      // normalize to winSize * 1.0
  195.     }
  196.  
  197.     float eukl(cv::Vec3f l, cv::Vec3f r) {
  198.         float a = l[0] - r[0];
  199.         float b = l[1] - r[1];
  200.         float c = l[2] - r[2];
  201.         return std::sqrt(a*a + b*b + c*c);
  202.     }
  203.  
  204.     ~GradientCost() {
  205.         l_grad.release();
  206.         r_grad.release();
  207.     }
  208. };
  209.  
  210. class CensusCost : public CostFunction {
  211. public:
  212.     int censusWindow;
  213.     int censusMargin;
  214.     CensusCost(cv::Mat left, cv::Mat right, int censusWindow, float lambda) : CostFunction(left, right, lambda) {
  215.         // census.... nimmt einen Block
  216.         this->censusWindow = censusWindow;
  217.         this->censusMargin = censusWindow / 2;
  218.  
  219.         this->normWin = censusWindow * censusWindow;
  220.         // nimmt einen Block
  221.     }
  222.  
  223.     bool imageType(cv::Mat left, cv::Mat right) {
  224.         assert(left.type() == right.type() && "imgL imgR types not equal");
  225.         assert(left.type() == CV_8UC1 && "img type not supported");
  226.  
  227.         return true;
  228.     }
  229.  
  230.     unsigned int census(int x1, int x2, int y, uchar c1, uchar c2) {
  231.         unsigned int diff = 0;
  232.  
  233.         for(int i = y - censusMargin; i <= y + censusMargin; ++i) {
  234.             uchar* lptr = left.ptr<uchar>(i);
  235.             uchar* rptr = right.ptr<uchar>(i);
  236.  
  237.             for(int j = -censusMargin; j <= censusMargin; ++j) {
  238.                 bool t1 = (c1 < lptr[x1 + j]);
  239.                 bool t2 = (c2 < rptr[x2 + j]);
  240.  
  241.                 if(t1 != t2) diff++;
  242.             }
  243.         }
  244.  
  245.         return diff; /// (censusWindow*censusWindow);
  246.     }
  247.  
  248.     float aggregate(int x1, int x2, int y) {
  249.         float sum = 0;
  250.         /*for(int i = y - margin; i <= y + margin; ++i) {
  251.             uchar *lptr = left.ptr<uchar>(i);
  252.             uchar *rptr = right.ptr<uchar>(i);
  254.             for(int j = -margin; j <= margin; ++j)
  255.                 sum += census(x1 + j, x2 + j, i, lptr[x1 + j], rptr[x2 + j]);
  256.         }*/
  257.         uchar *lptr = left.ptr<uchar>(y);
  258.         uchar *rptr = right.ptr<uchar>(y);
  259.         sum = census(x1, x2, y, lptr[x1], rptr[x2]);
  260.         return sum / normWin;
  261.     }
  262. };
  263.  
  264. class CensusFloatCost : public CostFunction {
  265. public:
  266.     int censusWindow;
  267.     int censusMargin;
  268.  
  269.     CensusFloatCost(cv::Mat left, cv::Mat right, int censusWindow, float lambda) : CostFunction(left, right, lambda) {
  270.         // census.... nimmt einen Block
  271.         this->censusWindow = censusWindow;
  272.         this->censusMargin = censusWindow / 2;
  273.     }
  274.  
  275.     bool imageType(cv::Mat left, cv::Mat right) {
  276.         assert(left.type() == right.type() && "imgL imgR types not equal");
  277.         assert(left.type() == CV_32F && "img type not supported");
  278.  
  279.         return true;
  280.     }
  281.  
  282.     unsigned int census(int x1, int x2, int y, float c1, float c2) {
  283.         unsigned int diff = 0;
  284.  
  285.         for(int i = y - censusMargin; i <= y + censusMargin; ++i) {
  286.             float* lptr = left.ptr<float>(i);
  287.             float* rptr = right.ptr<float>(i);
  288.  
  289.             for(int j = -censusMargin; j <= censusMargin; ++j) {
  290.                 bool t1 = (c1 < lptr[x1 + j]);
  291.                 bool t2 = (c2 < rptr[x2 + j]);
  292.  
  293.                 if(t1 != t2) diff++;
  294.             }
  295.         }
  296.  
  297.         return diff;
  298.     }
  299.  
  300.     float aggregate(int x1, int x2, int y) {
  301.         float sum = 0;
  302.         for(int i = y - margin; i <= y + margin; ++i) {
  303.             float *lptr = left.ptr<float>(i);
  304.             float *rptr = right.ptr<float>(i);
  305.  
  306.             for(int j = -margin; j <= margin; ++j)
  307.                 sum += census(x1 + j, x2 + j, i, lptr[x1 + j], rptr[x2 + j]);
  308.         }
  309.         float *lptr = left.ptr<float>(y);
  310.         float *rptr = right.ptr<float>(y);
  311.         //sum = census(x1, x2, y, lptr[x1], rptr[x2]);
  312.         return sum / (censusWindow*censusWindow*lambda);
  313.     }
  314. };
  315.  
  316. class RGBCensusCost : public CostFunction {
  317. public:
  318.     int censusWindow;
  319.     int censusMargin;
  320.  
  321.     RGBCensusCost(cv::Mat left, cv::Mat right, int censusWindow, float lambda) : CostFunction(left, right, lambda) {
  322.         // census.... nimmt einen Block
  323.         this->censusWindow = censusWindow;
  324.         this->censusMargin = censusWindow / 2;
  325.         normCost = censusWindow*censusWindow*3;
  326.         // nimmt einen Block
  327.     }
  328.  
  329.     bool imageType(cv::Mat left, cv::Mat right) {
  330.         assert(left.type() == right.type() && "imgL imgR types not equal");
  331.         assert(left.type() == CV_8UC3 && "img type not supported");
  332.  
  333.         return true;
  334.     }
  335.  
  336.     unsigned int census(int x1, int x2, int y, cv::Vec3b c1, cv::Vec3b c2) {
  337.         unsigned int diff = 0;
  338.  
  339.         for(int i = y - censusMargin; i <= y + censusMargin; ++i) {
  340.             cv::Vec3b* lptr = left.ptr<cv::Vec3b>(i);
  341.             cv::Vec3b* rptr = right.ptr<cv::Vec3b>(i);
  342.  
  343.             for(int j = -censusMargin; j <= censusMargin; ++j) {
  344.                 cv::Vec3b cl = lptr[x1 + j];
  345.                 cv::Vec3b cr = rptr[x2 + j];
  346.  
  347.                 for(int ch = 0; ch < 3; ++ch) {
  348.                     bool t1 = (c1[ch] < cl[ch]);
  349.                     bool t2 = (c2[ch] < cr[ch]);
  350.  
  351.                     if(t1 != t2) diff++;
  352.                 }
  353.             }
  354.         }
  355.  
  356.         return diff;
  357.     }
  358.  
  359.     float aggregate(int x1, int x2, int y) {
  360.         float sum = 0;
  361.         for(int i = y - margin; i <= y + margin; ++i) {
  362.             cv::Vec3b *lptr = left.ptr<cv::Vec3b>(i);
  363.             cv::Vec3b *rptr = right.ptr<cv::Vec3b>(i);
  364.  
  365.             for(int j = -margin; j <= margin; ++j)
  366.                 sum += census(x1 + j, x2 + j, i, lptr[x1 + j], rptr[x2 + j]);
  367.         }
  368.         //cv::Vec3b *lptr = left.ptr<cv::Vec3b>(y);
  369.         //cv::Vec3b *rptr = right.ptr<cv::Vec3b>(y);
  370.  
  371.         return sum / normCost;
  372.     }
  373. };
  374.  
  375. class RGBGradCensusCost : public CostFunction {
  376. public:
  377.     int censusWindow;
  378.     int censusMargin;
  379.  
  380.     float normCost;
  381.     float normWin;
  382.  
  383.     cv::Mat l_grad;
  384.     cv::Mat r_grad;
  385.  
  386.     RGBGradCensusCost(cv::Mat left, cv::Mat right, int censusWindow, float lambda) : CostFunction(left, right, lambda) {
  387.         // census.... nimmt einen Block
  388.         this->censusWindow = censusWindow;
  389.         this->censusMargin = censusWindow / 2;
  390.         normWin = censusWindow*censusWindow*3;
  391.         // nimmt einen Block
  392.         l_grad = getRGBGradientAngle(left);
  393.         r_grad = getRGBGradientAngle(right);
  394.     }
  395.  
  396.     bool imageType(cv::Mat left, cv::Mat right) {
  397.         assert(left.type() == right.type() && "imgL imgR types not equal");
  398.         assert(left.type() == CV_8UC3 && "img type not supported");
  399.  
  400.         return true;
  401.     }
  402.  
  403.     unsigned int census(int x1, int x2, int y, cv::Vec3f c1, cv::Vec3f c2) {
  404.         unsigned int diff = 0;
  405.  
  406.         for(int i = y - censusMargin; i <= y + censusMargin; ++i) {
  407.             cv::Vec3f* lptr = l_grad.ptr<cv::Vec3f>(i);
  408.             cv::Vec3f* rptr = r_grad.ptr<cv::Vec3f>(i);
  409.  
  410.             for(int j = -censusMargin; j <= censusMargin; ++j) {
  411.                 cv::Vec3f cl = lptr[x1 + j];
  412.                 cv::Vec3f cr = rptr[x2 + j];
  413.  
  414.                 for(int ch = 0; ch < 3; ++ch) {
  415.                     bool t1 = (c1[ch] < cl[ch]);
  416.                     bool t2 = (c2[ch] < cr[ch]);
  417.  
  418.                     if(t1 != t2) diff++;
  419.                 }
  420.             }
  421.         }
  422.  
  423.         return diff;
  424.     }
  425.  
  426.     float aggregate(int x1, int x2, int y) {
  427.         float sum = 0;
  428.         /*for(int i = y - margin; i <= y + margin; ++i) {
  429.             cv::Vec3f *lptr = l_grad.ptr<cv::Vec3f>(i);
  430.             cv::Vec3f *rptr = r_grad.ptr<cv::Vec3f>(i);
  432.             for(int j = -margin; j <= margin; ++j)
  433.                 sum += census(x1 + j, x2 + j, i, lptr[x1 + j], rptr[x2 + j]);
  434.         }*/
  435.         cv::Vec3f *lptr = l_grad.ptr<cv::Vec3f>(y);
  436.         cv::Vec3f *rptr = r_grad.ptr<cv::Vec3f>(y);
  437.  
  438.         sum = census(x1, x2, y, lptr[x1], rptr[x2]);
  439.         return sum / normWin;
  440.     }
  441. };

step３.视差空间的构建

DSI（Disparity Space Image）视差空间图像为一个三维矩阵，主要由横轴x,纵轴y，视差搜索范围d构成，传统的DP方法一般就是为了在某固定的Y（也就是某极线上）寻找一条从最左段到最右段的最小代价路径，每条路径的代价为

  • 视差空间的构建

输入参数：

  • imageSize:图像的大小
  • blocksize:块的大小
  • y:某条极线

输出参数：

  • map:某极线上左右图像两两领域相互的代价值
  • 左范围：[margin，imageSize.width - margin]
  • 右范围：[margin，imageSize.width - margin]
  1. cv::Mat BlockMatching::disparitySpace(Size imageSize, int blocksize, int y)
  2. {
  3.     int margin = blocksize / 2;
  4.     int start = margin;
  5.     int stopW = imageSize.width - margin;
  6.     int workSpace = stopW - start;
  7.  
  8.     // leave out the borders
  9.     //Mat map = Mat(workSpace, workSpace, CV_32F);        // not preinitialized.. // numeric_limits<float>::max());
  10.     Mat map = Mat(workSpace, workSpace, CV_32F, numeric_limits<float>::max());
  11.  
  12.     //int dmax = 101;
  13.     for(int x1 = start; x1 < stopW; x1++) {
  14.         float* ptr = map.ptr<float>(x1 - margin);       // [x1 - margin, x2 - margin]
  15.  
  16.         //ptr[max(x1 - 1, start) - margin] = numeric_limits<float>::max();              // fast borders
  17.         //ptr[min(x1 + dmax, stopW - 1) - margin] = numeric_limits<float>::max();
  18.         //for(int x2 = x1; x2 < min(x1 + dmax, stopW); x2++) {
  19.         for(int x2 = start; x2 < stopW; x2++) {
  20.  
  21.             // combine costs
  22.             float cost = 0;
  23.             for(size_t i = 0; i < functions.size(); ++i) {
  24.                 float val = functions[i]->aggregate(x1, x2, y);
  25.                 mins[i] = min(mins[i], val);                        // debug
  26.                 maxs[i] = max(maxs[i], val);                        // debug
  27.                 cost += val;
  28.             }
  29.  
  30.             // x1, x2. Das heißt x1 sind die Zeilen. Wir gehen jedes Mal die Zeilen runter.
  31.             // geht nur von 0 - workspace, deshalb margin abziehen
  32.             //map.at<float>(x1 - margin, x2 - margin) = greySad(leftRoi, rightRoi);
  33.             ptr[x2 - margin] = cost;
  34.         }
  35.     }
  36.     return map;
  37. }

step4:动态规划

主要分为四步：

  1. 设置初始位置的值（已知的值，这里设置的最后一行，最后一列的点为初始值）
  2. 计算边界上的代价（最后一行、最后一列）
  3. 从三个方向（向上、向左、斜向上）计算代价
  4. 记录每一步的方向
  5. 第一行的最小值即为视差点

 这里不容易理解，其实它和Leetcode的最短路径的题是相似的

题目：给定一个矩阵，从左上角开始每次只能向右或者向下移动，最后到达右下角的位置，路径上的所有的数字累加起来作为这条路径的路劲和。要求返回所有路径和中的最小路径和

131
151
421

路径1,3,1,1,1,1是所有路径中路径和最小的，所以返回其和7。

  • 第一步：起始点1
  • 第二步：如果一直向下或者向右：
145
251
621
  • 第三步：计算第二行第二列的值，比较第一行第二列和第二行第一列的数谁小就选谁，依次类推。
145
27
6
145
276
68
145
276
687

  •   最小路径和的代码

  1. int minPathSum1(int matrix[][col], int dp[][col], int path[][col])
  2. {
  3.     if(matrix == NULL)
  4.     {
  5.         return 0;
  6.     }
  7.  
  8.     dp[0][0] = matrix[0][0];
  9.  
  10.     //计算第一列的值
  11.     for(int i = 1; i < row; i ++)
  12.     {
  13.         dp[i][0] = dp[i - 1][0] + matrix[i][0];
  14.  
  15.         path[i][0] = 0;
  16.     }
  17.  
  18.     //计算第一行的值
  19.     for(int j = 1; j < col; j++)
  20.     {
  21.         dp[0][j] = dp[0][j- 1] + matrix[0][j];
  22.  
  23.         path[0][j] = 1;
  24.     }
  25.  
  26.     //计算其它的值
  27.     for(int i = 1; i < row; i++)
  28.     {
  29.         for(int j = 1; j < col; j++)
  30.         {
  31.             int direction = dp[i][j-1] < dp[i-1][j] ? 1 : 0;
  32.             dp[i][j] = (direction ?  dp[i][j-1] : dp[i-1][j]) + matrix[i][j];
  33.             path[i][j] = direction;
  34.         }
  35.     }//for
  36.  
  37.  
  38.     return dp[row - 1][col - 1];
  39. }

  •  dp视差空间代码
  1. // (1) set last entry sink in matrix (last value)
  2. // (2-3) Initializ Edges
  3. //    (2) initialize travelpath for last col (only south direction)
  4. //    (3) initialize travelpath for last row (only east direction)
  5. // (4) calculate paths till last sink (last entry) till xLast - 1, yLast - 1
  6. // (-) save all (chosen) directions along the way
  7. void DPmat::preCalc(Mat &matrix, Mat &sum, Mat &dirs) {
  8.     float occlusion_south = 1.0f;
  9.     float occlusion_east = 1.0f;
  10.     sum = Mat::zeros(matrix.rows, matrix.cols, matrix.type());         // not initialized with zero, should not be a problem,
  11.     dirs = Mat::zeros(matrix.rows, matrix.cols, CV_16U);               // because traversion is pre initialized with borders
  12.  
  13.     // dirs = (1, south), (2, south-east), (3, east)
  14.  
  15.     int rowLast = matrix.rows - 1;           // last index inclusive
  16.     int colLast = matrix.cols - 1;           // last index inclusive
  17.  
  18.     // (1) initialize sink (last Entry/ terminal point/ matrix exit value)
  19.     sum.at<float>(rowLast, colLast) = matrix.at<float>(rowLast, colLast);
  20.  
  21.     // (2-3) Initialize Edges
  22.  
  23.     // (2) calculate all last row entries down to exit value | only downward directions (so upward pre calculation)
  24.     for(int y = rowLast - 1; y >= 0; --y) {
  25.         // sum[y,x] = M[y,x] * occlusion_south + sum[y+1,x]
  26.         sum.at<float>(y, colLast) = matrix.at<float>(y, colLast) * occlusion_south + sum.at<float>(y + 1, colLast);   // add current value + successor min(value)
  27.         dirs.at<ushort>(y, colLast) = 1;    // south
  28.     }
  29.  
  30.     // (3) initialize last
  31.     for(int x = colLast - 1; x >= 0; --x) {
  32.         // sum[y,x] = M[y,x] * occlusion_east + sum[y+1,x]
  33.         sum.at<float>(rowLast, x) = matrix.at<float>(rowLast, x) * occlusion_east + sum.at<float>(rowLast, x + 1);   // add current value + successor min(value)
  34.         dirs.at<ushort>(rowLast, x) = 3;    // east
  35.     }
  36.  
  37.     // (4) Main calculation (3 way [south(s), east(e), south-east(se)])
  38.     for(int y = rowLast - 1; y >= 0; --y) {
  39.         float* sum_ptr = sum.ptr<float>(y);
  40.         float* sum_south_ptr = sum.ptr<float>(y+1);
  41.         float* mat_ptr = matrix.ptr<float>(y);
  42.         ushort* dirs_ptr = dirs.ptr<ushort>(y);
  43.  
  44.         for(int x = colLast - 1; x >= 0; --x) {
  45.             // dirs
  46.             //float s = sum.at<float>(y + 1, x);
  47.             //float se = sum.at<float>(y + 1, x + 1);
  48.             //float e = sum.at<float>(y, x + 1);
  49.             float s = sum_south_ptr[x] * occlusion_south;        // (y+1,x)     occlusion dir
  50.             float se = sum_south_ptr[x + 1];                     // (y+1,x+1)
  51.             float e = sum_ptr[x + 1] * occlusion_east;           // (y, x+1)    occlusion dir
  52.  
  53.             // lowest cost till current point
  54.             float p = min(s, min(se, e));
  55.  
  56.             //sum.at<float>(y, x) = p + matrix.at<float>(y, x);   // sum till current (cost + lowest path)
  57.             sum_ptr[x] = p + mat_ptr[x];        // sum[y,x] = p + mat[y, x]
  58.  
  59.             // selection for traversion direction
  60.             //if(p == s) dirs.at<ushort>(y, x) = 1;   // occlusion
  61.             //if(p == se) dirs.at<ushort>(y, x) = 2;   // math
  62.             //if(p == e) dirs.at<ushort>(y, x) = 3;   // occlusion
  63.  
  64.             if(p == s) dirs_ptr[x] = 1;   // occlusion
  65.             if(p == se) dirs_ptr[x] = 2;   // math
  66.             if(p == e) dirs_ptr[x] = 3;   // occlusion
  67.         }
  68.     }
  69. }
  1. void DPmat::disparityFromDirs(Mat &sum, Mat &dirs, Mat &disp, int line, int offset) {
  2.     assert(dirs.type() == CV_16U);
  3.  
  4.     // wir bekommen jetzt einen index x, y
  5.     int rowLast = dirs.rows - 1;
  6.     int colLast = dirs.cols - 1;
  7.  
  8.     int lastval = -1;
  9.     int x1 = 0;
  10.     int x2 = 0;
  11.  
  12.     float minVal = numeric_limits<float>::max();
  13.     int minIndex = 0;
  14.     // seek top entry
  15.     for(x2 = 0; x2 < sum.cols; ++x2) {
  16.         float val = sum.at<float>(x1, x2);
  17.         if(val > minVal) {
  18.             minIndex = x2;
  19.             minVal = val;
  20.         }
  21.     }
  22.  
  23.     x2 = minIndex;
  24.  
  25.     // safe x1, x2 as disparity match
  26.     ushort disparity = abs(x2 - x1);
  27.     ushort* disp_ptr = disp.ptr<ushort>(line);
  28.  
  29.     disp_ptr[x1 + offset] = disparity;
  30.  
  31.     while(x1 < rowLast && x2 < colLast) {
  32.         ushort d = dirs.at<ushort>(x1, x2);
  33.  
  34.         if(d == 1) {    // 1 = down, skipping left index, left got occloded (occlusion from right)
  35.             x1++;
  36.             if(lastval >= 0) disp_ptr[x1 + offset] = lastval;   // dips[line, x1 + offset] = lastval
  37.             //disp_ptr[x1 + offset] = 0;
  38.         }
  39.         if(d == 2) { // match
  40.             // next entry will be match
  41.             x1++;
  42.             x2++;
  43.             disparity = abs(x2 - x1);
  44.  
  45.             disp_ptr[x1 + offset] = disparity;
  46.             lastval = disparity;
  47.         }
  48.         if(d == 3) { // 2 = right, skipping right index, occlusion don't care..
  49.             x2++;
  50.             if(lastval >= 0) disp_ptr[x1 + offset] = lastval;   // dips[line, x1 + offset] = lastval
  51.             //disp_ptr[x1 + offset]= 0;
  52.         }
  53.     }
  54. }

原图像：

视差图像：


结论：

这种只考虑了图像左右相邻像素的视差约束，而忽略了上下领域像素之间的视差约束，这种方法得到的解因此也称为扫描线最优解，不能称为严格意义上的最优解，视差图像也出现了明显的横向条纹瑕疵。而且其计算时间也比较慢

声明：本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：【wpsshop博客】
推荐阅读
相关标签

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

  

闽ICP备14008679号