/
splineTransformation2D.m
93 lines (78 loc) · 2.82 KB
/
splineTransformation2D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [y,dy] = splineTransformation2D(w,x,varargin)
%
% computes y = dy*w and returns dy = kron(I2,Q2,Q1), where
% Q{i}(:,1) = spline(x(:,1));
%
% if no argumanets are given, the parameters for the identity map are returned.
%
% required inputs are: p (number of spline coefficients), omega and m
%
% see also transformations/contents.m, trafo.m
%==============================================================================
function [y,dy] = splineTransformation2D(w,x,varargin)
% the persitent variable stores the matrix
% Q(x) = kron( I_2 , Q2, Q1 );
% p is number of spline coefficients, m is size of grid, omega is domain
persistent Q p m omega
% p = []; m = []; omega = [];
for k=1:2:length(varargin), % overwrite default parameter
eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;
if nargin==0
help(mfilename)
runMinimalExample;
return;
else
y = mfilename('fullfile');
dy = zeros(2*prod(p),1); % parameterization of identity
if ischar(w), Q = []; w = []; end; % reset Q
if isempty(w), return; end;
end;
if isempty(w) || (size(Q,1) ~= numel(x)) || (size(Q,2) ~= numel(w)),
% it is assumed that x is a cell centered grid, extract xi1 and xi2
dim = size(omega,2)/2
n = numel(x)/dim
if n == prod(m)
q = m;
elseif n == prod(m+1)
q = m+1
else
error('can not handle this grid')
end;
x = reshape(x,[q,2]);
Q1 = getQ1d(omega(1:2),q(1),p(1),x(:,1,1));
Q2 = getQ1d(omega(3:4),q(2),p(2),x(1,:,2));
Q = kron(speye(2),kron(sparse(Q2),sparse(Q1)));
if nargout == 0, return; end;
end;
y = x(:) + Q*w;
dy = Q;
%------------------------------------------------------------------------------
function Q = getQ1d(omega,m,p,xi)
Q = zeros(m,p); xi = reshape(xi,[],1);
for j=1:p,
cj=zeros(p,1); cj(j) = 1;
Q(:,j) = splineInter(cj,omega,xi);
end;
%------------------------------------------------------------------------------
function runMinimalExample
fprintf('%s: minimal example\n',mfilename)
omega = [0,10,0,8]; m = [8,9]; p = [5,6];
w = zeros([p,2]); w(3,3,1) = 0.05; w(3,4,2) = -0.1;
x = getCellCenteredGrid(omega,m);
y = feval(mfilename,w(:),x,'omega',omega,'m',m,'p',p,'Q',[]);
FAIRfigure(1); clf;
plotGrid(x,omega,m,'color','r'); axis image; hold on;
plotGrid(y,omega,m,'color','b'); axis image; hold off;
fctn = @(w) feval(mfilename,w(:),x,'omega',omega,'m',m,'p',p,'Q',[]);
w = w + randn(size(w));
checkDerivative(fctn,w(:),'fig',2);
%==============================================================================