張智星 jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室 MATLAB 程式設計進階篇 內插法 張智星 jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室.

Presentation on theme: "張智星 jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室 MATLAB 程式設計進階篇 內插法 張智星 jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室."— Presentation transcript:

1 張智星 jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室
MATLAB 程式設計進階篇 內插法 張智星 台大資工系 多媒體檢索實驗室

2 9-1 一維內插法 一維內插（1-D interpolation） MATLAB 提供了兩種一維內插的基本方法：
9-1 一維內插法 一維內插（1-D interpolation） 根據一組已知的資料點（包含輸入及輸出，其中輸入是一維資料，輸出也是一維資料），建立一個連續的函數 算出任意輸入資料點所對應的輸出值 MATLAB 提供了兩種一維內插的基本方法： 基於多項式的內插法 Nearest, linear, spline, pchip, … 基於FFT（Fast Fourier Transform，快速傅立葉轉換）的內插法

3 一維內插：interp1 指令 interp1指令即是利用多項式來進行內插運算
使用語法為 yi = interp1(x, y, xi, method) 向量 x 是資料點的 x 座標（輸入值） 向量 y 是資料點的 y 座標（輸出值） 向量 xi 是內插點（輸入值，未知其對應輸出值） 字串 method 則指定使用的方法 ‘nearest’：鄰近點內插法 ‘linear’：線性內插法 ‘spline’：片段式的三次 Spline 內插法 ‘pchip’：保持形狀的片段式三次內插 ‘cubic’：和‘pchip’ 一樣 'v5cubic'：MATLAB 5 所用的三次內插法

4 一維內插：interp1 範例程式碼 範例9-1：interp101.m x = 0:1:4*pi;
y = sin(x).*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, 'pchip'); y4 = interp1(x, y, xi, 'spline'); plot(x, y, 'o', xi, y1, xi, y2, xi, y3, xi, y4); legend('Original', 'Nearest', 'Linear', 'Pchip', 'Spline');

5 一維內插：interp1 範例結果 由圖可看出，Spline 和 Pchip 所產生的曲線較平滑，但它們所需的計算時間也較久

6 四種內插法比較表 執行時間 1（短） 2 3 4（長） 曲線平滑度 1（差） 4（好） 記憶體使用量 1（少） 4（多） Nearest
Linear Pchip Spline 執行時間 1（短） 2 3 4（長） 曲線平滑度 1（差） 4（好） 記憶體使用量 1（少） 4（多）

7 Spline Spline: 雲尺、雲規、曲線板
Based on 3-order polynomials First and second derivatives are continuous Pchip: piecewise cubic Hermite interpolating polynomial First derivatives are continuous Local max. & min. are preserved

8 Piecewise Cubic Hermit Interp.

9 Formulas for Spline & Pchip
… Reference “Numerical Computing with MATLAB” by Cleve Moler

10 使用interp1指令：注意事項 使用 interp1指令 向量 x 必須是嚴格遞增或遞減 元素之間不必是等間隔
xi 的範圍必須落在 x 的範圍之內 如果 xi 的範圍在 x 的範圍之外， 可以使用yi = interp1(x, y, xi, method, ‘extrap’) 的方式來進行外插（Extrapolation），以求得在範圍外的 xi 元素所對應的 yi 值

11 一維內插：互動範例展示 可以讓使用者拖放取樣點，內插曲線就會跟著改變 範例 interpAnim01.m
使用者可以在上述圖形視窗中，拖放每一個取樣點，就可以立刻看到內插曲線的變化，也同時可以知道各種內插方法的特性

12 一維內插：interpft 指令 interpft指令，可進行基於 FFT（Fast Fourier Transform，快速傅立葉轉換）的內插法 先計算資料點的傅立葉轉換，再用更密集的內插點來進行反傅立葉轉換 使用語法為 y = interpft(yi, n) 向量 yi 是一個經由等距取點的函數值 n 則是等距取點的內插點數

13 一維內插：interpft 範例 範例9-2：interpft01.m
n = 11; % Number of the original data points factor = 4; % Increase the data by this factor x = linspace(0, 2*pi, n); y = sin(x).*exp(-x/5); xi = (0:factor*n-1)*(x(2)-x(1))/factor; yi = interpft(y, factor*n); plot(x, y, 'ro', xi, yi, '.-'); legend('Original', 'Curve by interpft');

