Re: [討論] fit平面問題

[tn00364361 (小氵斿)]

※ 引述《liwes5566 (wes5566)》之銘言：
: 目前是利用指令中的fit(2元2次多項式poly22)下去fit
: 再一般的情況下，都可以fit得很好
: 可是當我是fit一個平面(數值全部都一樣)，所得到的答案不是平面
: 這是我已1024*1024矩陣，數值為2^16的結果
: fitobject(x,y)=p00+p10x+p01y+p20x^2+p11xy+p02y^2
: p00=6.554e+04
: p10=-5.012e-12
: p01=3.291e-09
: p20=2.868e-14
: p11=-7.413e-16
: p02=-1.801e-11
: 有沒有甚麼方法可以讓p00的值為2^16，其他都為0呢?

任何方法fit的結果都不太可能會是真的零誤差

在一般電腦上MATLAB的machine epsilon大概是2.22e-16

通常可接受的數值誤差範圍大概是sqrt(machine epsilon)，也就是1.49e-8

----------------------------------------

Fit多項式曲面(或是曲線)可以看成是一個線性的least-squares問題

假設你有m個資料點，p = [p1; p2; ...; pn]為你想要fit的曲面的係數。

而曲面的方程式可以寫成 z = a1 * p1 + a2 * p2 + ... + an * pn = dot(a, p)

其中a是x跟y的函數

把所有的資料點"疊"起來，可以寫成
z1 ~ a11 * p1 + a12 * p2 + ... + a1n * pn
z2 ~ a21 * p1 + a22 * p2 + ... + a2n * pn
...
zm ~ am1 * p1 + am2 * p2 + ... + amn * pn

這可以重新寫成矩陣的形式: Z ~ A * p (這裡的"~"是指大約等於)

其中，[A]ij = aij，Z = [z1; z2; ...; zm]

這時候的最佳解(依least-squares)就是 p_fit = argmin(norm(A * p - Z))

也就是 p_fit = A \ Z (或是寫成 p_fit = pinv(A) * Z)

----------------------------------------

回到你的case。先把格點跟資料點都向量化，變成column vector
n = 1024;
x = reshape(x, n^2, 1); % grid point
y = reshape(y, n^2, 1); % grid point
z = reshape(z, n^2, 1); % data point

----------------------------------------

1. 如果用least-squares去fit所有的係數:
A = [ones(n^2, 1), x, y, x.^2, x .* y, y.^2];
p_fit = A \ z;  % = [p00; p10; p01; p20; p11; p02]

----------------------------------------

2. 如果只去fit除了p00以外的係數(已知p00 = 2^16):
z = z - 2^16;
A = [x, y, x.^2, x .* y, y.^2];
p_fit = A \ z; % = [p10; p01; p20; p11; p02]

----------------------------------------
這個方法可以很輕易的fit任意維度空間的任意多項式曲面(例如在三維空間中fit橢球)

希望這對你有幫助

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 68.43.178.167
※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1439899632.A.F90.html
※ 編輯: tn00364361 (68.43.178.167), 08/18/2015 22:09:36
推 liwes5566: 了解! 不過不過要處理的不只全平面，前面我還要再加個
→ liwes5566: 判斷，判斷其是不是為uniform平面
→ tn00364361: uniform是指？
→ liwes5566: 就是全部數值都相同的surface
→ tn00364361: 喔喔，那你要fit的面有雜訊嗎？
→ liwes5566: 我fit的目的就是要計算Noise Power Spectrum(NPS)，應
→ liwes5566: 該都有雜訊，只是遇到全部數值相同的影像就變這樣了