插值在面上分段常量的二維數據

我有一個不規則的網格，它由兩個變量描述 - 一個存儲構成每個面的頂點索引的faces數組，以及一個存儲每個頂點座標的vert數組。我還有一個函數，假定每個面都是分段常量，並以每個面的值數組形式存儲。插值在面上分段常量的二維數據

我正在尋找一種方法來從這些數據構造函數f。大意如下的內容：

faces = [[0,1,2], [1,2,3], [2,3,4] ...] 
verts = [[0,0], [0,1], [1,0], [1,1],....] 
vals = [0.0, 1.0, 0.5, 3.0,....] 

f = interpolate(faces, verts, vals) 

f(0.2, 0.2) = 0.0 # point inside face [0,1,2] 
f(0.6, 0.6) = 1.0 # point inside face [1,2,3]

評估f(x,y)是要找到相應的臉上那點x,y在於，返回存儲在該面價值的手工方式。有沒有一個函數已經在scipy（或matlab）中實現了這個功能？

來源

2010-02-19 dzhelil

有MATLAB中沒有內置功能，將你想要做什麼。您可以使用函數INPOLYGON作爲suggested by Jonas來構建自己的算法，但您可以使用某些標準算法自行創建更快的實現，以查找點是否在多邊形內。

前陣子我寫我自己的代碼爲一個線段並在3-d的一組三角形表面之間找出的交叉點，我發現this softsurfer link是用於實現該算法的最有幫助。我的情況比你的情況更復雜。由於您在2D中工作，因此您可以忽略關於找到線段與三角形平面相交的點的鏈接的第一部分。

我已經在下面包含了我的MATLAB代碼的簡化版本供您使用。函數interpolate將把您的faces,vertices和values矩陣作爲輸入，並返回一個函數句柄f，可以在給定的（x，y）點評估該函數句柄以獲取邊界三角形內的分段值。下面這段代碼中的幾個特點：

當你評估f將要完成的處理包含在nested functionevaluate_function。此功能可以訪問interpolate中的其他變量，因此需要預先計算三角形內測試所需的一些變量，以便evaluate_function儘可能快地運行。
如果你有很多三角形，測試你的觀點是否在所有三角形內可能是昂貴的。因此，代碼首先找到您的點的給定半徑（即三角形最長腿的長度）內的三角形。只有這些附近的三角形進行測試，看看它是否在其中。
如果一個點不落在任何三角形區域內，f將返回值NaN。

有一些事情沒有被包括在代碼中，你可能希望在這取決於你使用的是它什麼補充：

輸入檢查：代碼目前假定faces是一個N×3的矩陣，vertices是一個M乘2的矩陣，而values是一個長度爲N的向量。您可能希望添加錯誤檢查以確保輸入符合這些要求，並在出現錯誤時指出其中一個或多個錯誤。
退化三角形檢查：這有可能是一個或多個由您faces和vertices輸入所定義的三角形的可以是簡併的（即，它們可具有0的區域）。這發生在兩個或更多三角形頂點相同的精確點時，或三角形的所有三個頂點位於一條直線上時。當涉及到評估f時，您可能希望添加一個將忽略這些三角形的支票。
處理邊緣情況：點有可能最終位於兩個或更多三角形區域的邊緣。因此，您必須決定該點的價值（即最大的面值，面值的平均值等）。對於像這樣的邊緣情況，下面的代碼將自動選取更接近faces變量中面孔列表開頭的面的值。

最後，這裏是代碼：

function f = interpolate(faces,vertices,values) 

    %# Precompute some data (helps increase speed): 

    triVertex = vertices(faces(:,2),:);    %# Triangles main vertices 
    triLegLeft = vertices(faces(:,1),:)-triVertex; %# Triangles left legs 
    triLegRight = vertices(faces(:,3),:)-triVertex; %# Triangles right legs 
    C1 = sum(triLegLeft.*triLegRight,2); %# Dot product of legs 
    C2 = sum(triLegLeft.^2,2);   %# Squared length of left leg 
    C3 = sum(triLegRight.^2,2);   %# Squared length of right leg 
    triBoundary = max(C2,C3);    %# Squared radius of triangle boundary 
    scale = C1.^2-C2.*C3; 
    C1 = C1./scale; 
    C2 = C2./scale; 
    C3 = C3./scale; 

    %# Return a function handle to the nested function: 

    f = @evaluate_function; 

    %# The nested evaluation function: 

    function val = evaluate_function(x,y) 

    w = [x-triVertex(:,1) y-triVertex(:,2)]; 
    nearIndex = find(sum(w.^2,2) <= triBoundary); %# Find nearby triangles 
    if isempty(nearIndex) 
     val = NaN;   %# Return NaN if no triangles are nearby 
     return 
    end 
    w = w(nearIndex,:); 
    wdotu = sum(w.*triLegLeft(nearIndex,:),2); 
    wdotv = sum(w.*triLegRight(nearIndex,:),2); 
    C = C1(nearIndex); 
    s = C.*wdotv-C3(nearIndex).*wdotu; 
    t = C.*wdotu-C2(nearIndex).*wdotv; 
    inIndex = find((s >= 0) & (t >= 0) & (s+t <= 1),1); 
    if isempty(inIndex) 
     val = NaN;   %# Return NaN if point is outside all triangles 
    else 
     val = values(nearIndex(inIndex)); 
    end 

    end 

end

來源

2010-03-08 16:17:50 gnovice

對於我來說，這聽起來並不像我想象的那樣內插，因爲要找到哪個三角形面對的是內部的點。檢查this site瞭解測試每個三角形面的方法。你的函數只會確定它在哪裏，並返回相應的值。當然，如果你有很多面孔，或者如果你這麼做很多，那麼你會想找到方法來優化它（存儲每個三角形在x和y方向最遠的+和 - 點，並存儲例如，如果點不在這個邊界框內，那麼你也可以不檢查它是否在三角形內部）。

我真的懷疑你會發現一些內建於Matlab或scipy的東西來做你想做的事，但我可能是錯的。

來源

2010-02-19 23:22:19

您可能要使用橋接器CGAL-python模塊;如果我沒有記錯，CGAL提供三角形搜索的功能。但是，如果使用內置表示法遞增構造三角測量，它可能工作效率最高。對於快速而骯髒的情況，您可以通過Voronoi圖（在Matlab中的功能不是很好）查找離查詢點最近的網格頂點，或者對於單個查詢，計算所有距離並找到最小值，然後搜索具有該頂點的所有三角形。

來源

2010-02-20 02:01:02 shabbychef

Matlab內置函數inpolygon，它允許你測試你是否在一個三角形內。我不知道哪種功能可以識別你面對的是什麼。

如果您要編寫這樣一個函數，我會首先測試您的點最接近哪個頂點，然後評估共享該頂點的所有面的inpolygon，直到找到匹配。這應該是相當快的。

來源

2010-02-20 02:01:14 Jonas

看一看matplotlib.delaunay.interpolate，這是一個很好的C代碼包裝。
（但是class LinearInterpolator有說「目前，只有整齊的矩形網格進行插值支持。」）

來源

2010-03-07 20:08:49 denis

插值在面上分段常量的二維數據

回答

相關問題