14 內插與迴歸的差異 內插（Interpolation）與迴歸（Regression）的差異 因此： 內插曲線必須通過每一個已知的取樣點
迴歸曲線則不需要通過取樣點 因此： 內插法適用於雜訊很少的取樣資料，如電腦繪圖 迴歸法適用於雜訊較大的取樣資料，如氣象、股票。

15 任意曲線的內插 若要根據二維平面上任意資料點，來描繪出一個物體的外型，我們仍可使用一維內插，此時內插資料點會變成兩組
X座標對其索引（Index，或稱指標）的內插 Y座標對其索引（Index，或稱指標）的內插

16 任意曲線的內插範例：interp1 七個資料點，散佈在二度空間，我們可以使用 interp1 的 spline 方法，來畫出連續平滑的圖形
範例9-3：interp102.m x = [ ]; y = [ ]; 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');

17 任意曲線的內插：interp1 在上例中，原有的資料七點，使用一維內插後就可以產生平滑連續的圖形
這是在電腦圖學（Computer Graphics）中最常見的做法，亦即以少數控制點（Control Points）來代表一個物件，然後再使用內插來產生整個物件的細緻圖形

18 任意曲線的內插：各種方法 (I) 選用不同的方法，得到的平滑效果也會不同 使用 interp1 的四種方法，來對上一個範例的資料點進行內插
範例9-4：interp103.m

