---

前言

上一篇文章感覺對三角剖分問題沒有說清楚，這次專門對三角剖分問題再仔細說說。

---

一、問題引入

實際上這個問題是用來解決二維曲面插值問題的。

二維插值問題，用matlab的一些函數就可以方便操作，比如 interp2 。但 interp2函數是用在規則點數據集的情況下，比如已知“密度較稀疏”的一些數據點，進行插值，找到“密度適中”的點。下面舉個例子說明。

% 生成自定義的網格點
x = linspace(-3, 3, 10);
y = linspace(-3, 3, 15);
[X, Y] = meshgrid(x, y);% 計算相應的高度值
Z = peaks(X, Y);% 繪制原始網格圖
figure;
subplot(1, 2, 1);
surf(X, Y, Z);
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Original Peaks Surface');
axis tight;
grid on;% 指定插值后的網格點
xi = linspace(-3, 3, 30);
yi = linspace(-3, 3, 45);
[XI, YI] = meshgrid(xi, yi);% 使用插值方法計算插值后的高度值
ZI = interp2(X, Y, Z, XI, YI, 'spline');% 繪制插值后的網格圖
subplot(1, 2, 2);
surf(XI, YI, ZI);
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Interpolated Peaks Surface');
axis tight;
grid on;

得到圖如下：

但對于一些無序（或者說在空間排列沒有規律）的散點進行插值，這時要使用?三角剖分插值

具體概念介紹可以參考下面的wiki鏈接

https://en.wikipedia.org/wiki/Delauna4_triangulation

二、一個例子

1.生成散點圖

為了好說明，我們在1/4半球面上進行操作，隨機選取球面上的一些點作為散點，同時畫出它的xy面投影剖分（后面要用）

% 定義半徑和繪制分辨率
radius = 1;  % 半徑
resolution = 50;  % 分辨率% 生成球面上的坐標點
theta = linspace(-pi/2, 0, resolution);
phi = linspace(0, pi/2, resolution);
[THETA, PHI] = meshgrid(theta, phi);
x = radius * sin(PHI) .* cos(THETA);
y = radius * sin(PHI) .* sin(THETA);
z = radius * cos(PHI);% 隨機選擇50個點
numPoints = 50;
indices = randperm(resolution^2, numPoints);
selectedPoints = [x(indices(:)), y(indices(:)), z(indices(:))];% 進行三角剖分
tri = delaunay(selectedPoints(:, 1), selectedPoints(:, 2));% 繪制1/4半球面
figure;
surf(x, y, z);
hold on;% 繪制隨機選擇的點
scatter3(selectedPoints(:, 1), selectedPoints(:, 2), selectedPoints(:, 3), 'filled', 'r');

?

2.對數據進行剖分

代碼如下：

clear all
close all
clcrng(10)
% 定義半徑和繪制分辨率
radius = 1;  % 半徑
resolution = 50;  % 分辨率% 生成球面上的坐標點
theta = linspace(-pi/2, 0, resolution);
phi = linspace(0, pi/2, resolution);
[THETA, PHI] = meshgrid(theta, phi);
x = radius * sin(PHI) .* cos(THETA);
y = radius * sin(PHI) .* sin(THETA);
z = radius * cos(PHI);% 隨機選擇點
numPoints = 50;
indices = randperm(resolution^2, numPoints);
selectedPoints = [x(indices(:)), y(indices(:)), z(indices(:))];save selectedPoints.mat selectedPoints% 在xy平面上進行平面剖分
dt = delaunayTriangulation(selectedPoints(:, 1), selectedPoints(:, 2));
tri = dt.ConnectivityList;% 根據剖分點的坐標和對應的z值生成三維空間中的三角網格
tri3D = [tri, tri(:, 1) + size(selectedPoints, 1), tri(:, 2) + size(selectedPoints, 1)];x3D = [selectedPoints(:, 1); selectedPoints(:, 1); selectedPoints(:, 1)];
y3D = [selectedPoints(:, 2); selectedPoints(:, 2); selectedPoints(:, 2)];
z3D = [selectedPoints(:, 3); selectedPoints(:, 3); selectedPoints(:, 3)];% 繪制1/4半球面
figure;
h = surf(x, y, z);
set(h, 'FaceAlpha', 0.5);  % 設置表面的透明度
% set(h, 'FaceColor', 'green');  % 設置表面顏色為空
hold on;% 繪制隨機選擇的點
scatter3(selectedPoints(:, 1), selectedPoints(:, 2), selectedPoints(:, 3), 'filled', 'r');triplot(dt);% % 繪制三角網格
% trisurf(tri3D, x3D, y3D, z3D, 'FaceColor', 'none', 'EdgeColor', 'b', 'FaceAlpha', 0.5);% 繪制三角網格
patch('Faces', tri3D, 'Vertices', [x3D, y3D, z3D], 'FaceColor', 'red', 'EdgeColor', 'b', 'FaceAlpha', 0.5);% 設置坐標軸和標題
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Quarter Sphere with Random Points and Triangulation');% 設置坐標軸比例和網格
axis equal;
grid on;

