% Pass a test point through kernel splines to give
% comparative data for testing the ITK implementation
%
% Rupert Brooks, March 2007
d=3;  % set the dimension here to generate cases for 2,3 dimensions
mode=0; % 0 thin plate spline
        % 1 thin plate r2 logr
        % 2 elastic
        % 3 elastic reciprocal
        % 4 volume


% generate points 2D case
if(d==2)
testpt=[0.5 0.6];
pts_source=zeros(4,2);
pts_target=zeros(4,2);
temp=0.1;
step=0.01;
count=1;
for i=0:1
    for j=0:1
        pts_source(count,1)=j+temp;
        temp=temp+step;
        pts_source(count,2)=i+temp;
        temp=temp+count*step;
        pts_target(count,1)=j*3;
        pts_target(count,2)=i*3;
        count=count+1;
    end
end

elseif(d==3)
% generate points 3D case
testpt=[0.5 0.6 0.7];
pts_source=zeros(4,3);
pts_target=zeros(4,3);
temp=0.1;
step=0.01;
count=1;
for i=0:1
    for j=0:1
        for k=0:1
        pts_source(count,1)=k+temp;
        temp=temp+step;
        pts_source(count,2)=j+temp;
        temp=temp+step;
        pts_source(count,3)=i+temp;
        temp=temp+count*step;
        pts_target(count,1)=k*3;
        pts_target(count,2)=j*3;
        pts_target(count,3)=i*3;
        count=count+1;
    end
end
end
end

[p,J]=kernelspline(pts_source,pts_target,testpt,mode)

dim=size(pts_source,2);
jdim=dim*2^dim;
ptname=['answerPoint' num2str(dim) 'D'];
disp('//-------------------------------------------------------');
disp('// this code auto-generated using kernelspline_test.m');
disp('// answer point');
for i=1:dim
    disp([ptname '[' num2str(i-1) ']=' num2str(p(i),16) ';']);
end
disp('// answer jacobian');
for i=1:dim
for j=1:jdim
    disp(['answerJacobian[' num2str(i-1) '][' num2str(j-1) ']=' num2str(J(j,i),16) ';']);
end
end
disp('// End of Autogenerated code');
disp('//-------------------------------------------------------');