[关闭]
@TangWill 2019-07-26T12:10:29.000000Z 字数 8433 阅读 823

重建实验代码及结果

大创

一、条形图案生成代码

  1. # -*- coding: utf-8 -*-
  3. from PIL import Image
  4. import os
  7. def de_bruijn(k, n):
  8.     """
  9.     de Bruijn sequence for alphabet k
  10.     and subsequences of length n.
  11.     """
  12.     try:
  13.         # let's see if k can be cast to an integer;
  14.         # if so, make our alphabet a list
  15.         _ = int(k)
  16.         alphabet = list(map(str, range(k)))
  18.     except (ValueError, TypeError):
  19.         alphabet = k
  20.         k = len(k)
  22.     a = [0] * k * n
  23.     sequence = []
  25.     def db(t, p):
  26.         if t > n:
  27.             if n % p == 0:
  28.                 sequence.extend(a[1:p + 1])
  29.         else:
  30.             a[t] = a[t - p]
  31.             db(t + 1, p)
  32.             for j in range(a[t - p] + 1, k):
  33.                 a[t] = j
  34.                 db(t + 1, t)
  36.     db(1, 1)
  37.     return "".join(alphabet[i] for i in sequence)
  39. if __name__ == "__main__":
  40.     DB = de_bruijn(3, 3)
  41.     DB_list = list(DB)
  42.     print(DB_list)
  44.     x = 912
  45.     y = 1140
  47.     c = Image.new("RGB", (x, y))
  49.     skip = False
  50.     i = 0
  51.     while i < x:
  52.         for k in range(0, 15):
  53.             rgb = 255 - (5*(k - 7) ** 2)
  54.             temp = int(i / 20)
  55.             if temp < 27:
  56.                 for j in range(0, y):
  57.                     if i % 20 == 0:
  58.                         for m in range(0, 5):
  59.                             c.putpixel([i + m + 186, j], (0, 0, 0))
  60.                         i = i + 5
  61.                         skip = False
  62.                     elif DB_list[temp] == "0":
  63.                         c.putpixel([i + k + 186, j], (rgb, 0, 0))
  64.                         skip = True
  65.                     elif DB_list[temp] == "1":
  66.                         c.putpixel([i + k + 186, j], (0, rgb, 0))
  67.                         skip = True
  68.                     elif DB_list[temp] == "2":
  69.                         c.putpixel([i + k + 186, j], (0, 0, rgb))
  70.                         skip = True
  71.         if skip:
  72.             i = i + 15
  73.         else:
  74.             i = i + 1
  76.     c.show()
  77.     c.save("structured_light .bmp")

二、重建代码

