ITdaan
搜索本文相关内容

OpenCV學習筆記（19）雙目測距與三維重建的OpenCV實現問題集錦（四）三維重建與OpenGL顯示

本文转载自  rocky69  查看原文  2012/07/08  1    matlab/ float/ opengl/ opencv/ open           

http://blog.csdn.net/chenyusiyuan/article/details/5970799

五、三維重建與 OpenGL 顯示

.

在獲取到視差數據後，利用 OpenCV 的 reProjectImageTo3D 函數結合 Bouquet 校正方法得到的 Q 矩陣就可以得到環境的三維坐標數據，然後利用 OpenGL 來實現三維重構。 OpenCV 與 OpenGL 的編程范例，我在 學習筆記（ 15 ） 中有詳細的討論，這裡就不重復了，下面補充一些細節問題：

.

. 

1．             reProjectImageTo3D 是怎樣計算出三維坐標數據的？

 

圖 22

.

 

相信看過 OpenCV 第 12 章的朋友對上圖中的 Q 矩陣不會陌生，根據以上變換公式，按理說 OpenCV 應該也是通過矩陣運算的方式來計算出三維坐標數據的，但實際上仔細查看源代碼，會發現 cvReprojectImageTo3D 用了比較奇怪的方法來實現，主要代碼如下：

 

