1、MATLAB 程式設計進階篇 內插法,張智星 jangcs.nthu.edu.tw http:/www.cs.nthu.edu.tw/jang 清大資工系 多媒體檢索實驗室,9-1 一維內插法,一維內插（1-D Interpolation） 根據一組已知的資料點（包含輸入及輸出，其中輸入是一維資料，輸出也是一維資料） 建立一個連續的函數 算出任意輸入資料點所對應的輸出值 一維內插的方法很多，MATLAB 提供了兩種基本方法： 基於多項式的內插法 基於FFT（Fast Fourier Transform，快速傅立葉轉換）的內插法,一維內插範例-1 (I),interp1指令，其原理是利用多項式來
2、進行內插運算 使用語法為 yi = interp1(x, y, xi, method) 向量 x 是資料點的 x 座標（輸入值） 向量 y 是資料點的 y 座標（輸出值） 向量 xi 是內插點（輸入值，未知其對應輸出值） 字串 method 則指定使用的方法 nearest：鄰近點內插法 linear：線性內插法 spline：片段式的三次 Spline 內插法 pchip：保持形狀的片段式三次內插 cubic：和pchip 一樣 v5cubic：MATLAB 5 所用的三次內插法,一維內插範例-1 (II),範例9-1：interp101.m,x = 0:1:4*pi; y = sin(x)
3、.*exp(-x/5); xi = 0:0.1:4*pi; y1 = interp1(x, y, xi, nearest); y2 = interp1(x, y, xi, linear); y3 = interp1(x, y, xi, spline); y4 = interp1(x, y, xi, cubic); plot(x, y , o, xi ,y1 ,g- ,xi ,y2 ,r: ,xi ,y3 ,k-., xi, y4, b-); legend(Original, Nearest, Linear, Spline, Cubic);,一維內插範例-1 (III),由圖可看出，Spline
4、 和 Cubic 所產生的曲線較平滑，但它們所需的計算時間也較久,四種內插法比較表,提示,使用 interp1指令 向量 x 必須是嚴格遞增或遞減 元素之間不必是等間隔 xi 的範圍必須落在 x 的範圍之內 如果 xi 的範圍在 x 的範圍之外，可以使用yi = interp1(x, y, xi, method, extrap) 的方式來進行外插（Extrapolation），以求得在範圍外的 xi 元素所對應的 yi 值,一維內插範例-1 (IV),可以讓使用者拖放取樣點，內插曲線就會跟著改變 範例 interpAnim01.m 使用者可以在上述圖形視窗中，拖放每一個取樣點，就可以立刻看到內
5、插曲線的變化，也同時可以知道各種內插方法的特性,一維內插範例-2 (I),interpft指令，可進行基於 FFT（Fast Fourier Transform，快速傅立葉轉換）的內插法 先計算資料點的傅立葉轉換，再用更密集的內插點來進行反傅立葉轉換 使用語法為 y = interpft(yi, n) 向量 yi 是一個經由等距取點的函數值 n 則是等距取點的內插點數,一維內插範例-2 (II),範例9-2：interpft01.m,x = linspace(0, 2*pi, 11); y = sin(x).*exp(-x/5); xi = linspace(0, 2*pi, 21); yi
6、= interpft(y, 21); plot(x, y, o, xi, yi); legend(Original, Curve by interpft),一維內插範例-2 (III),內插曲線必須通過每一個已知的取樣點 迴歸曲線則不需要通過取樣點 內插法適用於雜訊很少的取樣資料,一維內插範例-3 (I),根據平面上有限資料點，來描繪出一個物體的外型，可使用一維內插，此時內插資料點會變成兩組 X座標對其索引（Index，或稱指標）的內插 Y座標對其索引（Index，或稱指標）的內插,一維內插範例-3 (II),七個資料點，散佈在二度空間，我們可以使用 interp1 的 spline 方法，來
7、畫出連續平滑的圖形 範例9-3：interp102.m,x = 0 2 4 3 1 2 1; y = 4 1 1 4 5 2 0; index = 1:length(x); index2 = linspace(1, length(x), 101); x2 = interp1(index, x, index2, spline); y2 = interp1(index, y, index2, spline); plot(x, y, o, x2, y2, -); legend(Origianl data, Interpolated data);,一維內插範例-3 (III),原有的資料七點，使用一維
8、內插後就可以產生平滑連續的圖形 在電腦圖學（Computer Graphics）中最常見的做法 以少數控制點（Control Points）來代表一個物件，然後再使用內插來產生整個物件的細緻圖形,一維內插範例-4 (I),選用不同的方法，得到的平滑效果也會不同 使用 interp1 的四種方法，來對上一個範例的資料點進行內插 範例9-4：interp103.m,一維內插範例-4 (II),x = 0 2 4 3 1 2 1; y = 4 1 1 4 5 2 0; index = 1:length(x); index2 = linspace(1, length(x), 101); x2 = in
9、terp1(index, x, index2, nearest); y2 = interp1(index, y, index2, nearest); x3 = interp1(index, x, index2, linear); y3 = interp1(index, y, index2, linear); x4 = interp1(index, x, index2, spline); y4 = interp1(index, y, index2, spline); x5 = interp1(index, x, index2, cubic); y5 = interp1(index, y, ind
10、ex2, cubic); plot(x, y, o, x2, y2, -); plot(x, y , o, x2 ,y2 ,m- ,x3 ,y3 ,r: ,x4 ,y4 ,k-., x5, y5, b-); legend(Original, Nearest, Linear, Spline, Cubic);,一維內插範例-4 (III),一維內插範例-4 (IV),可讓使用者拖放取樣點，內插曲線就會跟著改變 範例 interpAnim02.m,9-2 二維格子點內插法,interp2 指令可進行二維格子點內插 使用語法為 zi = interp2(x, y, z, xi, yi, method)
11、 z 是一矩陣，代表一函數的高度 矩陣 x 及 y 則是此函數在方格子點的 x 及 y 座標 矩陣 xi 及 yi 則是內插點的 x 及 y 座標 字串 method 則指定使用的內插法 nearest：鄰近點內插法 linear：二維線性內插法 spline：二維 Spline 內插法 cubic：二維三次內插法,二維格子點內插法範例,使用 interp2 時，x 及 y 都必須是嚴格遞增或遞減 x 及 y 都是由 meshgrid 所產生，以保證其格式的正確性 用不同的方法進行二維內插，會得到不同的效果,二維格子點內插法程式(I),範例9-5：interp201.m,x,y = meshg
12、rid(-3:1:3); % 產生 -3 至 3 的格子點 z = peaks(x,y); % 產生 peaks 曲面資料 xi, yi = meshgrid(-3:.25:3); % 生產內插點 zi1 = interp2(x, y, z, xi, yi, nearest); zi2 = interp2(x, y, z, xi, yi, linear); zi3 = interp2(x, y, z, xi, yi, cubic); subplot(2,2,1); surf(x, y, z); axis tight; title(Original); subplot(2,2,2); surf(
13、xi, yi, zi1); axis tight; title(Nearest); subplot(2,2,3); surf(xi, yi, zi2); axis tight; title(Linear); subplot(2,2,4); surf(xi,yi,zi3); axis tight; title(Cubic);,二維格子點內插法程式(II),% 以下畫出等高線 figure subplot(2,2,1); contour(x, y, z, 20); title(Original); subplot(2,2,2); contour(xi, yi, zi1, 20); title(Ne
14、arest); subplot(2,2,3); contour(xi, yi, zi2, 20); title(Linear); subplot(2,2,4); contour(xi, yi, zi3, 20); title(Cubic);,二維格子點內插法圖示,二維格子點內插法圖示,由上兩圖可看出，method = cubic 產生的曲面較其它兩種方法（nearest 及 linear）為平滑,9-3 二維散佈點內插法,對於散佈式資料（Scatter Data），使用 griddata 指令 使用語法為 zi = griddata(x, y, z, xi, yi) x、y、z 是散佈點的座標
15、 xi、yi 則是欲求內插值的點（通常是格子點）,二維散佈點內插法範例,範例9-6：griddata01.m,x = 6*rand(100, 1)-3; % -3, 3 之間的 100 個均勻分佈亂數 y = 6*rand(100, 1)-3; % -3, 3 之間的 100 個均勻分佈亂數 z = peaks(x, y); xi, yi = meshgrid(-3:0.2:3, -3:0.2:3); zi = griddata(x, y, z, xi, yi); mesh(xi, yi, zi); % 畫出曲面 hold on; plot3(x, y, z, o); hold off % 畫
16、出資料點 axis tight; hidden off;,二維散佈點內插法圖示,每一圓球代表取樣點（共 100 點），而曲面則是使用 griddata 指令的內插結果。（由於這100點是隨機產生，每次執行都不一樣，所以得到的曲面也會和上述曲面有所差異。）,二維散佈點內插法圖示,griddata 指令只進行內插，而不進行外插，因此所產生的曲面資料只定義在這些取樣點所形成的最小凸多邊型（Convex Hull） 驗證 : 在執行上述範例後，再輸入 view(2)，改由正上方俯視圖形,提示,在進行 griddata 的內插時，這些散佈的資料點要越多越好（但計算時間也會隨之拉長），而且資料的分佈也是要
17、越平均越好，才能逼近原來正確的 peaks 曲面 對於三維散佈資料點的內插，可以使用 griddata3 指令，對於更高維度的內插，則可以使用 griddatan 指令，其用法都和 griddata 類似，可由線上支援找到相關說明,9-4 三維格子點內插法,interp3指令可進行 3 維格子點的內插 語法為 vi = interp3(x, y, z, v, xi, yi, zi, method) x、y、z、v 是三維矩陣，前三者代表資料點的輸入部份 v 是資料點的輸出部份 xi、yi 及 zi 則是內插點 字串 method 則指定不同的內插方法 nearest：鄰近點內插 linear：
18、線性內插 cubic：三次內插 spline：Spline 內插法,注意,使用 interp3 時，矩陣 x、y、z 必須是嚴格遞增或遞減 x、y 及 z 都是由 ndgrid 指令 所產生，以保證其格式的正確性,三維格子點內插法範例-1 (I),使用 flow 指令產生三度空間中的資料點（這些資料是火箭噴射口的速度大小），並以切面圖（Slice）的方式畫出 範例9-7：slice01.m ans =10 20 10,x, y, z, v = flow(10); slice(x, y, z, v, 6 9.5, 2.5, -2 0); size(x),三維格子點內插法範例-1 (II),x、y
19、、z 的維度是 102010 slice 指令可對三度空間做“切片”，之後經由顏色的不同，就可以區別函數值的高低,提示,在MATLAB中 常用的一維輸入函數是 humps，用來測試曲線的零點或是極值 常用的二維輸入函數是 peaks，用來測試曲面圖或是等高線 常用的三維輸入函數是 flow，用來測試切片圖或三維內差,三維格子點內插法範例-2 (I),使上述範例圖形更為細緻，使用 meshgrid 指令產生三度空間的格子點，並用 interp3 來進行內插 範例9-8：interp301.m ans =25 40 25,x, y, z, v = flow(10); xi, yi, zi = me
20、shgrid(.1:.25:10, -3:.25:3, -3:.25:3); vi = interp3(x, y, z, v, xi, yi, zi); slice(xi, yi, zi, vi, 6 9.5, 2.5, -2 0); % 產生切片圖 colorbar; % 顯示顏色與函數值的對照表 size(xi),三維格子點內插法範例-2 (II),xi、yi、zi 的維度是 254025，比原先分佈的要稠密，畫出來的圖形也比較細緻 若不要顯示格線，可輸入shading interp，所得到的顏色就是進一步經由內插所得到平滑變化的顏色,9-5 高維格子點內插法,interpn 指令可進行高
21、維格子點內插 語法為 vi = interpn(x1, x2, v, y1, y2, method) x1、x2、 是資料點的輸入部份 v 是資料點的輸出部份 y1、y2、 是內插點 字串 method 則可指定內插方法 nearest：鄰近點內插 linear：線性內插 spline：Spline 內插 cubic：三次內插,注意,使用 interpn ，x1、x2、 必須是嚴格遞增或遞減 x1、x2、 最好是由 ndgrid 指令產生，以確保其格式正確性 二維格子點可由 meshgrid 指令產生 高維格子點可由 ndgrid 指令產生 語法為 X1, X2, X3, = ndgrid(x
22、1, x2, x3,) XI 的第 i 維元素是由 xi 反覆產生,高維格子點內插法程式(I),可計算此方程式 在格子點的值，並以顏色代表值的大小，然後再用更密的資料點來進行內插，得到更細緻的圖 範例9-9：interpn01.m,x1 = -2:0.4:2; x2 = -2:0.5:2; x3 = -2:0.3:2; x1, x2, x3 = ndgrid(x1, x2, x3); z = x2.*exp(-x1.2-x2.2-x3.2); subplot(2,1,1); slice(x2, x1, x3, z, -1, 1, , 0); shading interp % 刪除格線，並改用平
23、滑的顏色,高維格子點內插法程式(II),view(-20, 15); colorbar; % 顯示顏色與函數值的對照表 y1 = -2:0.1:2; y2 = -2:0.1:2; y3 = -2:0.1:2; y1, y2, y3 = ndgrid(y1, y2, y3); zz = interpn(x1, x2, x3, z, y1, y2, y3); subplot(2,1,2); slice(y2, y1, y3, zz, -1, 1, , 0); shading interp % 刪除格線，並改用平滑的顏色 view(-20, 15); colorbar; % 顯示顏色與函數值的對照表
24、,高維格子點內插法圖示,上方的圖比較粗糙，下方的圖則比較細密，顯示 interpn 指令的內插效果 使用 interpn 指令來對 3 個輸入的資料點進行內插，此功能和 interp3 相同的 interpn 指令可以對更高維度的資料點進行內插，所以功能更為強大,9-6 三角內插法與計算幾何,MATLAB 提供了一系列指令 ，用來解決 三角化（Triangulation） 鄰近點（Nearest Neighbors） 其它計算幾何（Computational Geometry）,三角內插法範例-1 (I),何謂 Delaunay 三角化（Trangulation） 給定一組 X-Y 平面上的點
25、，Delaunay 指令可傳回一組由這些點所形成的三角形 任一點均不會位於任一個三角形的外接圓之內 展示 Delaunay 的用法，讀入一組 X-Y平面上的資料點，以散佈圖的方式來表示 範例9-10：seamount01.m,load seamount.mat plot(x, y, .);,三角內插法範例-1 (II),三角內插法範例-2 (I),對上述資料點進行 Delaunary 三角化，並使用 triplot 指令畫出其結果，再將原先的資料疊上去 範例9-11：triplot01.m,load seamount.mat plot(x, y, .); tri = delaunay(x, y
26、); triplot(tri, x, y); hold on, plot(x, y, .r); hold off,三角內插法範例-2 (II),三角內插法範例-3 (I),也可以使用 trisurf 或是 trimesh 來畫出使用三角內插所產生的曲面 範例9-12：trisurf01.m把trisurf 改成 trimesh，就可以產生三角網狀曲面,load seamount.mat plot(x, y, .); tri = delaunay(x, y); trisurf(tri, x, y, z); axis tight; colorbar,三角內插法範例-3 (II),三角內插法範例-4
27、 (I),要產生等高線，必須使用 griddata 指令 範例9-13：tricontour01.m,load seamount.mat xi = linspace(min(x), max(x), 50); yi = linspace(min(y), max(y), 50); xi, yi = meshgrid(xi, yi); zi = griddata (x, y, z, xi, yi, cubic); c, h = contourf(xi, yi, zi, b-); clabel (c, h); colorbar; % 顯示顏色與函數值的對照表,三角內插法範例-4 (II),Vorono
28、i 圖形,完成 Delaunary 三角化之後，可以進行兩種搜尋： dsearch 指令可找出一給定點的最鄰近資料點 tsearch 指令可找出一給定點所在之三角形Voronoi 圖形 和 Delaunay 三角化關係非常密切 可將資料點所在的平面畫分成數個多邊形 每一個多邊形只包含一個資料點 給定任一點時，此點和其最鄰近的資料點會落於同一個多邊形內 將 Delaunary 所得的各個三角形中，對每一邊畫出垂直平方線，即可得到 Voronoi 圖形,Voronoi 圖形範例-1,voronoi 指令可用來畫出 Voronoi 圖形 範例9-14：voronoi01.m,load seamoun
29、t.mat voronoi(x, y);,Voronoi 圖形範例-2 (I),將 Voronoi 圖形和 Delaunay 三角形畫在同一張圖 範例9-15：voronoi02.m,x = rand(20,1); y = rand(20,1); voronoi(x, y, b:); tri = delaunay(x, y); hold on h = trimesh(tri, x, y, 0*x); hold off set(h, facecolor, none); axis equal square axis(0 1 0 1);,Voronoi 圖形範例-2 (II),實線是 Delauna
30、y 三角形，虛線則是 Voronoi 圖形,發現,每一條 Voronoi 圖形的線段，都是某一個三角形一邊的中垂線 每一個 Voronoi 圖形所形成的多邊形頂點，都是某一個三角形的外心 每一個 Voronoi 圖形所形成的多邊形，都只有包含一個取樣資料點，代表此取樣資料點的勢力範圍 ， 落於此多邊形之內的任何點，都會離此取樣資料點較近，而離其他取樣資料點較遠,提示,在樣式辨識（Pattern Recognition）學科的觀點來看 Voronoi 圖形就是經由 KNNR (K-Nearest Neighbor Rule) 分類方法在 K 等於1時所得到的結果 依照不同的類別來對 Vornoi
31、 圖形進行著色，就可以得到 KNNR 的決策邊界（Decision Boundary）,計算幾何範例-1 (I),convhull 指令可用於計算一組資料點的最小凸多邊形（Convex Hulls） 計算並顯示 seamount 資料點的最小凸多邊形 範例9-16：convhull01.m,load seamount.mat plot(x, y, .); k = convhull(x, y); hold on, plot(x(k), y(k), r), hold off % 畫出最小凸多邊,計算幾何範例-1 (II),計算幾何範例-2 (I),Inpolygon指令，用於計算一個資料點是否落於
32、一個封閉的多邊形 先產生一個六邊形，再用 inpolygon 找出位於此六邊形內部的資料點，用不同的顏色畫出這兩類資料 範例9-17：inpolygon01.m,theta = linspace(0, 2*pi, 7); xv = cos(theta); yv = sin(theta); x = randn(250,1); y = randn(250,1); in = inpolygon(x, y, xv, yv); plot(xv, yv, b, x(in), y(in), g., x(in), y(in),k.); axis image,計算幾何範例-2 (II),9-7 本章指令彙整,內插法相關的指令可列表如下：,本章指令彙整,計算幾何相關的指令可列表如下：,
