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

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

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

一維內插：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');

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

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

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

Piecewise Cubic Hermit Interp.

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

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

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

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

一維內插：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');

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

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

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

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

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

任意曲線的內插：各種方法 (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');

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

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’：二維三次內插法

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

二維格子點內插法：範例 (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');

二維格子點內插法：範例 (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');

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

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

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

二維散佈點內插法：範例 (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;

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

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

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

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 內插法

三維格子點內插法：注意事項 使用 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) slice 指令可對三度空間做“切片”，之後經由顏色的不同，就可以區別函數值的高低 x、y、z 的維度是 10×20×10 slice 指令可對三度空間做“切片”，之後經由顏色的不同，就可以區別函數值的高低 z=[-2, 0] y=[2.5] X=[6, 9.5]

提示：MATLAB常用內建函數 在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] = 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)

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

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 內插

高維格子點內插法：注意事項 使用 interpn 注意事項 提示 x1、x2、… 必須是嚴格遞增或遞減 x1、x2、… 最好是由 ndgrid 指令產生，以確保其格式正確性 提示 二維格子點可由 meshgrid 指令產生 高維格子點可由 ndgrid 指令產生 語法為 [X1, X2, X3,…] = ndgrid(x1, 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 % 刪除格線，並改用平滑的顏色

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

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

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

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

三角內插法：範例-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

三角內插法：範例-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

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

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

Voronoi 圖形：範例-1 voronoi 指令可用來畫出 Voronoi 圖形 範例9-14：voronoi01.m load seamount.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]); 實線是 Delaunay 三角形 虛線則是 Voronoi 圖形

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

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

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

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

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