[cpp] view plaincopy
  1. 02737     for( y = 0; y < rows; y++ )  
  2. 02738     {  
  3. 02739         const float* sptr = (const float*)(src->data.ptr + src->step*y);   // 視差矩陣指針  
  4. 02740         float* dptr0 = (float*)(dst->data.ptr + dst->step*y), *dptr = dptr0;   // 三維坐標矩陣指針  
  5. // 每一行運算開始時，用 當前行號y 乘以Q陣第2列、再加上Q陣第4列，作為初始值  
  6. // 記 qq=[qx, qy, qz, qw]』  
  7. 02741         double qx = q[0][1]*y + q[0][3], qy = q[1][1]*y + q[1][3];      
  8. 02742         double qz = q[2][1]*y + q[2][3], qw = q[3][1]*y + q[3][3];     
  9. …   
  10. // 每算完一個像素的三維坐標，向量qq 累加一次q陣第1列  
  11. // 即：qq = qq + q(:,1)  
  12. 02769         for( x = 0; x < cols; x++, qx += q[0][0], qy += q[1][0], qz += q[2][0], qw += q[3][0] )  
  13. 02770         {  
  14. 02771             double d = sptr[x];  
  15. // 計算當前像素三維坐標  
  16. // 將向量qq 加上 Q陣第3列與當前像素視差d的乘積，用所得結果的第4元素除前三位元素即可  
  17. // [X,Y,Z,W]』 = qq + q(:,3) * d;   iW = 1/W; X=X*iW; Y=Y*iW; Z=Z*iW;  
  18. 02772             double iW = 1./(qw + q[3][2]*d);  
  19. 02773             double X = (qx + q[0][2]*d)*iW;  
  20. 02774             double Y = (qy + q[1][2]*d)*iW;  
  21. 02775             double Z = (qz + q[2][2]*d)*iW;  
  22. 02776             if( fabs(d-minDisparity) <= FLT_EPSILON )  
  23. 02777                 Z = bigZ;   // 02713     const double bigZ = 10000.;  
  24. 02778   
  25. 02779             dptr[x*3] = (float)X;  
  26. 02780             dptr[x*3+1] = (float)Y;  
  27. 02781             dptr[x*3+2] = (float)Z;  
  28. 02782         }  

 

OpenCV 的這種計算方式比較令人費解，我的理解是可能這種方式的計算速度比較快。理論上，直接通過矩陣 Q 與向量[x,y,d,1]』 的乘積就可以得到相同的結果，下面用 Matlab 來驗證一下兩種方式是異曲同工的，用 Matlab 按照 OpenCV 計算方式得到的結果稱為「 OpenCV method 」，直接按公式計算得到的結果稱為「 Equation method 」，用 OpenCV 本身算出的三維坐標作為參考，程序代碼如下 ：  

 

 

[c-sharp] view plaincopy
  1. close all;clear all;clc   
  2. im = imread('C:/Stereo IO Data/lfFrame_01.jpg');  
  3. data = importdata('C:/Stereo IO Data/disparity_01.txt');  
  4. r = data(1);    % 行數  
  5. c = data(2);    % 列數  
  6. disp = data(3:end); % 視差  
  7. vmin = min(disp);  
  8. vmax = max(disp);  
  9. disp = reshape(disp, [c,r])'; % 將列向量形式的 disp 重構為 矩陣形式  
  10. %  OpenCV 是行掃描存儲圖像，Matlab 是列掃描存儲圖像  
  11. %  故對 disp 的重新排列是首先變成 c 行 r 列的矩陣，然後再轉置回 r 行 c 列  
  12. img = uint8( 255 * ( disp - vmin ) / ( vmax - vmin ) );  
  13. q = [1. 0. 0. -1.5690376663208008e+002;...  
  14.     0. 1. 0. -1.4282237243652344e+002;...      
  15.     0. 0. 0. 5.2004731331639300e+002;...  
  16.     0. 0. 1.0945105843175637e-002 0.]; % q(4,3) 原為負值，現修正為正值  
  17. big_z = 1e5;  
  18. pos1 = zeros(r,c,3);  
  19. pos2 = zeros(r,c,3);  
  20. for i = 1:r  
  21.     qq = q*[0 i 0 1]';  
  22.     for j = 1:c  
  23.         if disp(i,j)>0  
  24.         % OpenCV method  
  25.             vec = qq + q(:,3)*disp(i,j);  
  26.             vec = vec/vec(4);  
  27.             pos1(i,j,:) = vec(1:3);  
  28.         % Textbook method  
  29.             tmp = q*[j,i,disp(i,j),1]'; % j 是列數，i 是行數，分別對應公式中的 x 和 y  
  30.             pos2(i,j,:) = tmp(1:3)/tmp(4);  
  31.         else  
  32.             pos1(i,j,3) = big_z;  
  33.             pos2(i,j,3) = big_z;  
  34.         end  
  35.         qq = qq + q(:,1);  
  36.     end  
  37. end  
  38. subplot(221);  
  39. imshow(im); title('Left Frame');  
  40. subplot(222);  
  41. imshow(img); title('Disparity map');  
  42. % Matlab按OpenCV計算方式得到的三維坐標  
  43. x = pos1(:,:,1);   
  44. y = -pos1(:,:,2);  % 圖像坐標系Y軸是向下為正方向，因此需添加負號來修正  
  45. z = pos1(:,:,3);   
  46. ind = find(z>10000);  % 以毫米為量綱  
  47. x(ind)=NaN; y(ind)=NaN; z(ind)=NaN;  
  48. subplot(234);  
  49. mesh(x,z,y,double(im),'FaceColor','texturemap');  % Matlab 的 mesh、surf 函數支持紋理映射  
  50. colormap(gray);   
  51. axis equal;   
  52. axis([-1000 1000 0 9000 -500 2000]);  
  53. xlabel('Horizonal');ylabel('Depth');zlabel('Vertical'); title('OpenCV method');  
  54. view([0 0]);  % 正視圖  
  55. % view([0 90]);   % 俯視圖  
  56. % view([90 0]);   % 側視圖  
  57. % Matlab 按公式直接計算得到的三維坐標  
  58. x = pos2(:,:,1);   
  59. y = -pos2(:,:,2);   
  60. z = pos2(:,:,3);   
  61. ind = find(z>10000);  % 以毫米為量綱  
  62. x(ind)=NaN; y(ind)=NaN; z(ind)=NaN;  
  63. subplot(235);  
  64. mesh(x,z,y,double(im),'FaceColor','texturemap');   
  65. colormap(gray);   
  66. axis equal;   
  67. axis([-1000 1000 0 9000 -500 2000]);  
  68. xlabel('Horizonal');ylabel('Depth');zlabel('Vertical'); title('Equation method');  
  69. view([0 0]);  
  70. % 讀入OpenCV計算保存到本地的三維坐標作為參考  
  71. data=importdata('C:/Stereo IO Data/xyz.txt');  
  72. x=data(:,1); y=data(:,2); z=data(:,3);  
  73. ind=find(z>1000);  % 以釐米為量綱  
  74. x(ind)=NaN; y(ind)=NaN; z(ind)=NaN;  
  75. x=reshape(x,[352 288])'; % 數據寫入時是逐行進行的，而Matlab是逐列讀取  
  76. y=-reshape(y,[352 288])';   
  77. z=reshape(z,[352 288])';  
  78. subplot(236)  
  79. mesh(x,z, y,double(im),'FaceColor','texturemap');  
  80. colormap(gray);   
  81. axis equal;axis([-100 100 0 900 -50 200]);  
  82. xlabel('Horizonal');ylabel('Depth');zlabel('Vertical'); title('OpenCV result');  
  83. view([0 0]);  

 

 

圖 23

.

.  

2．             為什麼利用修正了的 Q 矩陣所計算得到的三維數據中， Y 坐標數據是正負顛倒的？

 

 

圖 24

.

 

這個問題我覺得可以從圖像坐標系與攝像機坐標系的關系這一角度來解釋。如上圖所示，一般圖像坐標系和攝像機坐標系都是以從左至右為 X 軸正方向，從上至下為 Y 軸正方向 ，攝像機坐標系的 Z 軸正方向則是從光心到成像平面的垂線方向。因此，我們得到的三維坐標數據中 Y 軸數據的正負與實際是相反的，在應用時要添加負號來修正。

.

. 

3．             如何畫出三維重建圖像和景深圖像？

 

.

利用 cvReprojectImageTo3D 計算出的三維坐標數據矩陣一般是三通道浮點型的，需要注意的是這個矩陣存儲的是三維坐標數據，而不是 RGB 顏色值，所以是不能調用 cvShowImage() 或者 OpenCV2.1 版的 imshow() 等函數來顯示這個矩陣，否則就會看到這種圖像：

.

 

圖 25

.

 

這裡出現的明顯的四個色塊，其實應該是由三維坐標數據中的 X 軸和 Y 軸數據造成，不同象限的數據形成相應的色塊。

要畫出正確的三維重建圖像，可以結合 OpenGL （可參考我的 學習筆記（ 15 ） ）或者 Matlab （例如保存三維數據到本地然後用 Matlab 的 mesh 函數畫出，例程見本文問題 1 ；也可以考慮在 OpenCV 中調用 Matlab 混合編程）來實現。

深度圖像的顯示相對比較簡單，只要從三維坐標數據中分離出來（可用 cvSplit() 函數），經過適當的格式轉換（例如轉換為 CV_8U 格式），就可用 cvShowImage() 或者 OpenCV2.1 版的 imshow() 等函數來顯示了，偽彩色的深度圖 也可以參考我的 學習筆記（ 18 ） 問題 6 給出的例程 稍作修改即可實現。

.

. 

4．             怎樣把 OpenGL 窗口的圖像復制到 OpenCV 中用 IplImage 格式顯示和保存？

 

.

在 學習筆記（ 15 ） 中詳細給出了將 OpenCV 生成的 IplImage 圖像和三維坐標數據復制到 OpenGL 中顯示的例程，而在應用中，我們有時候也需要把 OpenGL 實時顯示的三維圖像復制到 OpenCV 中，用 IplImage 格式保存，以便和其它圖像組合起來顯示或保存為視頻文件。這裡給出相應的例程以供參考：

首先在創建 OpenGL 窗口時，顯示模式要如下設置：

 

[c-sharp] view plaincopy
  1. //***OpenGL Window  
  2. glutInit(&argc, argv);  
  3. glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGB);  
  4. glutInitWindowPosition(10,420);  
  5. glutInitWindowSize(glWinWidth, glWinHeight);  
  6. glutCreateWindow("3D disparity image");  

 

