查看: 1093|回复: 19|关注: 0

[已解决] MATLAB提取点云圆柱的中心轴线

[复制链接]

新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
请教论坛里的大神们,去和用MATLAB提取点云圆柱的中心轴线呢?附件里给出的是一堆点云数据,为一个不完整的圆柱模型,我想得到其中心轴线的输出,不知道怎么来实现,希望有人可以尝试做下,讲下原理并分享下程序,小白在此先谢过了,感激!!!

圆柱.txt

140.69 KB, 下载次数: 27

新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
 楼主| 发表于 2018-5-16 09:50:10 | 显示全部楼层
可不可以不要沉:'(

论坛优秀回答者

入门

435 麦片

财富积分


50500


2

主题

432

帖子

95

最佳答案
  • 关注者: 24
发表于 2018-5-17 09:47:38 | 显示全部楼层
利用"圆柱柱面上的点到轴线的距离相等"这一特性来解决.
最简单,但不是通用型的方法,就是投影作图法.除了垂直于轴线的平面,圆柱柱面主体(即不考虑投影为椭圆的端口部分)在其他平面上的投影都是矩形,轴线必定位于该矩形的对称轴上.这个方法要求数据点足够多,能够覆盖住投影矩形,如果投影矩形不完整的话,这个方法就不能用了.
CylinderFitting01.png
比如在x-y投影面上,通过对较小区间Δx内取最小值yl和最大值yu即可确定投影矩形的边界,对一组x及 (yu-yl)/2用polyfit拟合即可得到轴线在x-y平面上的投影方程.在y-z平面上做相同处理,得到轴线的另一个投影方程,两方程联立,即所求轴线方程.
这个方法优点是直观易行,缺点除了上面说的对数据分布有要求之外,拟合中也可能引起误差,导致结果的精度略差.

另一个方法的原理是,只有在直线为轴线时,各散点到该直线的距离的方差才能达到最小值.用粒子群算法来优化直线参数,使距离的方差尽可能的小,就可以不断逼近轴线.

                               
登录/注册后可看大图

优点是精度较高,对散点投影是否具有规则形状也没有要求,缺点是如果初始参数设置得不好,可能需要算很久.
于是可以先用法一求出大致的轴线方程,将其参数带入法二作为初值来提高精度.图中两种方法的轴线大致重合略有偏差,法二距离的方差比法一要小很多.


CylinderFitting.png

论坛优秀回答者

15

主题

1851

帖子

115

最佳答案
  • 关注者: 149
发表于 2018-5-17 10:18:20 | 显示全部楼层
为利于显示图形,y=y+1415;
SSy =   1.3251e+06
fx1=@(b,x1,x2)b(1)+b(2)*x1+b(3)*x2+sqrt(b(4)+b(5)*x1+b(6)*x2+b(7)*x1.^2+b(8)*x2.^2+b(9)*x1.*x2);
fx2@(b,x1,x2)=b(1)+b(2)*x1+b(3)*x2-sqrt(b(4)+b(5)*x1+b(6)*x2+b(7)*x1.^2+b(8)*x2.^2+b(9)*x1.*x2);

b =[-21.70425549  -0.5997414942  0.2697712017  -206208.4142  425.4930561  975.4828918  -0.2217687332  -1.148980795  -1.001140106]
RSS = 469.400533058
MSe =  0.12487
R^2 =  0.99965
nh1480.png

新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
 楼主| 发表于 2018-5-18 16:49:10 | 显示全部楼层
TouAkira 发表于 2018-5-17 09:47
利用"圆柱柱面上的点到轴线的距离相等"这一特性来解决.
最简单,但不是通用型的方法,就是投影作图法.除了垂 ...

您好,请问方便分享下完整的代码吗,感谢!

新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
 楼主| 发表于 2018-5-18 16:50:19 | 显示全部楼层
