Rank: 186752 based on 0 downloads (last 30 days) and 0 file submitted
photo

Peter

E-mail

Personal Profile:

Boston Area, MA

Professional Interests:

 

Watch this Author's files

 

Comments and Ratings by Peter View all
Updated File Comments Rating
16 Apr 2012 Surface Fitting using gridfit Model 2-d surfaces from scattered data Author: John D'Errico

Oops. In my explanation above P = 1/sqrt(2), not sqrt(2).

13 Apr 2012 Surface Fitting using gridfit Model 2-d surfaces from scattered data Author: John D'Errico

Lanis, John, et al.

This is how to get consistent behavior with different grid sizes.

--------

Save a copy of gridfit.m to gridfit2.m.

In gridfit2.m make the following changes at the locations of the commented-out code (near Lines 550-600):

% Sorry, John D.
% dy1 = dy(j(:)-1)/params.yscale;
% dy2 = dy(j(:))/params.yscale;
dy1 = dy(j(:)-1);
dy2 = dy(j(:));


% Sorry, John D.
% Append the regularizer to the interpolation equations,
% scaling the problem first. Use the 1-norm for speed.
%NA = norm(A,1);
%NR = norm(Areg,1);
%A = [A;Areg*(smoothparam*NA/NR)];
FidelityEquationCount = size(A, 1);
RegularizerEquationCount = nx * (ny - 2) + ny * (nx - 2);
NewSmoothParam = smoothparam * sqrt(FidelityEquationCount / RegularizerEquationCount);
A = [A; Areg * NewSmoothParam];

--------

Run the following script to demonstrate:

x = [0, 1.333333333, 2.666666667, 4, 0, 0, 4, 4, 2];
y = [2, 2, 2.1, 1.8114, 0, 4, 0, 4, 0];
z = [0.5, 0.82, 0.4, 0.7, 0.5, 0.7, 0, 0.7, 0];

% 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

13 Apr 2012 Surface Fitting using gridfit Model 2-d surfaces from scattered data Author: John D'Errico

Lanis, John, et al.

How it works:

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!

12 Apr 2012 Surface Fitting using gridfit Model 2-d surfaces from scattered data Author: John D'Errico

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.

Here is an example.

x = [0, 1.333333333, 2.666666667, 4, 0, 0, 4, 4, 2];
y = [2, 2, 2.1, 1.8114, 0, 4, 0, 4, 0];
z = [0.5, 0.82, 0.4, 0.7, 0.5, 0.7, 0, 0.7, 0];

Smoothness = 10;

x_coarse = linspace(0, 4, 11);
y_coarse = x_coarse;
Index_coarse = find(y_coarse == 1.2);

x_fine = linspace(0, 4, 41);
y_fine = x_fine;
Index_fine = find(y_fine == 1.2);

z_coarse = gridfit(x, y, z, x_coarse, y_coarse, 'Smoothness', Smoothness);
z_fine = gridfit(x, y, z, x_fine, y_fine, 'Smoothness', Smoothness);

hold on
plot(x_coarse, z_coarse(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine(Index_fine, :), 'color', 'red')
hold off

10 Nov 2011 Surface Fitting using gridfit Model 2-d surfaces from scattered data Author: John D'Errico

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
];

% Create a 4x4 grid.
xNodes = linspace(0, 1, 4);
yNodes = xNodes;

% 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)), '.']);

Contact us