已知曲面点云数据，如何利用求交法计算交线轮廓点数据？算法原理我截图放在下面了曲面点云坐标我存放在一个"p. txt"文件中，也有Excel版的。

Question

不知道能否在MATLAB中编程实现？曲面点云数据太多，这里放不下，我只截取了少部分：
352.479248 89.48173523 162.313385
352.4787903 89.54194641 159.8517914
346.180603 88.32884216 159.8513947
352.4784546 89.58816528 157.3894806
352.4781799 89.62073517 154.9266205
346.1799622 88.40332794 154.9264069
352.481842 89.09365845 172.1489716
352.482666 88.95804596 174.6042633
358.782135 90.03437042 174.605011
387.1472778 92.69024658 181.9615936
387.1471558 92.91171265 179.5123901
383.9943237 92.69641876 179.5122681
339.887146 86.61011505 172.147522
339.8863525 86.71890259 169.6906891
333.5910339 85.32966614 169.6900635
321.0059204 81.862854 179.5063934
314.7128296 80.27714539 179.5062103
314.7137756 79.98086548 184.4034729
232.9160309 57.8214035 169.7000427
239.210556 59.49936676 169.6986084
239.212677 59.64931107 159.854599
289.5455322 73.65122986 174.6020203
289.5449829 73.45576477 179.5079803
302.1282959 76.94490051 179.5065613
333.5994263 83.85174561 191.7304382
333.5983276 84.08112335 189.2910004
327.3031311 82.66703796 189.2909241
327.3041077 82.44667053 191.730545
327.3050842 82.20880127 194.1673126
333.6005249 83.60454559 194.167038
327.3060608 81.95292664 196.6011047
321.0102234 80.73815155 194.1677704
308.4220886 77.60447693 194.169281
314.7159119 79.20067596 194.1684113
314.7148132 79.6234436 189.2913055
195.1248932 49.34513855 159.8592224
188.821106 48.2123642 159.8596802
188.8216553 48.10380936 164.7861481
122.6195831 43.46050644 150
119.470108 43.58943176 150
119.4702148 43.57605743 152.4655457
103.733139 44.06523895 167.2483978
100.5865555 44.40322113 167.2480774
100.5866089 44.18263245 169.7088165
141.5358276 42.5806694 172.1712036
138.3838348 42.47247314 172.1711731
138.3859253 42.27077866 174.629425
141.5379181 42.38530731 174.6294556
138.3881531 42.0486908 177.0860138
141.5476074 41.40615463 184.4442139
138.3956604 41.26166534 184.4442444
138.3983765 40.95960236 186.8925781
128.940567 41.03086472 184.443985
125.7891769 41.02352142 184.4437866
125.7941895 40.33057785 189.338089
103.7357407 42.37817383 181.9903564
103.7352295 42.7228775 179.5377197
100.5868988 43.03321075 179.5371552
113.1909027 40.18465424 191.7798767
113.1926727 39.750103 194.2200623
116.3442535 39.64498138 194.2205505
122.6476669 39.54827118 194.2212372
119.4959183 39.57793427 194.220932
119.5006332 38.68481445 199.0909729
147.8713531 39.46974564 201.5197906
144.7197113 39.195961 201.5201416
144.7227173 38.78204346 203.9449158
106.893013 38.59479141 201.5180664
106.8940811 38.06061935 203.9430389
110.0470657 37.91475296 203.9437103
100.5877533 39.02837753 201.5164185
100.5878601 38.47539902 203.941391
103.7410202 38.2472496 203.942276
157.3457031 36.99240494 220.7718811
157.3479767 36.48025131 223.1518555
163.6473694 37.46087265 223.1515656
169.9473419 38.04698563 225.525116
169.9488831 37.53359604 227.8920135
176.245636 38.7388649 227.8927002
141.5830383 36.65897751 213.5954437
138.4312286 36.39822769 213.5957184
138.4341583 35.87759018 215.9944458
94.28246307 39.63037491 201.514328
94.28136444 39.05936813 203.9392853
97.43462372 38.74583435 203.9403992
87.97750092 40.40711594 201.511795
87.97502136 39.81901169 203.9367371
91.12814331 39.41682053 203.9380646
75.37038422 42.51052094 201.5053406
75.36490631 41.89172363 203.9301605
78.51695251 41.30230713 203.9319916
69.06915283 43.85003281 201.501358
69.06219482 43.2179718 203.9260864
72.21330261 42.53000259 203.9281921
62.77009201 45.38999939 201.4968567
62.76174545 44.74627304 203.9214783
65.91165161 43.95645523 203.9238586
56.47364807 47.13683319 201.491806
56.46408844 46.48322296 203.9163055
59.61252975 45.58825684 203.9189606
43.8904686 51.27682877 201.4800262
43.87927246 50.60943222 203.9041748
47.02400208 49.49427795 203.9074249
37.60465622 53.68284607 201.4732361
37.59317398 53.01187897 203.897171
40.73564148 51.7816925 203.9007568
31.32330322 56.32143402 201.4658508
31.31198502 55.64934921 203.8895264
34.45193481 54.30082703 203.8934326
25.04687309 59.19880676 201.4577942
25.0362606 58.52814865 203.8812256
28.17340279 57.05821228 203.8854675
21.91065788 60.72885513 201.453537
21.90063095 60.05987549 203.8768158
352.5108948 79.63265228 233.8579559
352.5113831 79.29969788 235.0427246
358.8077698 80.13591766 235.0449524
352.5127563 78.26118469 238.5885162
352.5131531 77.90163422 239.7674408
358.8093872 78.71581268 239.7710724
349.3620605 79.49656677 232.6709595
346.2135925 79.00943756 232.6700745
346.2140808 78.68875885 233.8558807
384.0015564 75.19369507 254.967804
384.0016174 74.71365356 256.1222229
387.1521301 74.80090332 256.1239014
251.7717743 57.21540833 230.2720642
251.7727051 57.67344666 227.9024811
239.186142 54.19793701 227.9008026
270.6485596 54.62760544 258.2266541
264.3540039 52.99527359 258.2139893
264.3540649 52.16174316 260.4944458
270.6486206 53.78556824 260.5078125
264.3542175 51.29572678 262.7668152
346.218811 63.55766296 273.0898438
346.2184753 62.97745514 274.2033386
352.5154114 63.75004959 274.2161865
390.303894 54.81349945 292.7046204
390.3038635 54.07771683 293.7581482
396.6088257 54.04807281 293.7590942
396.6088562 53.30640411 294.8088074
390.3038635 53.33418274 294.8080139
390.303833 52.58288193 295.8540955
396.6088867 52.55715561 295.8547363
396.6089172 51.80031586 296.8968811
390.303833 51.82380676 296.8963318
390.3038025 51.05693054 297.9347839
396.6089172 51.0358696 297.9351196
396.6089478 50.26380539 298.9695435
358.7970581 49.40677643 297.9284973
365.096405 49.9527092 297.9302979
365.0980225 51.48841476 295.8455811
258.0837097 29.55899429 300
264.3756714 31.01039314 300
264.37323 32.52234268 297.8904724
207.7228851 23.21099091 291.4196167
214.0186615 24.67689133 291.4282227
214.0174866 25.9001236 289.2636719
245.4885101 32.25561142 291.4710388
239.1953888 30.72506523 291.4624634
239.1976624 29.39570236 293.6114197
195.1315765 42.34568405 230.2576904
195.1317596 41.84274292 232.6136475
201.4254303 43.35517883 232.6160583
122.6938934 26.24068642 246.5496979
125.8486786 26.46017265 246.5495148
125.8465805 27.94641113 241.9267578
195.1328125 35.60373688 255.7967529
188.84021 34.09359741 255.7851715
188.8401642 33.34037399 258.0643311
169.9585419 29.15918922 258.0297852
169.9582214 28.37870216 260.2999268
176.2530212 29.70764732 260.3114929
125.849411 25.69800377 248.8505859
125.8502045 24.13677979 253.4323883
125.8499146 24.9234581 251.1448517
122.6950684 24.68991852 251.1450653
103.7495346 25.7221241 246.5587311
100.5892792 25.77842522 246.5617065
100.5893173 24.12003708 251.1598511
103.7498093 24.08002472 251.1561584
100.589325 22.42566872 255.7307434
78.45560455 29.19773293 241.9620819
75.29368591 29.65155029 241.967392
75.28957367 27.88884354 246.603241
72.12487793 26.56020164 251.2207031
72.12763214 28.376688 246.6108551
68.96640778 28.91799736 246.6190796
68.96877289 29.82253265 244.3026733
67.38613892 29.20904732 246.6233978
3.11E-13 53.06700516 254.0130005
2.95E-13 51.8722496 256.3399658
3.091678619 49.94054031 256.3035583
2.77E-13 50.65631104 258.6569214
2.59E-13 49.41960144 260.9632874
3.09223628 47.51965714 260.9221191
188.8393097 30.08695602 267.1104126
188.8389282 29.21190834 269.3541565
201.4252777 32.24407959 269.3840637

