Summary

We have developed a very fast algorithm for image reconstruction from non-uniform samples. The solution is expressed as a polynomial spline and is computed efficiently using a multigrid solver.

Introduction

Thin-plate splines provide an elegant solution to the problem of interpolating or approximating non-uniform data. The method is optimal in a well-defined variational sense. Unfortunately, it tends to break down numerically when there are too many data points. The main difficulty is that the linear system to be solved is ill-conditioned. Moreover, the matrix is not sparse. Hence, solving the system soon becomes overly expensive.

Main Contribution

We have been successful in developing a much faster algorithm that essentially solves the same approximation problem. In the thin-plate-spline variational formulation, the reconstruction is formulated as the minimizer of a cost that is a weighted sum of two terms: (i) the sum of squared errors at the specified points; (ii) a quadratic regularizing functional that penalizes any lack of smoothness. Here, we discretize the problem in a uniform B-spline basis and show how the solution can be determined by solving a large, sparse system of linear equations. Using the two-scale relation for B-splines, we derive an algebraic relation that links together the linear systems of equations that specify the reconstructions at different levels of resolution. We use this relation to develop an ultra-fast multigrid algorithm. We have implemented the approach and have obtained impressive reconstruction results, both in terms of speed and quality.

Our algorithm has a number of advantageous features that should make it attractive for practical applications:

• The user can select the step size of the reconstruction grid. This parameter controls the tradeoff between computational complexity and reconstruction accuracy. As the step size gets smaller, the reconstruction converges to the solution of the thin-plate-spline problem.
• The algorithm allows for a tradeoff between smoothness and closeness of fit through the adjustment of the regularization parameter. As the regularization factor tends to zero, it tends to interpolate the data.
• The algorithm has the ability to handle arbitrary sample locations.
• It has a complexity that depends primarily on the number of reconstruction grid points. There is essentially no dependence on the location and number of samples, which is unusual in this type of problem.