[已解决] DIC获取的部分圆柱面将它拟合成圆柱,求其轴线

[复制链接]
之乎者也1 发表于 2022-10-21 15:24:45
球球大家了,求大家给个方法吧?:'(

matlab-2-2-dis.mat

54.31 KB, 下载次数: 3

圆柱面数据

最佳答案


TouAkira 发表于 2022-10-21 18:20:24
原理讲了好多遍了,照着写然后算就是了
拟合出的轴线端点参考值
  1. 6.99248        -3.55689         520.423
  2. 15.4961        -4.28875         517.894
复制代码

[attach]333761[/attach]

回复此楼

28 条回复


TouAkira 发表于 2022-10-21 18:20:24
原理讲了好多遍了,照着写然后算就是了
拟合出的轴线端点参考值
  1. 6.99248        -3.55689         520.423
  2. 15.4961        -4.28875         517.894
复制代码

CylinderFitting2022Oct21_1.png

回复此楼

之乎者也1 发表于 2022-10-21 18:56:50
TouAkira 发表于 2022-10-21 18:20
原理讲了好多遍了,照着写然后算就是了
拟合出的轴线端点参考值

老师,我没找到这个计算原理,方便分享一下嘛?谢谢老师

TouAkira 发表于 2022-10-21 20:06:05
本帖最后由 TouAkira 于 2022-10-21 08:07 编辑
之乎者也1 发表于 2022-10-21 06:56
老师,我没找到这个计算原理,方便分享一下嘛?谢谢老师

不方便
原理我讲过很多遍了
旋转体拟合轴线都是有套路的,去参考我之前的回复
提取点云圆柱的中心轴线
不完整圆柱点云提取中轴线并拟合圆柱
提取圆锥中轴线

当时我把代码发出来,甚至每一步都仔细注释为什么要这样算,结果原帖楼主别说确认为最佳答案了,连个最起码的道谢都没有。前一阵就有个论坛用户,他提问当天我费劲写了近千字回帖仔细分析讲解,结果他看都不看,等过了一个月后又想起来了,就把问题复制粘贴重问一遍。
类似这样的论坛用户,太多太多了。退役版主 kastin 说得对
kastin 发表于 2017-3-19 06:22:51
免费廉价往往被人们轻视

既然免费的好意得不到最起码的尊重,那我也不会再直接给出代码了。原理我都讲过不止一遍,或者自己按照原理写代码,或者私下咨询。不要打着“分享”旗号要求别人把东西给你,除非先把你自己的东西拿出来“分享”。

之乎者也1 发表于 2022-10-21 20:32:13
TouAkira 发表于 2022-10-21 20:06
不方便。
原理我讲过很多遍了
当时我把代码发出来,甚至每一步都仔细注释为什么要这样算,结果原帖楼主别 ...

ptCloud2=pointCloud(matlab_dis_2);
figure
pcshow(ptCloud2)
title('Original Point Cloud')

maxDistance = 15;%点到圆柱面的距离阈值
roi = [-inf,-inf,-inf,inf,inf,inf];%点云ROI选取。不需要ROI则设置为-inf到inf就可以了
sampleIndices = findPointsInROI(ptCloud2,roi);%ROI点索引
referenceVector = [1,0,0];%初始轴向。根据情况,设置为xyz轴都可以
[model,inlierIndices] = pcfitcylinder(ptCloud2,maxDistance,...
        referenceVector,'SampleIndices',sampleIndices);
pc = select(ptCloud2,inlierIndices);
figure
pcshow(pc)
title('Cylinder Point Cloud')
hold on
plot(model)
axis equal
zlim([507 511])

老师,真的很抱歉,您的感受我能理解,我开始认为您不分享代码的是因为秉着授人以鱼不如授人以渔的想法。我虽然不是什么圣人,但是我会很感谢帮助过我的人。我也觉得这是一个论坛,一个讨论的平台,把自己的想法分享出来,大家一起讨论。我一开始有问题,并没有把自己拟合的方法表达,当我想向老师私下请教的时候,网页提示我还只是一个新手,新手不能够私发信息给您,所以在这边提供我写的代码,希望老师能够帮忙分析一下,上述是我用matlab自带的圆柱拟合代码拟合,但是拟合的结果并没有老师您拟合的结果好,所以在此也想真诚希望老师能够给我一份代码,如果您还是不方便,我也同样非常感谢您的原理提示。十分感谢!希望老师以后心情都能愉悦~

之乎者也1 发表于 2022-10-21 21:45:41
TouAkira 发表于 2022-10-21 20:06
不方便。
原理我讲过很多遍了
当时我把代码发出来,甚至每一步都仔细注释为什么要这样算,结果原帖楼主别 ...

老师,我有个问题,向您请教一下,能不能指点一下呢?就是我现在只是知道这个圆柱面的点云坐标,但是中心轴线是不知道的?我们怎么实现点云到轴线的距离是一样的呢?

