斷層面擬合（Matlab實現）

斷層面擬合(Matlab)

一、算法目的

斷層面的可視化對於地質構造分析、油氣勘探開發具有重要作用。常用的斷層面重構方法是利用斷層數據進行曲面擬合，然而的斷層數據一般用離散不規則分布的數據散點雲來反映，需要基於這些散點雲數據進行三維地質建模。建模過程需要利用插值算法與計算機圖形學以獲得光滑、連續的斷層面。

算法將實現散點雲\(\Longrightarrow\)斷層面的過程。

二、算法流程

三、程序結構

flowchart LR; Z(main.m)-->A; A(normfaultdata.m)-->B(planeSurfaceFitting.m); B-->C(projectplane.m); C-->D(Convert3Dto2D.m); D-->E(inversenormfaultdata.m);

四、程序分析

0 讀取斷層數據

待處理的數據格式如下表所示，為文本格式，程序首先從origin_point.txt中獲取斷層散點

大地坐標X	大地坐標Y	時間T	線號	道號	所屬斷層線索引
598449.6	5132452	1473.383	786	1682	0
598470	5132452	1535.572	786	1685	0
598489.1	5132452	1560.448	786	1686	0
598619.9	5132532	1455.97	790	1699	1
598655.4	5132532	1516.916	790	1703	1
598698.9	5132532	1569.154	790	1707	1
598719.4	5132612	1454.726	794	1709	2
598749.4	5132612	1526.866	794	1712	2

faultorg=load("origin_point.txt");
X=faultorg(:,1);
Y=faultorg(:,2);
Z=-faultorg(:,3);
LineIndex=faultorg(:,4);
CdpIndex=faultorg(:,5);
Index=faultorg(:,6);

1 歸一化、中心化

歸一化是為了平衡各維度數據的量綱，減小方差大的特征的影響，並加快算法計算速度；

本次采用的歸一化方法為線性歸一化(Min-Max scaling)，將原始數據轉到\([0,1]\)的范圍內，公式為：

\[{X_{norm}} = \frac{{X - {X_{\min }}}}{{{X_{\max }} - {X_{\min }}}} \]

中心化的操作則是對各列數據減去列平均值，目的是將坐標原點放在數據的中心。

[Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z);

normfaultdata.m

function [Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z)

%{
*函數功能：* ***完成斷層數據中心化、歸一化***

*輸入：*

- X：斷層面的X坐標
- Y：斷層面的Y坐標
- Z：斷層面的Z坐標

*輸出：*

- Xnorm/Ynorm/Znorm：空間坐標的歸一化值
- Xmean/Ymean/Zmean：空間坐標的平均值
- Xmax/Ymax/Zmax：空間坐標的最大值
- Xmin/Ymin/Zmin：空間坐標的最小值
%}

Xmax=max(X); 
Xmin=min(X);
Ymax=max(Y);
Ymin=min(Y);
Zmax=max(Z);
Zmin=min(Z);

Xnorm1=(X-Xmin)/(Xmax-Xmin);   %歸一化
Ynorm1=(Y-Ymin)/(Ymax-Ymin);
Znorm1=(Z-Zmin)/(Zmax-Zmin);


Xmean=mean(Xnorm1);  % 數據中心化
Ymean=mean(Ynorm1);
Zmean=mean(Znorm1);

Xnorm=Xnorm1-Xmean;
Ynorm=Ynorm1-Ymean;
Znorm=Znorm1-Zmean;

end

2 擬合投影平面

投影的目的是將三維空間降低一個維度，從而簡化空間中的插值問題，投影的實現將采用主成分分析(PCA)算法。

首先簡要介紹PCA算法。

引用自https://www.bilibili.com/video/BV1E5411E71z?from=search&seid=17363694020924718907&spm_id_from=333.337.0.0

考慮這樣一種情況

相當於把二維的信息降到一維，所以PCA的目的就是找這樣一個坐標系，使得數據降維后信息損失最小，即信息保留最多。