Answer 1

同学，你好，看了你这个问题和步骤，我认为是很好实现的，但是我需要你的具体数据，点云实在是太少了。我只能根据你给的数据大致画一个曲面点云

另外，我需要确认的几条：
（1）PLHk的方程里ABCD四个参数是不是自己给定的，还是说我拟合出一个平面（我觉得切平面应该是局部域内拟合）
（2）阈值ε需要自己给定吧，所以ε代表的是局部大小的一个值吗？

代码如下，你如果自己求出了ABCD并且定义了epsil就可以输出相应的交点坐标集合M，M的每行都是一个交点，可以遍历求得其法向

function [P, P_Pairs, M] = getPairsCross(A,B,C,D,epsil)
%% 输入
% A B C D 切平面参数方程
% epsil切平面两边的厚度
%% 输出
% P 所有的点云的点
% P_Pairs 互为最近的点的编号
% M 互为最近的点的交点坐标
%% 文件夹【需要自己指定】
folder = '点云切片求交资料/';%储存的点云文件夹
P = load([folder,'曲面点云坐标.txt']); % 加载点云坐标
if(nargin<4)
    epsil = 2; % ε的大小
end
if(nargin==0)
    figure(1);clf
    scatter3(P(:,1),P(:,2),P(:,3),10,'r','filled')
    fP = 100;% 随便选了个点
    dx = 10; % 局部距离
    % 求局部的面
    D = sqrt(sum((P - P(fP,:)).^2,2))<dx;
    eP = P(D,:);
    QPM = @(a,x,y,z) a(1)*x + a(2)*y + a(3)*z + 1;%切面方程
    a = (eP'*eP)\(eP'*(-ones(size(eP(:,1))))); % 得到面的参数
