% Set the smoothness for the original GridFit. This smoothness produces inconsistent results if the grid size changes.
Smoothness = 10;
% Set up the grids with different spacings...
x_coarse = linspace(0, 4, 11);
y_coarse = x_coarse;
x_fine = linspace(0, 4, 41);
y_fine = x_fine;
% We're taking a slice through each surface at the same coordinate.
% They have different indexes but the same coordinate of 1.2.
Index_coarse = find(y_coarse == 1.2);
Index_fine = find(y_fine == 1.2);
% Run GridFit.
z_coarse = gridfit(x, y, z, x_coarse, y_coarse, 'Smoothness', Smoothness);
z_fine = gridfit(x, y, z, x_fine, y_fine, 'Smoothness', Smoothness);
% Try the same thing with the updated GridFit. Note that the smoothness parameter is on a different scale compared to
% the original GridFit. However, all choices of smoothness value should produce consistent results for any grid size.
UpdatedGridFitSmoothness = 1;
% Run the updated version of GridFit, GridFit2.
z_coarse2 = gridfit2(x, y, z, x_coarse, y_coarse, 'interp', 'bilinear', 'Smoothness', UpdatedGridFitSmoothness);
z_fine2 = gridfit2(x, y, z, x_fine, y_fine, 'interp', 'bilinear', 'Smoothness', UpdatedGridFitSmoothness);
% Show the results for GridFit...
% These lines only line up when you use a very large or very small smoothness value.
subplot(1, 2, 1)
hold on
plot(x_coarse, z_coarse(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine(Index_fine, :), 'color', 'red')
title('Profiles from the Original GridFit');
axis([1, 4, 0.3, 0.6]);
hold off
% Show the results for GridFit2...
% These lines have the same profile for all smoothness values.
subplot(1, 2, 2)
hold on
plot(x_coarse, z_coarse2(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine2(Index_fine, :), 'color', 'red')
title('Profiles from the Updated GridFit');
axis([1, 4, 0.3, 0.6]);
hold off
GridFit is a function that regularizes a dataset. It constructs a series of linear equations that are the "fidelity equations" (how closely the output matches your input data) and the "smoothness equations" (how close the surface's second derivatives are to zero). It assembles these equations and solves them in a way that minimizes the sum of the squared error (SSE) of the output surface.
The key is to recognize that the SSE is just a sum of squared residuals, the fidelity squared residuals (SSE_F) and the smoothness squared residuals (SSE_Sm).
SSE = SSE_F + SSE_Sm
If we have a certain number of input points, that tells us how many terms are in SSE_F. The number of grid points (more precisely, the number of second derivatives within the grid) tells us how many terms are in SSE_Sm.
Imagine that we are trying to fit this parabola: y = x^2 to our data. The fidelity terms are some finite value (which doesn't matter for this example), such as 100. Because we're using this parabola, the second derivative is always 2. If we have a grid such that there are 20 smoothness terms, then SSE_Sm = 20 * (2^2) = 80.
The ratio SSE_F / SSE_Sm = 100 / 80 = 1.25. This ratio determines how to weigh the smoothness against the goodness of fit.
Now suppose that we suddenly double the number of grid points. Now SSE_Sm = 2 * 20 * (2^2) = 160. That means SSE_F / SSE_Sm = 100 / 160 = 0.625. Uh-Oh! The ratio has changed, which means our fit will favor smoothness more than it did before!
The way to fix this is to reduce the values of the second derivatives by multiplying them by a constant, which we'll call P. In this example, if P = sqrt(2), then the individual terms in SSE_Sm will have half of the value they had before. That's ok because there are twice as many of them. In other words, 0.5 * 2 = 1.
Now using P, SSE_Sm = 2 * 20 * ((P * 2) ^ 2) = 80. This provides the same balance of fidelity vs. smoothness that we originally had.
What I did in GridFit2 is a little more complicated because it also considers the smoothness parameter and the number of fidelity equations. I hope this helps!
I'm not convinced that the smoothness parameter produces consistent results for different grid sizes. I think the problem comes from the use of the 1-norm instead of just using the numbers of fidelity equations vs. derivatives.
GridFit is pretty good, but I have a few questions.
1) The interpolation seems to be triangular, which causes some symmetry problems. Below is an example program to illustrate this. Does anybody know if there is an advantage to using triangular interpolation (three points) instead of rectangular interpolation (four points)?
2) If I change the grid spacing, I also have to adjust the smoothness parameter or GridFit gives different results. This appears to be related to the way GridFit sets up its derivative equations. It mixes horizontal and vertical second derivatives in the same equation. For example, if the residual in the horizontal direction is +0.1 and the residual in the vertical direction is -0.1, the regularizer wrongly deems that point a "perfect fit" because the derivatives are added together. Would it be better to separate the derivatives into separate equations (separate rows in Areg)? If so, is it possible to produce the same results (same smoothness or curve profile) for a single "smoothness" parameter, no matter what grid spacing is used? I think that would be really helpful and convenient to GridFit users so that we wouldn't have to fiddle with the smoothness number manually.
% Choose four points at each of four corners and one point in the center.
% Make sure everything is symmetric around the center.
SourceData = [
0, 0, 0
0, 1, 0
1, 0, 0
1, 1, 0
0.5, 0.5, 1
];
% Run GridFit. The smoothness parameter has little effect on the result
% for this demonstration.
[zGrid, xGrid, yGrid] = gridfit(SourceData(:, 1), SourceData(:, 2), SourceData(:, 3), xNodes, yNodes, 'smoothness', 0.1);
% Display the surface so that we can see the asymmetry.
surf(xGrid, yGrid, zGrid);
xlabel('x');
ylabel('y');
zlabel('z');
% Calculate the asymmetry.
disp(['GridFit asymmetry is ', num2str(zGrid(2, 2) - zGrid(3, 2)), '.']);
