1차원 polyfit을 이용한 2차원 곡면 fitting 방법을 생각해보았습니다.
최근 문젤 하나 풀다가 짠 루틴



f(x,y)=sin(x^2)^3 + 3*cos(y)  (왼쪽 그림)

를

10차 2변수 다항식 (36개 coefficient) 로 근사합니다.

a55x^5y^5 +a54x^5y^4+a53x^5y^3+..............+a00 (오른쪽 그림)


% Polynomial surface fitting

clear all;
clc;
close all;

Nx=101;
Ny=101;
X=[-1:2/(Nx-1):1];
Y=[-1:2/(Ny-1):1];

A=rand(6,6)-0.5;

% example function #1

for k=1:length(X)
    for l=1:length(Y)

        x=X(k);
        y=Y(l);
        
a5=A(6,6)*y^5+A(6,5)*y^4+A(6,4)*y^3+A(6,3)*y^2+A(6,2)*y+A(6,1);
a4=A(5,6)*y^5+A(5,5)*y^4+A(5,4)*y^3+A(5,3)*y^2+A(5,2)*y+A(5,1);
a3=A(4,6)*y^5+A(4,5)*y^4+A(4,4)*y^3+A(4,3)*y^2+A(4,2)*y+A(4,1);
a2=A(3,6)*y^5+A(3,5)*y^4+A(3,4)*y^3+A(3,3)*y^2+A(3,2)*y+A(3,1);
a1=A(2,6)*y^5+A(2,5)*y^4+A(2,4)*y^3+A(2,3)*y^2+A(2,2)*y+A(2,1);
a0=A(1,6)*y^5+A(1,5)*y^4+A(1,4)*y^3+A(1,3)*y^2+A(1,2)*y+A(1,1);

f(k,l)=a5*x^5   + a4* x^4  + a3* x^3  + a2*x^2 +a1*x + a0;

    end;
end;

figure(1); mesh(f);

% example function #2

for k=1:length(X)
    for l=1:length(Y)
    
        x=X(k);
        y=Y(l);
        f(k,l)=sin(x^2)^3+3*cos(y);
    
    end;
end;
figure(2); mesh(f);

% Inverse----------------------------

for l=1:length(Y)
    
    P=polyfit(X',f(:,l),5);
    
    c5(l)=P(1);
    c4(l)=P(2);
    c3(l)=P(3);
    c2(l)=P(4);
    c1(l)=P(5);
    c0(l)=P(6);
    
end;

P5=polyfit(Y,c5,5);
P4=polyfit(Y,c4,5);
P3=polyfit(Y,c3,5);
P2=polyfit(Y,c2,5);
P1=polyfit(Y,c1,5);
P0=polyfit(Y,c0,5);

for l=1:6
B(6,l)=P5(6-l+1);
B(5,l)=P4(6-l+1);
B(4,l)=P3(6-l+1);
B(3,l)=P2(6-l+1);
B(2,l)=P1(6-l+1);
B(1,l)=P0(6-l+1);
end;

for k=1:length(X)
    for l=1:length(Y)

        x=X(k);
        y=Y(l);
        
a5=B(6,6)*y^5+B(6,5)*y^4+B(6,4)*y^3+B(6,3)*y^2+B(6,2)*y+B(6,1);
a4=B(5,6)*y^5+B(5,5)*y^4+B(5,4)*y^3+B(5,3)*y^2+B(5,2)*y+B(5,1);
a3=B(4,6)*y^5+B(4,5)*y^4+B(4,4)*y^3+B(4,3)*y^2+B(4,2)*y+B(4,1);
a2=B(3,6)*y^5+B(3,5)*y^4+B(3,4)*y^3+B(3,3)*y^2+B(3,2)*y+B(3,1);
a1=B(2,6)*y^5+B(2,5)*y^4+B(2,4)*y^3+B(2,3)*y^2+B(2,2)*y+B(2,1);
a0=B(1,6)*y^5+B(1,5)*y^4+B(1,4)*y^3+B(1,3)*y^2+B(1,2)*y+B(1,1);

f(k,l)=a5*x^5   + a4* x^4  + a3* x^3  + a2*x^2 +a1*x + a0;

    end;
end;

figure(3); mesh(f);