假設我們手上的數據已服從正態分布，在經過中心化、歸一化之后，那它可以由白數據經過拉伸、旋轉變換得到。

可以看出求解\(\boldsymbol R\)矩陣是獲取新坐標系的關鍵。

概括PCA的求解過程：

flowchart LR; A(三維空間散點)-->B(歸一化<br/>中心化); B-->C(找投影坐標系<br/>使得數據盡可能分散) C-->D(找方差最大方向) D-->E(求旋轉矩陣R<br/>) E-->F(求協方差矩陣的<br/>特征向量) F-->G(投影面方程) G-->H(二維投影點)

(1)協方差矩陣

\[\operatorname{cov}(x, y)=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{n-1}\\ \operatorname{cov}(x, x)=\frac{\sum_{i=1}^{n}x_{i}^2}{n-1}(即方差)\\ \]

協方差表示的是兩個向量在變化過程中是同方向變化還是反方向變化，同向或反向的程度如何，協方差為正時，\(x\)，\(y\)為正相關。而協方差矩陣表示為：

\[\begin{align} C = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm cov}} (x,x)}&{{\mathop{\rm cov}} (x,y)}\\ {{\mathop{\rm cov}} (x,y)}&{{\mathop{\rm cov}} (y,y)} \end{array}} \right]& = \left[ {\begin{array}{*{20}{c}} {\frac{{\sum\limits_{i = 1}^n {x_i^2} }}{{n - 1}}}&{\frac{{\sum\limits_{i = 1}^n {{x_i}} {y_i}}}{{n - 1}}}\\ {\frac{{\sum\limits_{i = 1}^n {{x_i}} {y_i}}}{{n - 1}}}&{\frac{{\sum\limits_{i = 1}^n {y_i^2} }}{{n - 1}}} \end{array}} \right] \\ & =\frac{1}{{n - 1}}\left[ {\begin{array}{*{20}{l}} {{x_1}}&{{x_2}}&{{x_3}}&{{x_4}}\\ {{y_1}}&{{y_2}}&{{y_3}}&{{y_4}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{x_1}}&{{y_1}}\\ {{x_2}}&{{y_2}}\\ {{x_3}}&{{y_3}}\\ {{x_4}}&{{y_4}} \end{array}} \right]\\ &= \frac{1}{{n - 1}} D {D^T} \end{align} \]

白數據中\(x\)，\(y\)不相關，其協方差矩陣為\(C_0=\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]\)，於是我們手上的數據的協方差：

\[\begin{align} C & = \frac{1}{{n - 1}}D{D^T} = \frac{1}{{n - 1}}RSD_0{(RSD_0)}^T\\ & = \frac{1}{{n - 1}}RSD_0D_0^TS^TR^T = RS(\frac{1}{{n - 1}}D_0D_0^T)S^TR^T\\ & = RSC_0S^TR^T = RS\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]S^TR^T=RSS^TR^T \\ \end{align} \]

又有拉伸變換\(S^T=S\)，旋轉變換\(R^T=R^{-1}\)，令\(L=SS^T=\left[ {\begin{array}{*{20}{c}} a^2&0\\ 0&b^2 \end{array}} \right]\)，於是原式化為

\[C=RLR^{-1}\Rightarrow CR=RL \]

(2)協方差矩陣的特征值

\(L=\left[ {\begin{array}{*{20}{c}} \lambda_1&0\\ 0&\lambda_2 \end{array}} \right]=SS^T=\left[ {\begin{array}{*{20}{c}} a^2&0\\ 0&b^2 \end{array}} \right]\)，\(a^2\)和\(b^2\)是兩個軸方向的方差（拉伸倍數的平方），同時又是協方差矩陣的特征值。

(3)協方差矩陣的特征向量

\[C{{\vec v}} = \lambda {{\vec v}} \]