main.m

  1. clear all; clc;  %清空变量
  2. %--------------------------------------------------------------
  3. % 0:red,1:green,2:blue
  4. % num: 为debruijn序列，包含红绿蓝3种颜色，window size为3
  5. %      num序列总长度为3^3 = 27
  6. % 原投影图案生成方法: 
  7. % (1).原投影图片大小竖直方向1140像素，水平方向912像素
  8. % (2).每个条纹水平方向长度15个像素，第一个条纹中心位于第8个像素位置
  9. %     因此条纹中心像素坐标为8+(i-1)*20
  10. % map: map(i,j,k)表示连续三个颜色为(i,j,k)的三元组i颜色对应的条纹中心位置
  11. %      例如map(0,0,0)=8，因为三个连续红条纹在开始的位置，第一个红条纹中心为8
  12. %          map(1,0,0)=68,  绿红红位于第4个位置，68为绿条纹中心 
  13. %--------------------------------------------------------------
  14. num = [0, 0, 0, 1, 0, 0, 2, 0, 1, 1, 0, 1, 2, 0, 2, 1, 0, 2, 2, 1, 1, 1, 2, 1, 2, 2, 2];
  15. map = zeros(3,3,3);
  16. for i=1:25
  17.     map(num(i)+1,num(i+1)+1,num(i+2)+1) = 8 + (i-1)*20;
  18. end
  20. %--------------------------------------------------------------
  21. % R12: 由相机到投影仪的旋转矩阵 3*3
  22. % T12: 由相机到投影仪的平移矩阵 1*3
  23. % Kc1: 相机内参矩阵 3*3
  24. % Kp2: 投影仪内参矩阵 3*3
  25. % Hc1: 相机单应性矩阵 3*4
  26. % Hp2: 投影仪单应性矩阵 3*4
  27. %--------------------------------------------------------------
  28. R12 = [9.7516454746001247e-001 -7.7281024823287244e-003 2.2134674564838733e-001;
  29.        5.4950021604755926e-004 9.9947240635257240e-001 3.2474713075112072e-002;
  30.     -2.2148093242223318e-001 -3.1546558795230203e-002 9.7465430343355963e-001];
  31. T12 = [-1.4326112323874105e+002 2.0614986909838471e+001 -1.0037865572591169e+001];
  32. Kc1 = [2.1363192826189816e+003 0. 6.4087239901301598e+002;
  33.        0. 2.1377724542629630e+003 4.9766519824009470e+002;
  34.        0. 0. 1.];
  35. Kp2 = [1.7342859113112791e+003 0. 4.5022560667574703e+002;
  36.        0. 3.4614881021035653e+003 5.2438368036642009e+002;
  37.        0. 0. 1.];
  38. Hp2 = Kp2*[R12,T12'];
  39. Hc1 = Kc1*[eye(3,3),zeros(3,1)];
  41. img12 = imread('2.png'); %读取图片
  42. imshow(img12);
  43. [flag12,color_flag12] = feature_point(img12); %提取特征，详情见函数内部注释
  44. figure;
  45. imshow(color_flag12);
  46. [num_coord12, coord12] = reconstruction(color_flag12,Hc1,Hp2,map); %三维坐标重建，详情见函数内部注释
  47. figure;
  48. scatter3(coord12(:,1),coord12(:,2),coord12(:,3),2) %画三维离散点图

feature.m

  1. %--------------------------------------------------------------
  2. %条纹特征点提取，包括条纹中心点的位置和对应的颜色
  3. %1:red,2:green,3:blue
  4. %输入参数IMG: 大小n*m*3，表示输入三通道彩色图片
  5. %输出参数flag: 大小n*m, 表示对应像素是否是条纹中心
  6. %              flag(i,j)=1表示(i,j)位置是条纹中心，flag的值为0或1
  7. %输出参数color_flag: 大小n*m，表示对应条纹中心的颜色
  8. %                    color_flag(i,j)=2代表(i,j)位置条纹中心颜色为绿色
  9. %                    color_flag(i,j)=0代表该位置不是条纹中心
  10. %--------------------------------------------------------------
  11. function [flag, color_flag] = feature_point(IMG)
  12.     IMG_gray = rgb2gray(IMG); %彩色图转灰度图
  13.     flag = zeros(1024,1280); %初始化
  14.     color_flag = zeros(1024,1280); %初始化
  15.     linyu = 5; %一个阈值，可自取
  16.     %i和j的范围代表物体在图片的这个范围内，限制范围可以增加计算速度，可自取
  17.     for i=300:880
  18.         for j=500:1000
  19.             if IMG_gray(i,j)<10 %如果灰度图的对应像素值小于10认为该点为黑色背景
  20.                 continue;
  21.             end
  22.             % temp为以(i,j)为中心长度为2*linyu+1的图像灰度数组
  23.             temp = IMG_gray(i,j-linyu:j+linyu);
  24.             % 如果当前中心点的灰度值是这个邻域中的最大值，那么认为该点是条纹中心
  25.             [~,da] = max(temp);
  26.             if da == linyu+1
  27.                 flag(i,j)=1;
  28.             end
  29.         end
  30.     end
  31.     for i=300:880
  32.         for j=500:1000
  33.             if flag(i,j)==0 %如果不是条纹中心则跳过
  34.                 continue;
  35.             end
  36.             %灰度图片，如果灰度值大于150认为该点是白色条纹中心
  37. %             if IMG_gray(i,j)>150 
  38. %                color_flag(i,j)=4;
  39. %                continue;
  40. %             end
  41.             %彩色图片，如果第一个通道比另外两个通道值大，认为该位置是红色条纹中心
  42.             if IMG(i,j,1)>IMG(i,j,2) && IMG(i,j,1)>IMG(i,j,3)
  43.                 color_flag(i,j)=1;
  44.                 continue;
  45.             end
  46.             %彩色图片，如果第二个通道比另外两个通道值大，认为该位置是绿色条纹中心
  47.             if IMG(i,j,2)>IMG(i,j,1) && IMG(i,j,2)>IMG(i,j,3)
  48.                 color_flag(i,j)=2;
  49.                 continue;
  50.             end
  51.             %彩色图片，如果第三个通道比另外两个通道值大，认为该位置是蓝色条纹中心
  52.             if IMG(i,j,3)>IMG(i,j,2) && IMG(i,j,3)>IMG(i,j,1)
  53.                 color_flag(i,j)=3;
  54.                 continue;
  55.             end
  56.         end
  57.     end
  58. end

reconstruction.m

  1. %--------------------------------------------------------------
  2. %物体三维坐标点重建
  3. %输入参数color_flag: 对应feature_point函数的第二个输出
  4. %输入参数Hc: 相机单应性矩阵
  5. %输入参数Hp: 投影仪单应性矩阵
  6. %输入参数map: 主函数开始的那个map数组
  7. %输出参数num_coord: 输出三维点的数目
  8. %输出参数coord: 输出三维点的坐标
  9. %--------------------------------------------------------------
  10. function [num_coord, coord] = reconstruction(color_flag,Hc,Hp,map)
  11.     num_coord = 0;    %初始化
  12.     for i=300:880
  13.         ind = 0;
  14.         loc = zeros(1,100);
  15.         %将图片第i行的所有条纹中心临时存储在loc中
  16.         for j=500:1000
  17.             if color_flag(i,j)~=0
  18.                 ind = ind+1;
  19.                 loc(ind) = j;
  20.             end
  21.         end
  22.         if (ind<3) % 如果少于3个中心则无法重建
  23.             continue;
  24.         end
  25.         for j = 1:ind-2
  26.             A = zeros(3,3);
  27.             num_coord = num_coord+1;
  28.             A(1,:) = Hc(3,1:3).*loc(j)-Hc(1,1:3);
  29.             A(2,:) = Hc(3,1:3).*i-Hc(2,1:3);
  30.             A(3,:) = Hp(3,1:3).*map(color_flag(i,loc(j)),color_flag(i,loc(j+1)),color_flag(i,loc(j+2)))-Hp(1,1:3);
  31.             ydx = A\[Hc(1,4)-Hc(3,4)*loc(j);Hc(2,4)-Hc(3,4)*i;Hp(1,4)-Hp(3,4)*map(color_flag(i,loc(j)),color_flag(i,loc(j+1)),color_flag(i,loc(j+2)))];
  32.             % 设定一下和相机的距离，如果过近或者过远显然都是噪声点，可自取
  33. %             if ydx(3)>1000 || ydx(3)<500
  34. %                 continue;
  35. %             end
  36.             coord(num_coord,:) = ydx';
  37.         end
  38.         % debruijn序列最后两个条纹仍利用最后三个颜色重建
  39.         for j=ind-1:ind
  40.             if j==ind-1
  41.                 A = zeros(3,3);
  42.                 num_coord = num_coord+1;
  43.                 A(1,:) = Hc(3,1:3).*loc(j)-Hc(1,1:3);
  44.                 A(2,:) = Hc(3,1:3).*i-Hc(2,1:3);
  45.                 A(3,:) = Hp(3,1:3).*(map(color_flag(i,loc(j-1)),color_flag(i,loc(j)),color_flag(i,loc(j+1)))+15)-Hp(1,1:3);
  46.                 ydx = A\[Hc(1,4)-Hc(3,4)*loc(j);
  47.                          Hc(2,4)-Hc(3,4)*i;
  48.                          Hp(1,4)-Hp(3,4)*(map(color_flag(i,loc(j-1)),color_flag(i,loc(j)),color_flag(i,loc(j+1)))+15)];
  49. %                 if ydx(3)>1000 || ydx(3)<500
  50. %                     continue;
  51. %                 end
  52.                 coord(num_coord,:) = ydx';
  53.             elseif j==ind
  54.                 A = zeros(3,3);
  55.                 num_coord = num_coord+1;
  56.                 A(1,:) = Hc(3,1:3).*loc(j)-Hc(1,1:3);
  57.                 A(2,:) = Hc(3,1:3).*i-Hc(2,1:3);
  58.                 A(3,:) = Hp(3,1:3).*(map(color_flag(i,loc(j-2)),color_flag(i,loc(j-1)),color_flag(i,loc(j)))+30)-Hp(1,1:3);
  59.                 ydx = A\[Hc(1,4)-Hc(3,4)*loc(j);
  60.                          Hc(2,4)-Hc(3,4)*i;
  61.                          Hp(1,4)-Hp(3,4)*(map(color_flag(i,loc(j-2)),color_flag(i,loc(j-1)),color_flag(i,loc(j)))+30)];
  62. %                 if ydx(3)>1000 || ydx(3)<500
  63. %                     continue;
  64. %                 end
  65.                 coord(num_coord,:) = ydx';
  66.             end
  67.         end
  68.     end
  69. end

未注释部分参考公式：

式中，表示的是矩阵的第行和第列元素。将两个式子整理后得到关于、和的三元一次方程，

三、重建结果

3.1 平板

3.2 球面

附录

附录一 相机及投影仪标定参数

  1. <?xml version="1.0"?>
  2. <opencv_storage>
  3. <Kc type_id="opencv-matrix">
  4.   <rows>3</rows>
  5.   <cols>3</cols>
  6.   <dt>d</dt>
  7.   <data>
  8.     2.1363192826189816e+003 0. 6.4087239901301598e+002 0.
  9.     2.1377724542629630e+003 4.9766519824009470e+002 0. 0. 1.</data></Kc>
  10. <kc type_id="opencv-matrix">
  11.   <rows>5</rows>
  12.   <cols>1</cols>
  13.   <dt>d</dt>
  14.   <data>
  15.     -5.6850946659891023e-002 2.6322763056141931e+000
  16.     -1.1377015012654045e-003 1.6506197473016759e-003
  17.     -2.1654415041354078e+001</data></kc>
  18. <Kp type_id="opencv-matrix">
  19.   <rows>3</rows>
  20.   <cols>3</cols>
  21.   <dt>d</dt>
  22.   <data>
  23.     1.7342859113112791e+003 0. 4.5022560667574703e+002 0.
  24.     3.4614881021035653e+003 5.2438368036642009e+002 0. 0. 1.</data></Kp>
  25. <kp type_id="opencv-matrix">
  26.   <rows>5</rows>
  27.   <cols>1</cols>
  28.   <dt>d</dt>
  29.   <data>
  30.     9.6270780483633647e-002 -8.0172849754942543e+000
  31.     -2.6246483366578217e-003 2.3503112604238339e-003
  32.     2.2313956040160841e+002</data></kp>
  33. <Rp type_id="opencv-matrix">
  34.   <rows>3</rows>
  35.   <cols>3</cols>
  36.   <dt>d</dt>
  37.   <data>
  38.     9.8817697193050202e-001 -2.7706101018562532e-003
  39.     1.5329251731877641e-001 1.1226239513894419e-002
  40.     9.9846035682073886e-001 -5.4322071055690473e-002
  41.     -1.5290599626123413e-001 5.5400718199914163e-002
  42.     9.8668663552836944e-001</data></Rp>
  43. <Tp type_id="opencv-matrix">
  44.   <rows>3</rows>
  45.   <cols>1</cols>
  46.   <dt>d</dt>
  47.   <data>
  48.     -1.5168118625605635e+002 1.7918489397033632e+001
  49.     -5.3101094723810185e+001</data></Tp>
  50. <cam_error>6.0649322922523355e-002</cam_error>
  51. <proj_error>1.5337621234636628e-001</proj_error>
  52. <stereo_error>1.3013498108534971e-001</stereo_error>
  53. <frameWidth>1280</frameWidth>
  54. <frameHeight>1024</frameHeight>
  55. <screenResX>0</screenResX>
  56. <screenResY>0</screenResY>
  57. </opencv_storage>
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注