這個顯示的有點復雜了，它是幾個數據圖像的結合 包括 原始數據（網格圖）、散點圖（紅色）、空間和平面三角剖分圖（藍色）

我們去掉原始網格圖，只留平面剖分和對應曲面的映射，看看如下

待插值的點會在這些紅色三角面上找到對應的z值，因為散點插值可不同于interp2插值，根本沒有"可以依賴"現成附近的方形網格點用來估算，需要借助剖分的找到它附近的點。好了，思路有了，流程化的東西如下：

三角剖分的流程

1、對空間散點的xy平面投影進行三角剖分（注意：并不是直接在空間曲面上進行三角剖分，而是對平面進行，因為使用delaunayTriangulation會對xyz三維數據直接給出的四面體立體剖分！即它會認為是立體剖分）

2、對待插值點的xy平面投影點，找到它屬于xy平面剖分的哪個三角形（注意，是在xy平面）

3、在空間對對應三角形建立平面方程，然后使用點法式方式對待插值點求出z的值

平面和立體對應關系如下圖（共同的xy，z當然不同了）

3.點法式分析

參考代碼，還是沿用上一次提到的線性三角剖分的matlab代碼

https://www.mathworks.com/matlabcentral/fileexchange/38925-linearly-interpolate-triangulation

?這代碼核心的部分在這里：

其中點法式大家估計印象不是很深刻了，需要復習下空間解析幾何的一點知識 ，兩頁ppt幫大家回疑

?法向量怎么求呢，相當于在平面中兩個矢量的叉乘！我們翻一下matlab cross的內容

比如兩個矢量 V1 = [x2-x1,y2-y1,z2-z1],V2=[x3-x1,y3-y1,z3-z1],，記為 (a1,a2,a3)? (b1,b2,b3)

叉乘的結果是

A=a2*b3 - a3*b2=(y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)

B = a3*b1-a1*b3=(z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)

C = a1*b2-a2*b1 = (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)?

A*(xi-x1)+B*(yi-y1)+C*(zi-z1) = 0

zi =??(-((y3 - y1) * (z2 - z1) - (z3 - z1) * (y2 - y1)) * (x - x1) - ((z3 - z1) * (x2 - x1) - (x3 - x1) * (z2 - z1)) * (y - y1)) / ((x3 - x1) * (y2 - y1) - (y3 - y1) * (x2 - x1)) + z1

z1變到分子上去，然后分子變成 z1*XXX+z2*YYY+z3*CCC ，XXX，YYY，CCC就是代碼中N1 N2 N3的分子

這樣按照待插值的網格的點的順序，依次計算即可得到全部的插值數據。

---

三、最后結果

簡單對幾個點進行插值，插值之后的空間點是黃色：

load selectedPoints.mat%散點
x = selectedPoints(:,1);
y = selectedPoints(:,2);
z = selectedPoints(:,3);% 定義插值點的網格
n_points = 5; % 插值點個數xi = linspace(min(x), max(x), n_points); % x 坐標范圍
yi = linspace(min(y), max(y), n_points); % y 坐標范圍
[XI, YI] = meshgrid(xi, yi); % 插值點的網格%x y z數據同前% 構建三角剖分
DT = delaunayTriangulation(x, y);% Get the connectivity table
tri = DT.ConnectivityList;
tri = tri(:, [1, 2, 3]);ZI=interptri(tri,x,y,z,XI,YI);%  繪制插值結果
figure(1)
hold on
% surf(XI, YI, ZI);scatter3(XI(:), YI(:), ZI(:), 'filled', 'y');xlabel('X');
ylabel('Y');
zlabel('Z');

很顯然，原始插值點密集的話，插出來的曲面會更理想。

總結：空間曲面散點的三角剖分線性插值是一種常用的方法，用于從離散的散點數據中構建連續的曲面模型。該方法基于三角剖分技術，將散點分布的空間曲面劃分為一系列三角形，然后利用線性插值來估計每個三角形內部的數據點。