\(R\)矩陣由\({\vec v}_1\)、\({\vec v}_2\)組成，\({\vec v}_1\)為新坐標系的\(X\)軸方向，\({\vec v}_2\)為新坐標系的\(Y\)軸方向，或者說這兩個向量就是PCA的主成分。

(4)PCA與空間變換

從上文的介紹可以看出，PCA並未對數據做任何改變，它僅切換了坐標軸。其主成分（特征向量）則恰好表示新的基向量，例如從一開始的\(\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right]\)變為\(\left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right]\)，則表示X軸由通常理解的方向變為向量\([-0.63,0.35,-0.68]\)指向的方向。但這一過程僅變動了坐標軸指向，並未發生線性變換，即並未帶動數據一起轉動。

當僅討論三維數據時，PCA的第一、第二主成分即為方差最大（\(X\)軸方向）、次大（\(Y\)軸方向）的方向，而最小特征值對應的特征向量則是\(Z\)軸，這個新的\(Z\)軸是降維平面\(X-Y\)的法向量。設其為\(v_1=[a,b,c]\)，中心點(Xmean,Ymean,Zmean)過平面，易得截距\(d\)，於是PCA投影所成的平面可表示為：

\[ax+by+cz+d=0 \]

本節中將用到PCA擬合投影面和獲取投影坐標點兩個函數：

[a,b,c,d]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean);

planeSurfaceFitting.m

function [a,b,c,d,V2]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean)

%{
*函數功能：* ***斷層面數據擬合斷層平面***

*輸入：*

- Xnorm/Ynorm/Znorm：空間坐標的歸一化值
- Xmean/Ymean/Zmean：空間坐標的平均值

*輸出：*

- a/b/c/d：斷層面多項式系數
- V2：得到的特征向量
%}

C=[Xnorm,Ynorm,Znorm];		% 輸入的矩陣為歸一化后的空間坐標
Cov=C'*C;					% CT*C得到方陣
[eigenV,eigenU]=eig(Cov);	% 求出特征值
[~,IndexValue]=min(diag(eigenU)); %最小的特征值對應的向量即為法向量
[~,Index]=sort(diag(eigenU),'descend');

V2=eigenV(:,Index);		%V2中存放了所有的特征向量
              
V1=eigenV(:,IndexValue);
a=V1(1);
b=V1(2);
c=V1(3);
d=-dot([a b c],[Xmean,Ymean,Zmean]); %均值點過投影平面，算出截距
 
end

接下來根據原坐標點和投影平面，獲取投影點坐標。

[x_projection,y_projection,z_projection]=projectplane(Xnorm,Ynorm,Znorm,[a,b,c,d]);

projectplane.m

function [x_projection,y_projection,z_projection]=projectplane(X,Y,Z,cfitt)

%{
*函數功能 * ***獲取投影點坐標***

*輸入：*

- X/Y/Z：原始空間坐標
- cfitt：斷層面多項式系數列表

*輸出：*

- x_projection/y_projection/z_projection：擬合面上投影點的空間坐標
%}

    x_projection=zeros(1,length(X));
    y_projection=zeros(1,length(Y));
    z_projection=zeros(1,length(Z));
    A0=cfitt(1);
    B0=cfitt(2);
    C0=cfitt(3);
    D0=cfitt(4);

   for i=1:length(X)
       temp=-(A0*X(i)+B0*Y(i)+C0*Z(i)+D0)/(C0^2+A0^2+B0^2);		%點到平面的距離公式
       
      x_projection(i)=X(i)+A0*temp;		%投影之后的新的x,y,z
      y_projection(i)=Y(i)+B0*temp;
      z_projection(i)=Z(i)+C0*temp;
   end
end

隨后匯總和調整新坐標系下的數據

data=[x_projection;y_projection;z_projection];	%新坐標存入data中
ori = (data(:,1)+data(:,end))/2;  %取始末點的坐標中點為新的坐標原點
xorder=data(:,1)-data(:,end);%選取原始始末點相減為投影后x軸方向
xor=xorder/norm(xorder);%2維坐標下的x軸方向向量