stats01 发表于 2018-5-17 10:18
为利于显示图形,y=y+1415;
SSy =   1.3251e+06
fx1=@(b,x1,x2)b(1)+b(2)*x1+b(3)*x2+sqrt(b(4)+b(5)*x1+b( ...

老师您好,请问方便分享下完整的代码吗?这样看看不太懂呢

论坛优秀回答者

入门

435 麦片

财富积分


50500


2

主题

432

帖子

95

最佳答案
  • 关注者: 24
发表于 2018-5-19 09:22:11 | 显示全部楼层
云妹y 发表于 2018-5-18 04:49
您好,请问方便分享下完整的代码吗,感谢!
  1. clear all;
  2. data = importdata('圆柱.txt');
  3. xx = data(:,1); yy = data(:,2); zz = data(:,3);
  4. [sx,idx] = sort(xx); [sz,idz] = sort(zz);
  5. %% view(2) 对x 在[40,90],z为任意值,搜索y的最大最小值,构成直线x-y关系
  6. Xst = 40; Xed = 90; Xstep = 1;% 起点,终点,步长
  7. n = 50;xylx = zeros(n,1);xyly = zeros(n,1);
  8. lsytx = yy(idx);%循环找出y的最小值构成的边界
  9. for k = 1:1:n
  10.     xylx(k) = Xst + (k-0.5)*Xstep;
  11.     ct = [];
  12.     for i =1:1:length(sx)
  13.         if sx(i) > Xst + (k-2)*Xstep && sx(i) < Xst + k*Xstep
  14.             ct = [ct;i];
  15.         else
  16.             continue;
  17.         end
  18.     end
  19.     xyly(k) = min(lsytx(ct));
  20. end

  21. xyux = zeros(n,1);xyuy = zeros(n,1);
  22. usytx = yy(idx);%循环找出y的最大值构成的边界
  23. for k = 1:1:n
  24.     xyux(k) = Xst + (k-0.5)*Xstep;
  25.     ct = [];
  26.     for i =1:1:length(sx)
  27.         if sx(i) > Xst + (k-2)*Xstep && sx(i) < Xst + k*Xstep
  28.             ct = [ct;i];
  29.         else
  30.             continue;
  31.         end
  32.     end
  33.     xyuy(k) = max(usytx(ct));
  34. end
  35. para1 = polyfit(xylx,(xyly+xyuy)/2,1);%拟合的轴线应在两边界正中
  36. Cy = @(x) subs(para1(1))*x + subs(para1(2));

  37. %% view([90 0]) 对z 在[-1388,-1358],x为任意值,搜索y的最大最小值,构成直线y-z关系
  38. Zst = -1388; Zed = -1358; Zstep = 1;
  39. n = 30;yzlz = zeros(n,1);yzly = zeros(n,1);
  40. lsytz = yy(idz);
  41. for k = 1:1:n
  42.     yzlz(k) = Zst + (k-0.5)*Zstep;
  43.     ct = [];
  44.     for i =1:1:length(sz)
  45.         if sz(i) > Zst + (k-2)*Zstep && sz(i) < Zst + k*Zstep
  46.             ct = [ct;i];
  47.         else
  48.             continue;
  49.         end
  50.     end
  51.     yzly(k) = min(lsytz(ct));
  52. end

  53. yzuz = zeros(n,1);yzuy = zeros(n,1);
  54. usytz = yy(idz);
  55. for k = 1:1:n
  56.     yzuz(k) = Zst + (k-0.5)*Zstep;
  57.     ct = [];
  58.     for i =1:1:length(sz)
  59.         if sz(i) > Zst + (k-2)*Zstep && sz(i) < Zst + k*Zstep
  60.             ct = [ct;i];
  61.         else
  62.             continue;
  63.         end
  64.     end
  65.     yzuy(k) = max(usytz(ct));
  66. end
  67. para2 = polyfit((yzly+yzuy)/2,yzlz,1);
  68. Czy = @(y) subs(para2(1))*y + subs(para2(2));
  69. %% 粒子群算法求轴线参数
  70. % X的六参量构成空间两点连成的直线,dscalc求原始数据到该直线的距离的方差
  71. dscalc=@(X) var( sqrt( (( (X(1)-X(4))*yy + (X(5)-X(2))*xx + (X(4)*X(2) - X(1)*X(5)) ).^2 + ...
  72.     ( (X(1)-X(4))*zz + (X(6)-X(3))*xx + (X(4)*X(3) - X(1)*X(6)) ).^2 + ...
  73.     ( (X(2)-X(5))*zz + (X(6)-X(3))*yy + (X(5)*X(3) - X(2)*X(6)) ).^2 )./...
  74.     ( (X(1)-X(4))^2 + (X(2)-X(5))^2 + (X(3)-X(6))^2 ) ) );
  75. % 由投影法得到的方程得到两点( t1,Cy(t1),Czy(Cy(t1)) )及( t2,Cy(t2),Czy(Cy(t2)) ),参照它们设置S的上下边界.
  76. S = particleswarm(dscalc,6,[-1;-10;-1500;-10;-1;-2500],[1;500;-1000;1500;1;-1500]);
  77. xa = S(1);ya = S(2);za = S(3);xb = S(4);yb = S(5);zb = S(6);
  78. fprintf('投影法:轴线在x-y平面投影的方程为: y = %f * x + %f \n',para1(1),para1(2));
  79. fprintf('投影法:轴线在y-z平面投影的方程为: z = %f * y + %f \n',para2(1),para2(2));
  80. fprintf('粒子群法:轴线通过空间上两点: (%.4f,%.4f,%.4f) 以及 (%.4f,%.4f,%.4f) \n',xa,ya,za,xb,yb,zb);
  81. fprintf('粒子群法:轴线的两点式方程为: (%.4f - x)/%.4f = (%.4f - y)/%.4f = (%.4f - z)/%.4f \n',xa,(xa-xb),ya,(ya-yb),za,(za-zb));
  82. Xps = (xa+1*(xb-xa)/100):(xb-xa)/100:(xa+13*(xb-xa)/100);
  83. Yps = (ya+1*(yb-ya)/100):(yb-ya)/100:(ya+13*(yb-ya)/100);
  84. Zps = (za+1*(zb-za)/100):(zb-za)/100:(za+13*(zb-za)/100);
  85. t = 20:1:120;
  86. figure(1);
  87. subplot(121)
  88. scatter3(xx,yy,zz,'g.');xlabel('x');ylabel('y');zlabel('z');hold on;plot3(t,Cy(t),Czy(Cy(t)),'bs',Xps,Yps,Zps,'ko-');
  89. title('在x-y平面的投影');legend('实验数据','投影拟合','粒子群拟合','Location','northeast');view(2);
  90. subplot(122)
  91. scatter3(xx,yy,zz,'g.');xlabel('x');ylabel('y');zlabel('z');hold on;plot3(t,Cy(t),Czy(Cy(t)),'bs',Xps,Yps,Zps,'ko-');
  92. title('在y-z平面的投影');legend('实验数据','投影拟合','粒子群拟合','Location','northwest');view([90 0]);
复制代码

新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
 楼主| 发表于 2018-5-23 19:38:39 | 显示全部楼层

非常感谢您的回复和分享,我有些地方不是很明白,为了方便,我直接附上代码在上面标注我的问题了。问题有些多,还望您能耐心解答下:loveliness:
clear all;
data = importdata('圆柱.txt');
xx = data(:,1); yy = data(:,2); zz = data(:,3);
[sx,idx] = sort(xx); [sz,idz] = sort(zz);
%% view(2) 对x 在[40,90],z为任意值,搜索y的最大最小值,构成直线x-y关系
Xst = 40; Xed = 90; Xstep = 1;% 起点,终点,步长
n = 50;xylx = zeros(n,1);xyly = zeros(n,1);
lsytx = yy(idx);%循环找出y的最小值构成的边界
%问题1:for循环里的计算公式是怎么来的。
for k = 1:1:n
    xylx(k) = Xst + (k-0.5)*Xstep;
    ct = [];
    for i =1:1:length(sx)
        if sx(i) > Xst + (k-2)*Xstep && sx(i) < Xst + k*Xstep
            ct = [ct;i];
        else
            continue;
        end
    end
    xyly(k) = min(lsytx(ct));
end
xyux = zeros(n,1);xyuy = zeros(n,1);
usytx = yy(idx);%循环找出y的最大值构成的边界
for k = 1:1:n
    xyux(k) = Xst + (k-0.5)*Xstep;
    ct = [];
    for i =1:1:length(sx)
        if sx(i) > Xst + (k-2)*Xstep && sx(i) < Xst + k*Xstep
            ct = [ct;i];
        else
            continue;
        end
    end
    xyuy(k) = max(usytx(ct));
end
para1 = polyfit(xylx,(xyly+xyuy)/2,1);%拟合的轴线应在两边界正中
Cy = @(x) subs(para1(1))*x + subs(para1(2));
%% view([90 0]) 对z 在[-1388,-1358],x为任意值,搜索y的最大最小值,构成直线y-z关系
Zst = -1388; Zed = -1358; Zstep = 1;
n = 30;yzlz = zeros(n,1);yzly = zeros(n,1);
lsytz = yy(idz);
for k = 1:1:n
    yzlz(k) = Zst + (k-0.5)*Zstep;
    ct = [];
    for i =1:1:length(sz)
        if sz(i) > Zst + (k-2)*Zstep && sz(i) < Zst + k*Zstep
            ct = [ct;i];
        else
            continue;
        end
    end
    yzly(k) = min(lsytz(ct));
end
yzuz = zeros(n,1);yzuy = zeros(n,1);
usytz = yy(idz);
for k = 1:1:n
    yzuz(k) = Zst + (k-0.5)*Zstep;
    ct = [];
    for i =1:1:length(sz)
        if sz(i) > Zst + (k-2)*Zstep && sz(i) < Zst + k*Zstep
            ct = [ct;i];
        else
            continue;
        end
    end
    yzuy(k) = max(usytz(ct));
end
para2 = polyfit((yzly+yzuy)/2,yzlz,1);
Czy = @(y) subs(para2(1))*y + subs(para2(2));
%% 粒子群算法求轴线参数
% X的六参量构成空间两点连成的直线,dscalc求原始数据到该直线的距离的方差
dscalc=@(X) var( sqrt( (( (X(1)-X(4))*yy + (X(5)-X(2))*xx + (X(4)*X(2) - X(1)*X(5)) ).^2 + ...
    ( (X(1)-X(4))*zz + (X(6)-X(3))*xx + (X(4)*X(3) - X(1)*X(6)) ).^2 + ...
    ( (X(2)-X(5))*zz + (X(6)-X(3))*yy + (X(5)*X(3) - X(2)*X(6)) ).^2 )./...
    ( (X(1)-X(4))^2 + (X(2)-X(5))^2 + (X(3)-X(6))^2 ) ) );

%下面这句话不懂,投影法得到的两点怎么解释?是自己随便取的两点吗还是怎么?
% 由投影法得到的方程得到两点( t1,Cy(t1),Czy(Cy(t1)) )及( t2,Cy(t2),Czy(Cy(t2)) ),参照它们设置S的上下边界.
S = particleswarm(dscalc,6,[-1;-10;-1500;-10;-1;-2500],[1;500;-1000;1500;1;-1500]);%S上下边界,参数的范围是怎么确定的?如果给定的点云数据变了这里的边界又该如何改?
xa = S(1);ya = S(2);za = S(3);xb = S(4);yb = S(5);zb = S(6);
fprintf('投影法:轴线在x-y平面投影的方程为: y = %f * x + %f \n',para1(1),para1(2));
fprintf('投影法:轴线在y-z平面投影的方程为: z = %f * y + %f \n',para2(1),para2(2));
fprintf('粒子群法:轴线通过空间上两点: (%.4f,%.4f,%.4f) 以及 (%.4f,%.4f,%.4f) \n',xa,ya,za,xb,yb,zb);
fprintf('粒子群法:轴线的两点式方程为: (%.4f - x)/%.4f = (%.4f - y)/%.4f = (%.4f - z)/%.4f \n',xa,(xa-xb),ya,(ya-yb),za,(za-zb));
Xps = (xa+1*(xb-xa)/100):(xb-xa)/100:(xa+13*(xb-xa)/100);
Yps = (ya+1*(yb-ya)/100):(yb-ya)/100:(ya+13*(yb-ya)/100); %这里的公式又是怎么来的呢?为什么只用13个点就可以呢?
Zps = (za+1*(zb-za)/100):(zb-za)/100:(za+13*(zb-za)/100);

t = 20:1:120;%这里的t是什么含义呢?为什么取在20-120之间,换个区间,轴线的位置好像会变呢
figure(1);
subplot(121)
scatter3(xx,yy,zz,'g.');xlabel('x');ylabel('y');zlabel('z');hold on;plot3(t,Cy(t),Czy(Cy(t)),'bs',Xps,Yps,Zps,'ko-');
title('在x-y平面的投影');legend('实验数据','投影拟合','粒子群拟合','Location','northeast');view(2);
subplot(122)
scatter3(xx,yy,zz,'g.');xlabel('x');ylabel('y');zlabel('z');hold on;plot3(t,Cy(t),Czy(Cy(t)),'bs',Xps,Yps,Zps,'ko-');
title('在y-z平面的投影');legend('实验数据','投影拟合','粒子群拟合','Location','northwest');view([90 0]);

论坛优秀回答者

入门

435 麦片

财富积分


50500


2

主题

432

帖子

95

最佳答案
  • 关注者: 24
发表于 2018-5-24 02:40:47 | 显示全部楼层 |此回复为最佳答案
云妹y 发表于 2018-5-23 07:38
非常感谢您的回复和分享,我有些地方不是很明白,为了方便,我直接附上代码在上面标注我的问题了。问题有 ...

第一个问题
[sx,idx] = sort(xx);%实验数据是随机的,先对x轴数据按照大小进行排序
%% 此前观察散点图可知 对x 在[40,90]区间,散点图的边界比较明显容易判别,因此可取z为任意值,x 在[40,90],搜索y的最大最小值,构成直线x-y关系
Xst = 40; Xed = 90; Xstep = 1;% x值的起点,终点,步长
n = 50;xylx = zeros(n,1);xyly = zeros(n,1);%n表示执行50次,n=(Xst- Xed)/Xstep,可以令deltaX = abs(sx(1:(end-1))-sx(2:end)) 即x轴相邻两数据间的差值,Xstep的设置可参考差值的最大值,以保证每一步的区间内都有数据点,若Xstep太小,可能在一些循环的步骤中找不到数据点.
for k = 1:1:n
    xylx(k) = Xst + (k-0.5)*Xstep;%搜集位于边界的数据,xylx(k)为其x值,每一步搜集区间是以 xylx(k) 为参考中心的一个区间,不是必须如此设置,只是参考
    ct = [];% ct储存位于区间之内的数据点对应的序号,每一步搜索时都设为空
    for i =1:1:length(sx)
        if sx(i) > Xst + (k-2)*Xstep && sx(i) < Xst + k*Xstep % 判断数据点是否落在搜集区间里面,此处搜集区间为 [参考中心左侧1.5个Xstep,参考中心右侧0.5个Xstep],这个跨度保证区间数据数量不为0
            ct = [ct;i]; % 位于区间之内的数据点对应的序号
        else
            continue;
        end
    end
    xyly(k) = min(lsytx(ct)); % 根据序号,对位于区间之内的数据求最小值,即y值的下界,得到区间最小y值数据点构成的集合(xylx,xyly).注意,xyly是该区间最小y值数据点的y坐标,但xylx的取值是人为规定的"区间中心"的值,并不是该区间最小y值数据点的x坐标,这样拟合轴线更方便,但也因此引入了误差.
end
% 之后类似地处理,得到相同区间最大y值数据点构成的集合(xyux,xyuy).其中xyux=xylx,因此区间内y的最大最小值连线是平行于y轴的,轴线必经过这些平行线的中点,拟合xylx与(xyly+xyuy)/2即可. 如果之前xylx的不取"区间中心"的值,而取区间最小y值数据点的x坐标,则区间内y最大最小值的连线不再一定平行于y轴,轴线也不再一定经过其中点了.即,无论采用区间中心值,还是y最值对应的x坐标,都会有误差,这里就采用更方便的前者了.

第二个问题
两点不是随便找的,空间上任一平面可表示为A*x+B*y+C*z+D=0,那么轴线在x-y和y-z两个面上的投影平面
Cy = @(x) subs(para1(1))*x + subs(para1(2));
Czy = @(y) subs(para2(1))*y + subs(para2(2));
联立,就是面交线,即轴线,的方程.两个方程三个自变量,所以解( t,Cy(t),Czy(Cy(t)) )里面会带有一个参量t,这个t可以是任意实数值.
为了直观,我选取的是,令x坐标近于0的t1,以及令y坐标近于0的t2,这样得到的两点分别为(0,424.5913,-1320.8752)以及(983.3593, 0, -2036.7256),以它们为参考值,设置S的边界为同数量级且包含两点的区间即可,因为目的只是提高参数"已知"的轴线方程的精度而不是寻找"未知"的参数,区间不必太大,这样可缩短运算时间.另,求出的S值并非唯一,因为轴线上有无数点,S求的只是其中距离初值较近的任意两点.

第三个问题
Xps = (xa+1*(xb-xa)/100):(xb-xa)/100:(xa+13*(xb-xa)/100);
Yps = (ya+1*(yb-ya)/100):(yb-ya)/100:(ya+13*(yb-ya)/100);
Zps = (za+1*(zb-za)/100):(zb-za)/100:(za+13*(zb-za)/100);
t = 20:1:120;
轴线的关键是方程本身,不是参量.标粗的分别与起点,步长,终点有关,都是可以自己调整的参量.
比如x-y平面上的直线y=2x+1,参数形式x=t,y=2t+1,直线本身是无限长的,但为了直观做图通常只取它的一部分线段,比如t1=10:20或者t2=-100:-90,它们是同一条直线的不同部分,但做出图来,t1参数的图在第一象限,t2做出来的图在第三象限,不看整体只看局部的话,会产生这两者完全不同的错觉.
修改标粗的那些参数只是改变轴线所显示出来的线段的起点和终点而已,类似于把"申"改成"甲" "由" "田".


新手

16 麦片

财富积分


050


3

主题

38

帖子

0

最佳答案
 楼主| 发表于 2018-5-29 08:36:40 | 显示全部楼层
TouAkira 发表于 2018-5-24 02:40
第一个问题
[sx,idx] = sort(xx);%实验数据是随机的,先对x轴数据按照大小进行排序
%% 此前观察散点图可知 ...

嗯嗯,真的非常感谢您的慷慨解答!!!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表