19 任意曲線的內插：各種方法 (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 = interp1(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, 'pchip'); y4 = interp1(index, y, index2, 'pchip'); x5 = interp1(index, x, index2, 'spline'); y5 = interp1(index, y, index2, 'spline'); plot(x, y, 'o', x2, y2, '-'); plot(x, y , 'o', x2 ,y2, x3 ,y3 ,x4 ,y4, x5, y5); legend('Original', 'Nearest', 'Linear', 'Pchip', 'Spline');

20 任意曲線的內插：互動展示 可讓使用者拖放取樣點，內插曲線就會跟著改變 範例 interpAnim02.m

21 9-2 二維格子點內插法 interp2 指令可進行二維格子點內插
9-2 二維格子點內插法 interp2 指令可進行二維格子點內插 使用語法為 zi = interp2(x, y, z, xi, yi, method) z 是一矩陣，代表一函數的高度 矩陣 x 及 y 則是此函數在方格子點的 x 及 y 座標 矩陣 xi 及 yi 則是內插點的 x 及 y 座標 字串 method 則指定使用的內插法 ‘nearest’：鄰近點內插法 ‘linear’：二維線性內插法 ‘spline’：二維 Spline 內插法 ‘cubic’：二維三次內插法

22 二維格子點內插：注意事項 使用 interp2 之注意事項 x 及 y 都必須是嚴格遞增或遞減
x 及 y 都是由 meshgrid 所產生，以保證其格式的正確性 用不同的方法進行二維內插，會得到不同的效果

23 二維格子點內插法：範例 (I) 範例9-5：interp201.m
[x,y] = meshgrid(-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'); zi4 = interp2(x, y, z, xi, yi, 'spline'); subplot(2,3,1); surf(x, y, z); axis tight; title('Original'); subplot(2,3,3); surf(xi, yi, zi1); axis tight; title('Nearest'); subplot(2,3,4); surf(xi, yi, zi2); axis tight; title('Linear'); subplot(2,3,5); surf(xi, yi, zi3); axis tight; title('Cubic'); subplot(2,3,6); surf(xi, yi, zi4); axis tight; title('Spline');

24 二維格子點內插法：範例 (II) % 以下畫出等高線 figure
subplot(2,3,1); contour(x, y, z, 20); title('Original'); subplot(2,3,3); contour(xi, yi, zi1, 20); title('Nearest'); subplot(2,3,4); contour(xi, yi, zi2, 20); title('Linear'); subplot(2,3,5); contour(xi, yi, zi3, 20); title('Cubic'); subplot(2,3,6); contour(xi, yi, zi4, 20); title('Spline');

25 二維格子點內插法：範例 (III)

26 二維格子點內插法：範例 (IV) 由上兩圖可看出，'cubic' and 'spline' 產生的曲面較其它兩種方法（'nearest' 及 'linear'）為平滑

27 9-3 二維散佈點內插法 對於散佈式資料（Scatter Data），我們可以使用 griddata 指令
9-3 二維散佈點內插法 對於散佈式資料（Scatter Data），我們可以使用 griddata 指令 使用語法為 zi = griddata(x, y, z, xi, yi) x、y、z 是散佈點的座標 xi、yi 則是欲求內插值的點（通常是格子點）

28 二維散佈點內插法：範例 (I) 範例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 % 畫出資料點 axis tight; hidden off;

29 二維散佈點內插法：範例 (II) 每一圓球代表取樣點（共 100 點），而曲面則是使用 griddata 指令的內插結果。（由於這100點是隨機產生，每次執行都不一樣，所以得到的曲面也會和上述曲面有所差異。）

30 二維散佈點內插法：範例 (III) griddata 指令只進行內插，而不進行外插，因此所產生的曲面資料只定義在這些取樣點所形成的最小凸多邊型（Convex Hull） 驗證 : 在執行上述範例後，再輸入 view(2)，改由正上方俯視圖形

31 二維散佈點內插法：注意事項 在進行 griddata 的內插時，這些散佈的資料點要越多越好（但計算時間也會隨之拉長），而且資料的分佈也是要越平均越好，才能逼近原來正確的 peaks 曲面 對於三維散佈資料點的內插，可以使用 griddata3 指令，對於更高維度的內插，則可以使用 griddatan 指令，其用法都和 griddata 類似，可由線上支援找到相關說明

32 9-4 三維格子點內插法 interp3指令可進行 3 維格子點的內插
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’：線性內插 ‘cubic’：三次內插 'spline'：Spline 內插法

33 三維格子點內插法：注意事項 使用 interp3 時的注意事項： 矩陣 x、y、z 必須是嚴格遞增或遞減
因此 x、y 及 z 最好都是由 ndgrid 指令 所產生，以保證其格式的正確性

34 三維格子點內插法：範例-1 (I) 使用 flow 指令產生三度空間中的資料點（這些資料是火箭噴射口的速度大小），並以切片圖（Slice）的方式畫出 範例9-7：slice01.m ans = [x, y, z, v] = flow(10); slice(x, y, z, v, [6 9.5], 2.5, [-2 0]); size(x)

35 三維格子點內插法：範例-1 (II) slice 指令可對三度空間做“切片”，之後經由顏色的不同，就可以區別函數值的高低
x、y、z 的維度是 10×20×10 slice 指令可對三度空間做“切片”，之後經由顏色的不同，就可以區別函數值的高低 z=[-2, 0] y=[2.5] X=[6, 9.5]

36 提示：MATLAB常用內建函數 在MATLAB中 常用的一輸入函數是 humps，用來測試曲線的零點或是極值
常用的二輸入函數是 peaks，用來測試曲面圖或是等高線 常用的三輸入函數是 flow，用來測試切片圖或三維內插

37 三維格子點內插法：範例-2 (I) 使上述範例圖形更為細緻，使用 meshgrid 指令產生三度空間的格子點，並用 interp3 來進行內插 範例9-8：interp301.m ans = [x, y, z, v] = flow(10); [xi, yi, zi] = meshgrid(.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)

38 三維格子點內插法：範例-2 (II) xi、yi、zi 的維度是 25×40×25，比原先分佈的要稠密，畫出來的圖形也比較細緻
原始資料 xi、yi、zi 的維度是 25×40×25，比原先分佈的要稠密，畫出來的圖形也比較細緻 若不要顯示格線，可輸入「shading interp」，所得到的顏色就是進一步經由內插所得到平滑變化的顏色 內插後資料

39 9-5 高維格子點內插法 interpn 指令可進行高維格子點內插
9-5 高維格子點內插法 interpn 指令可進行高維格子點內插 語法為 vi = interpn(x1, x2,.., v, y1, y2,.., method) x1、x2、… 是資料點的輸入部份 v 是資料點的輸出部份 y1、y2、… 是內插點 字串 method 則可指定內插方法 ‘nearest’：鄰近點內插 ‘linear’：線性內插 'cubic'：三次內插 ‘spline’：Spline 內插

40 高維格子點內插法：注意事項 使用 interpn 注意事項 提示 x1、x2、… 必須是嚴格遞增或遞減
x1、x2、… 最好是由 ndgrid 指令產生，以確保其格式正確性 提示 二維格子點可由 meshgrid 指令產生 高維格子點可由 ndgrid 指令產生 語法為 [X1, X2, X3,…] = ndgrid(x1, x2, x3,…) Xi 的第 i 維元素是由 xi 反覆產生

41 高維格子點內插法：範例 (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 % 刪除格線，並改用平滑的顏色

42 高維格子點內插法：範例 (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 % 刪除格線，並改用平滑的顏色

43 高維格子點內插法：範例 (III) 上方的圖比較粗糙，下方的圖則比較細密，顯示 interpn 指令的內插效果
使用 interpn 指令來對 3 個輸入的資料點進行內插，此功能和 interp3 相同 interpn 指令可以對更高維度的資料點進行內插，所以功能更為強大 原始資料 內插後資料

44 9-6 三角內插法與計算幾何 MATLAB 提供了一系列指令 ，用來解決下列問題： 三角化（Triangulation）
9-6 三角內插法與計算幾何 MATLAB 提供了一系列指令 ，用來解決下列問題： 三角化（Triangulation） 鄰近點（Nearest Neighbors） 計算幾何（Computational Geometry）

45 三角內插法：範例-1 (I) Delaunay 三角化（Triangulation）
欲展示 delaunay 的用法，我們先讀入一組 X-Y平面上的資料點，以散佈圖的方式來顯示其分佈： 範例9-10：seamount01.m load seamount.mat plot(x, y, '.');

46 三角內插法：範例-2 (I) 對上述資料點進行 Delaunary 三角化，並使用 triplot 指令畫出其結果，再將原先的資料疊上去
範例9-11：triplot01.m load seamount.mat tri = delaunay(x, y); triplot(tri, x, y); hold on, plot(x, y, '.r'); hold off

47 三角內插法：範例-3 (I) 也可以使用 trisurf 或是 trimesh 來畫出使用三角內插所產生的曲面
範例9-12：trisurf01.m 把trisurf 改成 trimesh，就可以產生三角網狀曲面 load seamount.mat tri = delaunay(x, y); trisurf(tri, x, y, z); axis tight; colorbar

48 三角內插法：範例-4 (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; % 顯示顏色與函數值的對照表

49 Voronoi 圖形 完成 Delaunary 三角化之後，可以進行兩種搜尋： Voronoi 圖形
dsearch 指令可找出一給定點的最鄰近資料點 tsearch 指令可找出一給定點所在之三角形 Voronoi 圖形 可將資料點所在的平面畫分成數個多邊形 每一個多邊形只包含一個資料點 給定任一點時，此點和其最鄰近的資料點會落於同一個多邊形內 將 Delaunary 所得的各個三角形中，對每一邊畫出垂直平方線，即可得到 Voronoi 圖形

50 Voronoi 圖形：範例-1 voronoi 指令可用來畫出 Voronoi 圖形 範例9-14：voronoi01.m
load seamount.mat voronoi(x, y);

51 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([ ]); 實線是 Delaunay 三角形 虛線則是 Voronoi 圖形

52 Voronoi 圖形特性 Voronoi 圖形特性 從樣式辨識（Pattern Recognition）的觀點來看
Voronoi 圖形就是經由 KNNR (K-Nearest Neighbor Rule) 分類方法在 K 等於1時所得到的結果 依照不同的類別來對 Vornoi 圖形進行著色，就可以得到 KNNR 的決策邊界（Decision Boundary）

53 計算幾何：範例-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

54 計算幾何：範例-2 (I) 我們可用 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

55 計算幾何：範例-2 (II)

56 9-7 本章指令彙整 內插法相關的指令可列表如下： 函數 功能 interp1 一維內插法（查表法） interpft
一維內插法（快速傅立葉轉換法） interp2 二維內插法（查表法） interp3 三維內插法（查表法） interpn n 維內插法（查表法） griddata 產生二維格子點 griddata3 產生三維格子點 griddatan 產生 n 維格子點 spline 三次 Spline 內插法

57 本章指令彙整 計算幾何相關的指令可列表如下： 指令 功能 delaunay 執行 Delaunay 三角化 delaunayn
執行 N-D Delaunay 三角化（第六版以後才支援） dsearch 搜尋 Delaunay 三角化的最近點 tsearch 搜尋給定一點所在之三角形 convhull 計算最小凸多邊形（Convex Hull） convhulln 計算 n 維空間中的最小凸多邊形（Convex Hull） voronoi 計算 Voronoi 圖形 voronoin 計算 n 維空間中的 Voronoi 圖形（第六版以後才支援） inpolygon 決定一點是否位於多邊形之內部 rectint 兩個（或多個） 長方形的交集面積 polyarea 多邊形的面積