在循環中的調用：

 

[c-sharp] view plaincopy
  1.     //////////////////////////////////////////////////////////////////////////  
  2. // OpenGL顯示  
  3.         img3dIpl = img3d;  
  4. load3dDataToGL(&img3dIpl);      // 載入需要顯示的圖像（視差數據）  
  5. loadTextureToGL(&img1roi);      // 顯示紋理  
  6. glutReshapeFunc (reshape);          // 窗口變化時重繪圖像  
  7. glutDisplayFunc(renderScene);       // 顯示三維圖像  
  8. glutPostRedisplay();                // 刷新畫面（不用此語句則不能動態更新圖像）  
  9. loadPixel2IplImage(imgGL);          // 將 OpenGL 生成的像素值存儲到 IplImage 中  
  

 

loadGLPixelToIplImage 函數定義：

 

 

[c-sharp] view plaincopy
  1. //////////////////////////////////////////////////////////////////////////  
  2. // 將OpenGL窗口像素存儲到 IplImage 中  
  3. void  loadGLPixelToIplImage(IplImage* img)  
  4. {  
  5.     const int n = 3*glWinWidth*glWinHeight;   
  6.     float *pixels = (float *)malloc(n * sizeof(GL_FLOAT));   
  7.     IplImage *tmp = cvCreateImage(cvSize(glWinWidth, glWinHeight), 8, 3);  
  8.     tmp->origin = CV_ORIGIN_BL;  
  9. /* 後台緩存的圖像數據才是我們需要復制的，若復制前台緩存會把可能的疊加在OpenGL窗口上的對象（其它窗口或者鼠標指針）也復制進去*/  
  10.     glReadBuffer(GL_BACK);        
  11.     glReadPixels(0, 0, glWinWidth, glWinHeight, GL_RGB, GL_FLOAT, pixels);  
  12.     int k = 0;  
  13.     for(int i = 0 ; i < glWinHeight; i++)  
  14.     {  
  15.         for(int j = 0 ; j < glWinWidth; j++,k+=3)  
  16.         {  
  17.             CvPoint pt = {j, glWinHeight - i - 1};  
  18.             uchar* temp_ptr = &((uchar*)(tmp->imageData + tmp->widthStep*pt.y))[pt.x*3];  
  19.             //OpenGL采用的是BGR格式，所以,讀出來以後,還要換一下R<->B,才能得到RGB  
  20.             temp_ptr[0] = pixels[k+2] * 255; //blue   
  21.             temp_ptr[1] = pixels[k+1] * 255; //green  
  22.             temp_ptr[2] = pixels[k] * 255;   //red  
  23.         }  
  24.     }     
  25.     cvResize(tmp, img);  
  26.     // 釋放內存  
  27.     free(pixels);  
  28.     cvReleaseImage(&tmp);  
  29. }  

 

 

顯示效果如下：

 

圖26 

 

 

 

（待續）


智能推荐

注意！

本站转载的文章为个人学习借鉴使用，本站对版权不负任何法律责任。如果侵犯了您的隐私权益，请联系我们删除。



猜您在找
OpenCV學習筆記（17）雙目測距與三維重建的OpenCV實現問題集錦（二）雙目定標與雙目校正 OpenCV學習筆記（16）雙目測距與三維重建的OpenCV實現問題集錦（一）圖像獲取與單目定標 OpenCV學習筆記（18）雙目測距與三維重建的OpenCV實現問題集錦（三）立體匹配與視差計算 OpenCV學習筆記（7）圖像的通道（channels）問題 OpenCV學習筆記（13）基於OpenCV2.1的MFC雙目視覺平台構建
 
© 2014-2019 ITdaan.com 粤ICP备14056181号  

赞助商广告