顾世梁 发表于 2022-10-22 16:36:57
fx=b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
SSy =   3014.8
b =[522.6881326  -0.345238879  -0.01403124636  0.008061419546  142.3318403  -1.788796121  -7.083479003  -1.039043508];
RSS = 4.82345524398
MSe =  0.0018431
R^2 =    0.9984
nh1771.png

之乎者也1 发表于 2022-10-22 18:05:17
顾世梁 发表于 2022-10-22 16:36
fx=b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
SSy =   3014.8
b =[522.6881326   ...

拟合结果挺好的,但是我实现不了:'(
辛苦顾老师了
可以出代码嘛?:'(

顾世梁 发表于 2022-10-22 19:02:32
用lsqcurvefit and/or nlinfit, cftool 等,以我给出的模型和初值(以我的结果作为初值),即可实现该组数据的拟合。

顾世梁 发表于 2022-10-22 19:17:49
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
fx2=@(b,X)b(1)+b(2)*X(:,1)+b(3)*X(:,2)+b(4)*X(:,1).*X(:,2)-sqrt(b(5)+b(6)*X(:,1)+b(7)*X(:,2)+b(8)*X(:,2).^2);
x1=[...]; %输入数据
x2=[...]; %
y=[...]; %
X=[x1,x2];
b=[522 -0.34524  -0.014  0.00806  142  -1.7888  -7.08  -1.039];
for l=1:5
    b=lsqcurvefit(fx2,b,X,y);
    b=nlinfit(X,y,fx2,b);
end
figure(1),clf
plot3(x1,x2,y,'o','markersize',5,'markerfacecolor','k')
hold on
[x11,x22]=meshgrid(min(x1):range(x1)/101:max(x1),min(x2):range(x2)/101:max(x2));
y1=fx1(b,x11,x22);
surf(x11,x22,y1)
alpha(.8)
axis tight
yhat=fx1(b,x1,x2);
b
RSS=(y-yhat)'*(y-yhat)
n=length(y); SSy=var(y)*(n-1);
R2=(SSy-RSS)/SSy

之乎者也1 发表于 2022-10-22 21:12:16
顾世梁 发表于 2022-10-22 19:17
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
fx2=@(b,X)b(1 ...

谢谢顾老师,您的代码能给出来,真是犹如久旱逢甘露,十分感谢!我先学习您的代码,等我有不懂的地方再向您请教可以嘛?我有一个问题,如果我们一开始没有初值,该怎么去拟合这组数据呀?

顾世梁 发表于 2022-10-22 22:10:22
我自编了拟合的程序,可在随机初值的情况下实现一般问题的最优拟合。

之乎者也1 发表于 2022-10-22 22:30:41
顾世梁 发表于 2022-10-22 22:10
我自编了拟合的程序,可在随机初值的情况下实现一般问题的最优拟合。

十分感谢顾老师能够回复!顾老师真厉害!我想更深入的拟合出整个圆柱面,拟合出其轴线,该用matlab怎么实现呢?顾老师有没有推荐的专业书籍呢?因为我觉得我现在在这方面的知识比较欠缺,需要恶补一下。

顾世梁 发表于 2022-10-22 23:21:07
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
fx2=@(b,X)b(1)+b(2)*X(:,1)+b(3)*X(:,2)+b(4)*X(:,1).*X(:,2)-sqrt(b(5)+b(6)*X(:,1)+b(7)*X(:,2)+b(8)*X(:,2).^2);
fx=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y+sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
x1=xlsread('nh1771','a1:a2626');
x2=xlsread('nh1771','b1:b2626');
y=xlsread('nh1771','c1:c2626');
X=[x1,x2];
b=[522 -0.34524  -0.014  0.00806  142  -1.7888  -7.08  -1.039];
for l=1:5
    b=lsqcurvefit(fx2,b,X,y);
    b=nlinfit(X,y,fx2,b);
end
figure(1),clf
plot3(x1,x2,y,'o','markersize',5,'markerfacecolor','k')
hold on
r1=range(x1); r2=range(x2);
[x11,x22]=meshgrid(min(x1)-r1/10:r1/150:max(x1)+r1/10,min(x2)-r2/2:r2/150:max(x2)+r2/2);
y1=fx1(b,x11,x22); y2=fx(b,x11,x22);
[n1,n2]=size(y1);
for i=1:n1
    for j=1:n2
        if ~isreal(y1(i,j))
            y1(i,j)=nan;
        end
        if ~isreal(y2(i,j))
            y2(i,j)=nan;
        end
    end
end
surf(x11,x22,y1), surf(x11,x22,y2)
alpha(.68)
axis tight
yhat=fx1(b,x1,x2);
b
RSS=(y-yhat)'*(y-yhat)
n=length(y); SSy=var(y)*(n-1);
R2=(SSy-RSS)/SSy

顾世梁 发表于 2022-10-23 06:56:37
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
fx2=@(b,X)b(1)+b(2)*X(:,1)+b(3)*X(:,2)+b(4)*X(:,1).*X(:,2)-sqrt(b(5)+b(6)*X(:,1)+b(7)*X(:,2)+b(8)*X(:,2).^2);
fx=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y+sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
X=[...]; %输入数据
x1=X(:,1); x2=X(:,2); y=X(:,3);
X=[x1,x2];
b=[522 -0.34524  -0.014  0.00806  142  -1.7888  -7.08  -1.039];
for l=1:5
    b=lsqcurvefit(fx2,b,X,y);
    b=nlinfit(X,y,fx2,b);
end
figure(1),clf
plot3(x1,x2,y,'o','markersize',5,'markerfacecolor','k')
hold on
r1=range(x1); r2=range(x2);
[x11,x22]=meshgrid(min(x1)-r1/10:r1/180:max(x1)+r1/10,min(x2)-r2/2:r2/180:max(x2)+r2/2);
y1=fx1(b,x11,x22); y2=fx(b,x11,x22);
[n1,n2]=size(y1);
for i=1:n1
    for j=1:n2
        if ~isreal(y1(i,j))
            y1(i,j)=nan;
        end
        if ~isreal(y2(i,j))
            y2(i,j)=nan;
        end
    end
end
surf(x11,x22,y1), surf(x11,x22,y2)
shading interp
alpha(.8)
axis tight
yhat=fx1(b,x1,x2);
b
RSS=(y-yhat)'*(y-yhat)
n=length(y); SSy=var(y)*(n-1);
Rsquare=(SSy-RSS)/SSy

之乎者也1 发表于 2022-10-23 16:45:02
顾世梁 发表于 2022-10-23 06:56
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
fx2=@(b,X)b(1 ...

很感谢顾老师能够回复!让小生受益匪浅,我拟合了之后发现,圆柱并没有完全贴合,中心轴线的拟合我还是没有实现,还请顾老师不吝赐教:'(

之乎者也1 发表于 2022-10-23 22:22:48
TouAkira 发表于 2022-10-21 20:06
不方便。
原理我讲过很多遍了
当时我把代码发出来,甚至每一步都仔细注释为什么要这样算,结果原帖楼主别 ...

YouAkira老师您好!我还是没有拟合出轴线,老师能提供您的代码学习嘛?实在是不会拟合:'(:'(,十分感谢!

顾世梁 发表于 2022-10-24 14:54:44
之乎者也1 发表于 2022-10-23 16:45
很感谢顾老师能够回复!让小生受益匪浅,我拟合了之后发现,圆柱并没有完全贴合,中心轴线的拟合我还是没 ...

fx1的前三项即中心平面,取x2的中值即为中心线。

之乎者也1 发表于 2022-10-24 23:04:07
十分感谢顾老师的指教!我能够实现轴线,但是对代码还有一些疑惑不知顾老师能不能再指导一下?
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y+b(4)*x.*y-sqrt(b(5)+b(6)*x+b(7)*y+b(8)*y.^2);
上面式子b这个参数有一些疑惑,b代表的参数是圆柱的7个未知参数:轴线上的点坐标(x0,y0,z0),轴线方向向量(l,m,n),半径r嘛?但是b怎么是8个参数呢?

您理解的圆柱面方程是这个方程嘛?感谢顾老师指点!

clear,clc
fx1=@(b,x,y,z)(x-b(1)).^2+(y-b(2)).^2+(z-b(3)).^2-b(7).^2-(b(4)*(x-b(1))+b(5)*(y-b(2))+b(6)*(z-b(3))).^2;
%b是个我理解的圆柱面方程参数,共7个,但是程序跑不起来
load('matlab-2-2-dis');
X=matlab_dis_2; % 输入数据
x=X(:,1); y=X(:,2); z=X(:,3);
X=[x,y];
figure(1),clf
plot3(x,y,z,'o','markersize',5,'markerfacecolor','k')
hold on

b=[500 -0.34524  -0.014  0.00806  142  -1.7888  -7.08  -1.039];%拟合所求方程系数初值
for l=1:5
    b=lsqcurvefit(fx1,b,X,z);
%     b=nlinfit(X,z,fx2,b);
end
r1=range(x); r2=range(y);
[x11,x22]=meshgrid(min(x)-r1/10:r1/180:max(x)+r1/10,min(y)-r2/2:r2/180:max(y)+r2/2);
y1=fx1(b,x11,x22); %y2=fx(b,x11,x22);
%去掉y中非实数
[n1,n2]=size(y1);
for i=1:n1
    for j=1:n2
        if ~isreal(y1(i,j))
            y1(i,j)=nan;
        end
%         if ~isreal(y2(i,j))
%             y2(i,j)=nan;
%         end
    end
end
%画图
surf(x11,x22,y1);
shading interp
alpha(.8)
axis tight
axis equal
%求系数结果b及拟合误差
yhat=fx1(b,x,y);
b   
RSS=(z-yhat)'*(z-yhat)%离差平方和:各项与平均项之差的平方的总和
n=length(z); SSy=var(z)*(n-1);
Rsquare=(SSy-RSS)/SSy


顾世梁 发表于 2022-10-24 23:32:13
主要有b(4)*x.*y这一项,该组数据并非标准的圆柱,而是稍有扭曲。当加入该项时,RSS得到有效减小。你若不想要该项,我可重新拟合成没有该项的圆柱方程。稍后给你答案。

顾世梁 发表于 2022-10-24 23:57:25
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y-sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y.^2);
b =[575.3153892  -1.129110086  -0.3639628717  4094.816527  -96.29183048  -85.59047694  -5.421967143];

之乎者也1 发表于 2022-10-25 21:27:57
顾世梁 发表于 2022-10-24 23:57
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y-sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y.^2);
b =[575.3153892  -1.129110086  -0 ...

感谢顾老师回复!我修改了您的代码吼,未能实现复现,顾老师能给出您的代码吗?感谢~
以下是我修改的代码,希望顾老师能指点迷津。
clear,clc
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y-sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y.^2);
sqrt(b(5)+b(6)*X(:,1)+b(7)*X(:,2)+b(8)*X(:,2).^2);
fx4=@(b,x,y)b(1)+b(2)*x+b(3)*y;
load('matlab_2_2_dis');
X=matlab_dis_2; % 输入数据
x=X(:,1); y=X(:,2); z=X(:,3);
X=[x,y];
figure(1),clf
plot3(x,y,z,'o','markersize',5,'markerfacecolor','k')
hold on

b =[575.3153892  -1.129110086  -0.3639628717  4094.816527  -96.29183048  -85.59047694  -5.421967143];
for l=1:5
    b=lsqcurvefit(fx1,b,X,z);
end
r1=range(x); r2=range(y);
[x11,x22]=meshgrid(min(x)-r1/10:r1/180:max(x)+r1/10,min(y)-r2/2:r2/180:max(y)+r2/2);
y1=fx1(b,x11,x22);
y3=fx4(b,x11,median(x22));
[n1,n2]=size(y1);
for i=1:n1
    for j=1:n2
        if ~isreal(y1(i,j))
            y1(i,j)=nan;
        end
%         if ~isreal(y2(i,j))
%             y2(i,j)=nan;
%         end
    end
end
surf(x11,x22,y1);
% surf(x11,x22,y2);
% plot3(x11,median(x22),y3,'.','markersize',5);
% plot3([6.99248 15.4961],[-3.55689 -4.28875],[520.423 517.894],'*');
shading interp
alpha(.8)
axis tight
yhat=fx1(b,x,y);
b   
RSS=(z-yhat)'*(z-yhat)
n=length(z); SSy=var(z)*(n-1);
Rsquare=(SSy-RSS)/SSy





顾世梁 发表于 2022-10-25 21:42:39
错在代码的第3~4两行,3~5行请照此修改:
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y-sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y.^2);
fx2=@(b,X)b(1)+b(2)*X(:,1)+b(3)*X(:,2)-sqrt(b(4)+b(5)*X(:,1)+b(6)*X(:,2)+b(7)*X(:,2).^2);
fx=@(b,x,y)b(1)+b(2)*x+b(3)*y+sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y.^2);

之乎者也1 发表于 2022-10-25 22:11:36
顾世梁 发表于 2022-10-25 21:42
错在代码的第3~4两行,3~5行请照此修改:
fx1=@(b,x,y)b(1)+b(2)*x+b(3)*y-sqrt(b(4)+b(5)*x+b(6)*y+b(7)*y ...

感谢顾老师的及时回复!学生对这个参数b还是有疑惑,我们拟合出的圆柱,圆柱的参数是通过b知道的吗?还是我们在做的只是将圆柱形状拟合出来,圆柱的参数未可知?因为我对b这个参数有比较大的疑惑,希望顾老师能指点一下。十分感谢!

顾世梁 发表于 2022-10-26 06:10:08
圆柱的方程即为fx1(下), fx(上)两面。参数b决定了圆柱表面的具体形态。根据数据拟合方程即为估计参数b,怎能说圆柱的参数未知?对于一般非线性拟合问题,不存在公式解求算b,只能通过给定初值的条件下不断迭代调整,得到b的优化估计。matlab自带的拟合命令,对于比较简单的问题,可在随机初值条件下达到最优解,而对于比较复杂的方程和数据,优化估计的能力较弱,需提供合适的初值才能有较好的效果。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表