else
    a = [A,B,C]/(D+eps);
end


hold on
scatter3(P(fP,1),P(fP,2),P(fP,3),10,'b','filled');
hold on
fimplicit3(@(x,y,z)QPM(a,x,y,z)) % 画切面

D_P = abs(QPM(a,P(:,1),P(:,2),P(:,3)))/sqrt(sum(a.*a))<epsil;
% 一侧的点
PP = find(D_P);
P_left = PP(QPM(a,P(PP,1),P(PP,2),P(PP,3))>0); % 面左边的点
P_right = PP(QPM(a,P(PP,1),P(PP,2),P(PP,3))<0); % 面右边的点
scatter3(P(P_left,1),P(P_left,2),P(P_left,3),10,'b','filled')
scatter3(P(P_right,1),P(P_right,2),P(P_right,3),10,'g','filled')
left = false(size(P_left));%用作标记左边点是否已经判断完了
right = false(size(P_right));%用作标记右边点是否已经判断完了
Pairs = []; %寻找点对
for i = 1:numel(P_left)
    epleft = P(P_left(i),:);
    disij = sqrt(sum((P(P_right,:) - epleft).^2,2));
    [mindis1,j] = min(disij);
    if(right(j)) %假设j已经被找过了，就不用找了
        left(i) = true;
        continue;
    end
    epright = P(P_right(j),:);
    disji = sqrt(sum((P(P_left,:) - epright).^2,2));
    [mindis2,i1] =  min(disji);
    if(i==i1) % 假如i和i1是同一个，那么i和j互为最近
        Pairs = [Pairs; i, j];
    end
    left(i) = true;
    right(j) = true;
end
figure(2); clf
fimplicit3(@(x,y,z)QPM(a,x,y,z),'facealpha',1,'edgecolor','k','facecolor','none') % 画切面
hold on;
P_Pairs = [P_left(Pairs(:,1)), P_right(Pairs(:,2))]; % 点的近距离配对
scatter3(P(P_Pairs(:,1),1),P(P_Pairs(:,1),2),P(P_Pairs(:,1),3),10,'b','filled')
scatter3(P(P_Pairs(:,2),1),P(P_Pairs(:,2),2),P(P_Pairs(:,2),3),10,'g','filled')
title('切平面两侧互为最近的点演示')

% 遍历两侧的互为最近点求交点集合M
M = zeros(size(P_Pairs,1),3);
for i = 1:1:size(P_Pairs,1)
    j1 = P_Pairs(i,1);% 第i对配对点的第1个点
    j2 = P_Pairs(i,2);% 第i对配对点的第2个点
    Pj1 = P(j1,:);% 第一个点的坐标
    Pj2 = P(j2,:);% 第二个点的坐标
    dP = Pj2-Pj1; % 两个点的距离
    t = -QPM(a,Pj1(1),Pj1(2),Pj1(3))/sum(a(:).*dP(:));% 直线参数方程与切平面的交点的参数t
    arr = Pj1+t*dP;% 获取交点坐标
    M(i,:) = arr;%存入集合M（M每一行代表一个交点坐标）
end
scatter3(M(:,1),M(:,2),M(:,3),10,'r','filled')
legend('切平面','左侧蓝色点','右侧红色点','连线与切平面的交点')
end