Two dimensional interpolation

% Solving peak position for 2d grid values
%   by 2nd order polynomial surface fitting
% Coded by funMV. 12/2013
%
% v(x,y)=v(row,col)=a1+a2*x+a3*y+a4*x*y+a5*x^2+a6*y^2
% parameters ai is solved by linear equation system:
% v(-1,-1)=a1-a2-a3+a4+a5+a6
% v(-1, 0)=a1-a2+a5
% v(-1, 1)=a1-a2+a3-a4+a5+a6
% v( 0,-1)=a1-a3+a6
% v( 0, 0)=a1
% v( 0, 1)=a1+a3+a6
% v( 1,-1)=a1+a2-a3-a4+a5+a6
% v( 1, 0)=a1+a2+a5
% v( 1, 1)=a1+a2+a3+a4+a5+a6 
%
% Max value is at extreme position:
%   dv(x,y)/dx = 0
%   dv(x,y)/dy = 0
%
%Two dimensional interpolation for finding max value on sub-pixel. Given values are on %grid and 2nd order polynomial surface is fit to grid values.

clear all; clc;

A=[1 -1 -1  1  1  1;
   1 -1  0  0  1  0;
   1 -1  1 -1  1  1;
   1  0 -1  0  0  1;
   1  0  0  0  0  0;
   1  0  1  0  0  1;
   1  1 -1 -1  1  1;
   1  1  0  0  1  0;
   1  1  1  1  1  1];
pinv = inv(A'*A)*A';  % 미리 계산해 둘 수 있음

b=[21 22 21   24 25 24   21 22 21]'; % 9개 격자 값. (중앙 값이 최대)
x= pinv*b;

%row = -1; col = 1;
%c1 = [1 row col row*col row*row col*col];
%val = c1*x;

m1 = [2*x(5) x(4); x(4) 2*x(6)];
b1 = [x(2); x(3)];

c2 = -inv(m1)*b1;
c2




References
[1] 김성완, 김남식, Digital Image Correlation 기법을 이용한 구조물의 다중 동적변위 응답측정, 한국지진공학회 논문집, 13권 3호, 2009.