3 旋轉投影平面

矩陣與線性變換

回到線性變換與矩陣的本質，一個二維線性變換可以僅由四個數字完全確定，變換后的\(i\)的兩個坐標和變換后的\(j\)的兩個坐標，將這兩個坐標封裝在2×2的矩陣中。現在有一個描述線性變換的2×2矩陣\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]\)，以及任意初始向量\(\vec v=(x,y)\)，我們將很容易得到變換之后的新向量，並通過矩陣乘法描述

\[x\left[ {\begin{array}{*{20}{c}} a\\ c \end{array}} \right] + y\left[ {\begin{array}{*{20}{c}} b\\ d \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y \end{array}} \right] \]

這時的\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]\)就可以看作一種函數，其逆矩陣\(\left[ {\begin{array}{*{20}{c}} a&b\\ c&d \end{array}} \right]^{-1}\)則相當於“重置”操作，這一思想將在后續程序中加以運用。

參考線性代數的本質 - 03 - 矩陣與線性變換

基向量的改變引起空間變換，而矩陣記錄了變換之后基向量的落點。在上一步PCA中得到的基向量為\(\left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right]\)，而接下來要執行的旋轉投影平面，則需要將這一基向量“校正”至“初始狀態”，即使用拉伸(\(S\))、旋轉(\(R\))變換將基向量\(\vec i\)重新指向\([1,0,0]\)，基向量\(\vec j\)重新指向\([0,1,0]\)，基向量\(\vec k\)重新指向\([0,0,1]\)。反投影則只需對變換矩陣\(RS\)求逆，反變換即可。

\[RS\left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right]\\ \left[ {\begin{array}{*{20}{c}} { - 0.63}&{0.43}&{0.64}\\ {0.35}&{0.90}&{ - 0.25}\\ { - 0.68}&{0.07}&{ - 0.72} \end{array}} \right] = (RS)^{-1}\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ 0&0&1 \end{array}} \right] \]

將三維空間中的投影平面通過旋轉、平移至水平

[data_xyzor,R1]=Convert3Dto2D([a,b,c,d],xor,ori,norm_data);

function [data_xyzor,R1]=Convert3Dto2D(cfitt,xor,ori,data)

%{
*函數功能：*將三維空間中的X-Y平面旋轉、平移至水平面

*輸入：*

- cfitt：斷層面多項式系數列表
- xor：投影面 x軸方向向量
- ori：投影面原點坐標
- data：待旋轉、投影的數據/散點


*輸出：*

- data_xyzor：旋轉、平移后的散點
- R1：變換矩陣（旋轉*平移）
%}


aall=cfitt(1);ball=cfitt(2); call=cfitt(3); 	%從輸入列表中提取出a,b,c
zor=[aall,ball,call]/norm([aall,ball,call]);	%匯總到法向量z，並標准化z方向
yor=cross(xor,zor);yor=yor/norm(yor);			%現在有了x方向和z方向的方向向量，用叉乘算出y方向向量

move=[1,0,0,0;0,1,0,0;0,0,1,0;-ori',1];			%平移矩陣，以新的坐標原點去描述平移的尺度
%旋轉矩陣，以xor為例，其為新的空間中的一個坐標軸，上面已經標准化為單位長度，即xor只剩下旋轉的分量，表示從最初的x軸--->xor所需要的變換向量
rotate=[xor(1),yor(1),zor(1),0;xor(2),yor(2),zor(2),0;xor(3),yor(3),zor(3),0;0,0,0,1];

R1=move*rotate;  %R1為整個變換矩陣

data(4,:)=1;	% 將數據補充至四維，增加一個平移維度
data_xyzor=(data')*R1;%  乘變換矩陣，所有散點跟隨坐標系旋轉到X、Y軸水平的狀態

end

4 限制插值區域

這一步主要是為了裁剪掉后續多余的插值點，使用Matlab中的凸包convhull工具，根據現有的基礎散點，作出其凸包絡線，對插值范圍進行約束。

Xr=data_xyzor(:,1);		%獲得數據點新的x,y,z坐標
Yr=data_xyzor(:,2);
Zr=data_xyzor(:,3);
[k,a1]=convhull(Xr,Yr);
xv=Xr(k);
yv=Yr(k);

此時已經獲得二維平面的散點，\(T\)值將是待插屬性。

經此步驟，得到斷層散點的邊界線，如圖；

5 二維平面插值

現在手上的數據變成了(Xr，Yr)以及對應坐標上的屬性值Tr，目標是在圈定的范圍內插出新的、過渡的屬性值。這一環節已簡化為二維的插值任務。

本次插值采用的是雙調和樣條插值(biharmonic spline interpolation)算法（matlab中的V4算法），是不規則間隔的二維數據點的插值，簡要介紹該算法：

插值平面是以每個數據點為中心的格林函數的線性組合。 通過求解線性方程組可以找到格林函數的幅度。

插值平面\(S(\boldsymbol{x})\)的表達式為：

\[S(\boldsymbol{x})=\sum_{j=1}^{n} w_{j} g(\boldsymbol{x},\boldsymbol{x}_{j})\\ g\left(\boldsymbol{x}_{i}, \boldsymbol{x}_{j}\right)=r_{ij}^{2} \ln \left(r_{ij}\right)\quad i=1,2, \ldots, n,\quad j=1,2, \ldots, n\\ r_{ij}=\left|\boldsymbol{x}_{i}, \boldsymbol{x}_{j}\right|=\sqrt{\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}}, \quad i=1,2, \ldots, n,\quad j=1,2, \ldots, n \]

其中\(\boldsymbol x\)是待插值的樣本點坐標，\(g(\boldsymbol{x},\boldsymbol{x}_{j})\)為格林函數，\(\boldsymbol{x}_{i}\)、\(\boldsymbol{x}_{j}\)分別是第個\(i、j\)個已知數據，\(n\)是平面的總已知點數，\(\boldsymbol{w}_{j}\)是\(\boldsymbol{x}與\boldsymbol{x}_{j}\)的相關性的權重矩陣。

算法步驟如下：

①假設插值平面所需的輸出節點是\(p_0\)，最接近\(p_0\)的插值點有\(k\)個，\(p_0\)到第\(j\)個已知點的距離按下式計算：

\[r_{0j}=\left|\boldsymbol{x}_{0}, \boldsymbol{x}_{j}\right|, \quad j=1,2, \ldots, n \]

其中\(\boldsymbol{x}_{0}\)是\(p_0\)的位置向量；所有距離都使用合並排序算法按升序排序，在排序之后，獲取前\(k\)個數據點的索引和坐標；

②計算這\(k\)個數據點的格林函數矩陣\(\boldsymbol G\)：

\[\mathbf{G}=\left[\begin{array}{cccc} d_{11} & d_{12} & \cdots & d_{1 k} \\ d_{21} & d_{22} & \cdots & d_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ d_{k 1} & d_{k 2} & \cdots & d_{k k} \end{array}\right]\\ d_{i j}=\left\{\begin{array}{ll} 0, & \text { when } r_{i j}=0, \\ {\left[\ln \left(r_{i j}\right)-1.0\right] r_{i j}^{2},} & \text { when } r_{i j} \neq 0, \end{array} {i=1, \ldots, k, j=1, \ldots, k}\right. \]

根據格林函數矩陣和屬性矩陣，權重矩陣可以計算：

\[\mathbf{W} = \mathbf{G}^{-1}\mathbf{Z} \]

③計算點\(p_0\)的\(1×k\)的格林函數矩陣$ \mathbf{G}_{p_0}$：

\[\mathbf{G}_{p_0}=\left[\begin{array}{llll} d_{01} & d_{02} & \ldots & d_{0 k} \end{array}\right] \]

易得，點\(p_0\)的屬性值\(z_p\)：

\[z_{p_0}=\mathbf{G}_{p_0}\mathbf{W} \]

④重復步驟"①~③"，計算其他點\(p_i\)插值后的屬性值\(z_{p_i}\).

程序將使用Matlab中的griddata函數中的v4算法，在平面上對\(T\)值進行插值。

Xrmax=max(Xr);		%最大最小值坐標
Xrmin=min(Xr);
Yrmax=max(Yr);
Yrmin=min(Yr);
h=0.01;     %計算步長
[Xin,Yin]=meshgrid(Xrmin:h:Xrmax,Yrmin:h:Yrmax);  %搭建插值網格，假設X方向上的采樣點為100個，Y方向上50個，那么Xin和Yin都將是100×50的矩陣，Xin與Yin組合，方能表示矩陣中任一位置的坐標
Zin=griddata(Xr,Yr,Zr,Xin,Yin,'v4'); %進行插值，Zin的大小將和Xin、Yin相同，Zin也可延申為在這個平面網格上的其他屬性值

插值完成后，剔除范圍外的插值點：

% 先將Xin,Yin,Zin展成一維，因為要配合裁剪的xv與yv也是一維的
[ny,nx]=size(Zin);
Xin_row1=zeros(1,ny*nx);
Yin_row1=zeros(1,ny*nx);
Zin_row1=zeros(1,ny*nx);

for iy=1:ny
    for ix=1:nx
        index=(iy-1)*nx+ix;
        Xin_row1(index)=Xin(iy,ix);
        Yin_row1(index)=Yin(iy,ix);
        Zin_row1(index)=Zin(iy,ix);
    end
end

% 同樣是Matlab工具箱中的函數
index_cut=inpolygon(Xin_row1,Yin_row1,xv,yv);

% 得到插值、裁剪后的一維數據
Xin_row=Xin_row1(index_cut);
Yin_row=Yin_row1(index_cut);
Zin_row=Zin_row1(index_cut);

插值后的效果為：

6 散點雲三角剖分

插值后的數據間隔已經非常小，但依然是空間中的三維散點，需要將這些散點連接成面，才算達到構建斷層面的效果。此時將用到點集的三角剖分技術(Triangulation)，三角剖分研究的是如何把散點集合剖分成不均勻的三角形網格。其中Delaunay三角剖分算法較為成熟，廣泛應用與各領域圖像和點雲的處理上。

使用matlab中的庫函數delaunay實現

%進行Delauny三角剖分算法 
%輸入為點坐標P，以向量/矩陣形式指定
tri=delaunay(Xin_row,Yin_row); % 輸出tri的每行都包含P的行號（索引），它們構成三角剖分中的單個三角形

三角剖分后的平面如圖，此時只需按照Z值進行渲染，便能夠得到最終的插值斷層面。

7 反投影

pr_in=[Xin_row',Yin_row',Zin_row',ones(length(Zin_row),1)];  %同樣增加一維
Faultnorm_in=pr_in*inv(R1);		%乘變換矩陣的逆，返回原坐標系

8 反歸一化

Xin_norm=Faultnorm_in(:,1);
Yin_norm=Faultnorm_in(:,2);
Zin_norm=Faultnorm_in(:,3);
[Xorg,Yorg,Zorg]=inversenormfaultdata(Xin_norm,Yin_norm,Zin_norm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin);

function [Xorg,Yorg,Zorg]=inversenormfaultdata(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin)

%{
*函數功能：*斷層數據反歸一化

*輸入：*

- Xnorm/Ynorm/Znorm：斷層面上的歸一化空間坐標
- Xmean/Ymean/Zmean：斷層面上空間坐標的平均值
- Xmax/Ymax/Zmax：斷層面上空間坐標的最大值
- Xmin/Ymin/Zmin：斷層面上空間坐標的最小值

*輸出：*

- Xorg/Yorg/Zorg：反歸一化后的空間坐標，即原始坐標系下的坐標值
%}

    Xorg1=Xnorm+Xmean;
    Yorg1=Ynorm+Ymean;
    Zorg1=Znorm+Zmean;

    Xorg=Xorg1*(Xmax-Xmin)+Xmin;
    Yorg=Yorg1*(Ymax-Ymin)+Ymin;
    Zorg=Zorg1*(Zmax-Zmin)+Zmin;  
    
end

得到插值、擬合后的最終斷層面

main函數完整代碼

clc
clear all
close all


%% 導入斷層散點 
faultorg=load("origin_point.txt");
X=faultorg(:,1);
Y=faultorg(:,2);
Z=-faultorg(:,3);
LineIndex=faultorg(:,4);
CdpIndex=faultorg(:,5);
Index=faultorg(:,6);

%% 數據歸一化、中心化
[Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin]=normfaultdata(X,Y,Z);

%% Step three: 斷層面PCA算法擬合及投影
[a,b,c,d]=planeSurfaceFitting(Xnorm,Ynorm,Znorm,Xmean,Ymean,Zmean);
%坐標旋轉投影 
[x_projection,y_projection,z_projection]=projectplane(Xnorm,Ynorm,Znorm,[a,b,c,d]);

data=[x_projection;y_projection;z_projection];


ori=(data(:,1)+data(:,end))/2;%投影后的坐標原點，即轉換為2維坐標后的原點 
xorder=data(:,1)-data(:,end);%選取原始始末點相減為投影后x軸方向
xor=xorder/norm(xorder);%2維坐標下的x軸方向向量

norm_data = [Xnorm';Ynorm';Znorm']
[data_xyzor,R1]=Convert3Dto2D([a,b,c,d],xor,ori,norm_data);

Xr=data_xyzor(:,1);
Yr=data_xyzor(:,2);
Zr=data_xyzor(:,3);

%% 形成凸多邊形
[k,a1]=convhull(Xr,Yr);
xv=Xr(k);
yv=Yr(k);

%%  插值
Xrmax=max(Xr);
Xrmin=min(Xr);
Yrmax=max(Yr);
Yrmin=min(Yr);
h=0.01;     %計算步長
[Xin,Yin]=meshgrid(Xrmin:h:Xrmax,Yrmin:h:Yrmax);  %搭建插值網格
Zin=griddata(Xr,Yr,Zr,Xin,Yin,'v4'); %進行插值 

[ny,nx]=size(Zin);
Xin_row1=zeros(1,ny*nx);
Yin_row1=zeros(1,ny*nx);
Zin_row1=zeros(1,ny*nx);

for iy=1:ny
    for ix=1:nx
        index=(iy-1)*nx+ix;
        Xin_row1(index)=Xin(iy,ix);
        Yin_row1(index)=Yin(iy,ix);
        Zin_row1(index)=Zin(iy,ix);
    end
end

%剔除部分點
index_cut=inpolygon(Xin_row1,Yin_row1,xv,yv);

Xin_row=Xin_row1(index_cut);
Yin_row=Yin_row1(index_cut);
Zin_row=Zin_row1(index_cut);

%進行Delauny三角剖分算法 
tri=delaunay(Xin_row,Yin_row); 

%反向旋轉投影變化
pr_in=[Xin_row',Yin_row',Zin_row',ones(length(Zin_row),1)];
Faultnorm_in=pr_in*inv(R1);

%反歸一化操作恢復原始坐標
Xin_norm=Faultnorm_in(:,1);
Yin_norm=Faultnorm_in(:,2);
Zin_norm=Faultnorm_in(:,3);
[Xorg,Yorg,Zorg]=inversenormfaultdata(Xin_norm,Yin_norm,Zin_norm,Xmean,Ymean,Zmean,Xmax,Ymax,Zmax,Xmin,Ymin,Zmin);