SURFACE_FITTING
Weighted least-squares fitting by tensor product B-splines to discrete two-dimensional data is performed. Constraints on the spline or its partial derivatives are optional. The spline function
its derivatives, or the square root of its variance function are evaluated after the fitting.
Required Arguments
DATA = DATA(1:4,:) (Input/Output)
An assumed-shape array with size(data,1) = 4. The data are placed in the array:
data(1,i) = xi ,
data(2,i) = yi ,
data(3,i) = zi ,
data(4,i) = σi , i = 1,…ndata .
If the variances are not known, but are proportional to an unknown value, use
data(4,i) = 1,i = 1,…ndata .
KNOTSX = KNOTSX (Input)
A derived type, ?_SPLINE_KNOTS, that defines the degree of the spline and the breakpoints for the data fitting domain, in the first dimension.
KNOTSY = KNOTSY (Input)
A derived type, ?_SPLINE_KNOTS, that defines the degree of the spline and the breakpoints for the data fitting domain, in the second dimension.
Optional Arguments
CONSTRAINTS = SURFACE_CONSTRAINTS (Input)
A rank-1 array of derived type ?_SURFACE_CONSTRAINTS that defines constraints the tensor product spline is to satisfy.
COVARIANCE = G (Output)
An assumed-shape rank-2 array of the same precision as the data. This output is the covariance matrix of the coefficients. It is optionally used to evaluate the square root of the variance function.
IOPT = IOPT(:) (Input/Output)
Derived type array with the same precision as the input array; used for passing optional data to SURFACE_FITTING. The options are as follows:
|
Packaged Options for SURFACE_FITTING |
||
|
Prefix = None |
Option Name |
Option Value |
|
|
SURFACE_FITTING_SMALLNESS |
1 |
|
|
SURFACE_FITTING_FLATNESS |
2 |
|
|
SURFACE_FITTING_TOL_EQUAL |
3 |
|
|
SURFACE_FITTING_TOL_LEAST |
4 |
|
|
SURFACE_FITTING_RESIDUALS |
5 |
|
|
SURFACE_FITTING_PRINT |
6 |
|
|
SURFACE_FITTING_THINNESS |
7 |
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_SMALLNESS, ?_VALUE)
This resets the square root of the regularizing parameter multiplying the squared integral of the unknown function. The argument ?_value is replaced by the default value. The default is ?_value = 0.
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_FLATNESS, ?_VALUE)
This resets the square root of the regularizing parameter multiplying the squared integral of the partial derivatives of the unknown function. The argument ?_VALUE is replaced by the default value.
The default is ?_VALUE = SQRT(EPSILON(?_VALUE))*SIZE, where
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_TOL_EQUAL, ?_VALUE)
This resets the value for determining that equality constraint equations are rank-deficient. The default is ?_VALUE = 10-4.
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_TOL_LEAST, ?_VALUE)
This resets the value for determining that least-squares equations are rank-deficient. The default is ?_VALUE = 10-4.
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_RESIDUALS, DUMMY)
This option returns the residuals = surface - data, in data(4,:). That row of the array is overwritten by the residuals. The data is returned in the order of cell processing order, or left-to-right in x and then increasing in y. The allocation of a temporary for data(1:4,:) is avoided, which may be desirable for problems with large amounts of data. The default is to not evaluate the residuals and to leave data(1:4,:) as input.
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_PRINT, DUMMY)
This option prints the knots or breakpoints for x and y, and the count of data points in cell processing order. The default is to not print these arrays.
IOPT(IO) = ?_OPTIONS(SURFACE_FITTING_THINNESS, ?_VALUE)
This resets the square root of the regularizing parameter multiplying the squared integral of the second partial derivatives of the unknown function. The argument ?_VALUE is replaced by the default value. The default is ?_VALUE = 10-3 × SIZE, where
FORTRAN 90 Interface
Generic: CALL SURFACE_FITTING (DATA, KNOTSX, KNOTSX, KNOTSY[, …])
Specific: The specific interface names are S_SURFACE_FITTING and D_SURFACE_FITTING.
Description
The coefficients are obtained by solving a least-squares system of linear algebraic equations, subject to linear equality and inequality constraints. The system is the result of the weighted data equations and regularization. If there are no constraints, the solution is computed using a banded least-squares solver. Details are found in Hanson (1995).
Fatal and Terminal Error Messages
See the messages.gls file for error messages for SURFACE_FITTING. These error messages are numbered 1151-1152, 1161-1162, 1370-1393.
Examples
Example 1: Tensor Product Spline Fitting of Data
The function
is least-squares fit by a tensor product of cubic splines on the square
There are ndata random pairs of values for the independent variables. Each datum is given unit uncertainty. The grid of knots in both x and y dimensions are equally spaced, in the interior cells, and identical to each other. After the coefficients are computed a check is made that the surface approximately agrees with g(x,y) at a tensor product grid of equally spaced values.
USE surface_fitting_int
USE rand_int
USE norm_int
implicit none
! This is Example 1 for SURFACE_FITTING, tensor product
! B-splines approximation. Use the function
! exp(-x**2-y**2) on the square (0, 2) x (0, 2) for samples.
! The spline order is "nord" and the number of cells is
! "(ngrid-1)**2". There are "ndata" data values in the square.
integer :: i
integer, parameter :: ngrid=9, nord=4, ndegree=nord-1, &
nbkpt=ngrid+2*ndegree, ndata = 2000, nvalues=100
real(kind(1d0)), parameter :: zero=0d0, one=1d0, two=2d0
real(kind(1d0)), parameter :: TOLERANCE=1d-3
real(kind(1d0)), target :: spline_data (4, ndata), bkpt(nbkpt), &
coeff(ngrid+ndegree-1,ngrid+ndegree-1), delta, sizev, &
x(nvalues), y(nvalues), values(nvalues, nvalues)
real(kind(1d0)), pointer :: pointer_bkpt(:)
type (d_spline_knots) knotsx, knotsy
! Generate random (x,y) pairs and evaluate the
! example exponential function at these values.
spline_data(1:2,:)=two*rand(spline_data(1:2,:))
spline_data(3,:)=exp(-sum(spline_data(1:2,:)**2,dim=1))
spline_data(4,:)=one
! Define the knots for the tensor product data fitting problem.
delta = two/(ngrid-1)
bkpt(1:ndegree) = zero
bkpt(nbkpt-ndegree+1:nbkpt) = two
bkpt(nord:nbkpt-ndegree)=(/(i*delta,i=0,ngrid-1)/)
! Assign the degree of the polynomial and the knots.
pointer_bkpt => bkpt
knotsx=d_spline_knots(ndegree, pointer_bkpt)
knotsy=knotsx
! Fit the data and obtain the coefficients.
coeff = surface_fitting(spline_data, knotsx, knotsy)
! Evaluate the residual = spline - function
! at a grid of points inside the square.
delta=two/(nvalues+1)
x=(/(i*delta,i=1,nvalues)/); y=x
values=exp(-spread(x**2,1,nvalues)-spread(y**2,2,nvalues))
values=surface_values((/0,0/), x, y, knotsx, knotsy, coeff)-&
values
! Compute the R.M.S. error:
sizev=norm(pack(values, (values == values)))/nvalues
if (sizev <= TOLERANCE) then
write(*,*) 'Example 1 for SURFACE_FITTING is correct.'
end if
end
Example 1 for SURFACE_FITTING is correct.
Example 2: Parametric Representation of a Sphere
From Struik (1961), the parametric representation of points (x,y,z) on the surface of a sphere of radius a > 0 is expressed in terms of spherical coordinates,
The parameters are radians of latitude (u) and longitude (ν). The example program fits the same ndata random pairs of latitude and longitude in each coordinate. We have covered the sphere twice by allowing:
for latitude. We solve three data fitting problems, one for each coordinate function. Periodic constraints on the value of the spline are used for both u and ν. We could reduce the computational effort by fitting a spline function in one variable for the z coordinate. To illustrate the representation of more general surfaces than spheres, we did not do this. When the surface is evaluated we compute latitude, moving from the South Pole to the North Pole,
Our surface will approximately satisfy the equality
These residuals are checked at a rectangular mesh of latitude and longitude pairs. To illustrate the use of some options, we have reset the three regularization parameters to the value zero, the least-squares system tolerance to a smaller value than the default, and obtained the residuals for each parametric coordinate function at the data points.
USE surface_fitting_int
USE rand_int
USE norm_int
USE Numerical_Libraries
implicit none
! This is Example 2 for SURFACE_FITTING, tensor product
! B-splines approximation. Fit x, y, z parametric functions
! for points on the surface of a sphere of radius “A”.
! Random values of latitude and longitude are used to generate
! data. The functions are evaluated at a rectangular grid
! in latitude and longitude and checked to lie on the surface
! of the sphere.
integer :: i, j
integer, parameter :: ngrid=6, nord=6, ndegree=nord-1, &
nbkpt=ngrid+2*ndegree, ndata =1000, nvalues=50, NOPT=5
real(kind(1d0)), parameter :: zero=0d0, one=1d0, two=2d0
real(kind(1d0)), parameter :: TOLERANCE=1d-2
real(kind(1d0)), target :: spline_data (4, ndata, 3), bkpt(nbkpt), &
coeff(ngrid+ndegree-1,ngrid+ndegree-1, 3), delta, sizev, &
pi, A, x(nvalues), y(nvalues), values(nvalues, nvalues), &
data(4,ndata)
real(kind(1d0)), pointer :: pointer_bkpt(:)
type (d_spline_knots) knotsx, knotsy
type (d_options) OPTIONS(NOPT)
! Get the constant "pi" and a random radius, > 1.
pi = DCONST("pi"); A=one+rand(A)
! Generate random (latitude, longitude) pairs and evaluate the
! surface parameters at these points.
spline_data(1:2,:,1)=pi*(two*rand(spline_data(1:2,:,1))-one)
spline_data(1:2,:,2)=spline_data(1:2,:,1)
spline_data(1:2,:,3)=spline_data(1:2,:,1)
! Evaluate x, y, z parametric points.
spline_data(3,:,1)=A*cos(spline_data(1,:,1))*cos(spline_data(2,:,1))
spline_data(3,:,2)=A*cos(spline_data(1,:,2))*sin(spline_data(2,:,2))
spline_data(3,:,3)=A*sin(spline_data(1,:,3))
! The values are equally uncertain.
spline_data(4,:,:)=one
! Define the knots for the tensor product data fitting problem.
delta = two*pi/(ngrid-1)
bkpt(1:ndegree) = -pi
bkpt(nbkpt-ndegree+1:nbkpt) = pi
bkpt(nord:nbkpt-ndegree)=(/(-pi+i*delta,i=0,ngrid-1)/)
! Assign the degree of the polynomial and the knots.
pointer_bkpt => bkpt
knotsx=d_spline_knots(ndegree, pointer_bkpt)
knotsy=knotsx
! Fit a data surface for each coordinate.
! Set default regularization parameters to zero and compute
! residuals of the individual points. These are returned
! in DATA(4,:).
do j=1,3
data=spline_data(:,:,j)
OPTIONS(1)=d_options(surface_fitting_thinness,zero)
OPTIONS(2)=d_options(surface_fitting_flatness,zero)
OPTIONS(3)=d_options(surface_fitting_smallness,zero)
OPTIONS(4)=d_options(surface_fitting_tol_least,1d-5)
OPTIONS(5)=surface_fitting_residuals
coeff(:,:,j) = surface_fitting(data, knotsx, knotsy,&
IOPT=OPTIONS)
end do
! Evaluate the function at a grid of points inside the rectangle of
! latitude and longitude covering the sphere just once. Add the
! sum of squares. They should equal "A**2" but will not due to
! truncation and rounding errors.
delta=pi/(nvalues+1)
x=(/(-pi/two+i*delta,i=1,nvalues)/); y=two*x
values=zero
do j=1,3
values=values+&
surface_values((/0,0/), x, y, knotsx, knotsy, coeff(:,:,j))**2
end do
values=values-A**2
! Compute the R.M.S. error:
sizev=norm(pack(values, (values == values)))/nvalues
if (sizev <= TOLERANCE) then
write(*,*) "Example 2 for SURFACE_FITTING is correct."
end if
end
Example 2 for SURFACE_FITTING is correct.
Example 3: Constraining Some Points using a Spline Surface
This example illustrates the use of discrete constraints to shape the surface. The data fitting problem of Example 1 is modified by requiring that the surface interpolate the value one at x = y = 0. The shape is constrained so first partial derivatives in both x and y are zero at x = y = 0. These constraints mimic some properties of the function g(x,y). The size of the residuals at a grid of points and the residuals of the constraints are checked.
USE surface_fitting_int
USE rand_int
USE norm_int
implicit none
! This is Example 3 for SURFACE_FITTING, tensor product
! B-splines approximation, f(x,y). Use the function
! exp(-x**2-y**2) on the square (0, 2) x (0, 2) for samples.
! The spline order is "nord" and the number of cells is
! "(ngrid-1)**2". There are "ndata" data values in the square.
! Constraints are put on the surface at (0,0). Namely
! f(0,0) = 1, f_x(0,0) = 0, f_y(0,0) = 0.
integer :: i
integer, parameter :: ngrid=9, nord=4, ndegree=nord-1, &
nbkpt=ngrid+2*ndegree, ndata = 2000, nvalues=100, NC = 3
real(kind(1d0)), parameter :: zero=0d0, one=1d0, two=2d0
real(kind(1d0)), parameter :: TOLERANCE=1d-3
real(kind(1d0)), target :: spline_data (4, ndata), bkpt(nbkpt), &
coeff(ngrid+ndegree-1,ngrid+ndegree-1), delta, sizev, &
x(nvalues), y(nvalues), values(nvalues, nvalues), &
f_00, f_x00, f_y00
real(kind(1d0)), pointer :: pointer_bkpt(:)
type (d_spline_knots) knotsx, knotsy
type (d_surface_constraints) C(NC)
LOGICAL PASS
! Generate random (x,y) pairs and evaluate the
! example exponential function at these values.
spline_data(1:2,:)=two*rand(spline_data(1:2,:))
spline_data(3,:)=exp(-sum(spline_data(1:2,:)**2,dim=1))
spline_data(4,:)=one
! Define the knots for the tensor product data fitting problem.
delta = two/(ngrid-1)
bkpt(1:ndegree) = zero
bkpt(nbkpt-ndegree+1:nbkpt) = two
bkpt(nord:nbkpt-ndegree)=(/(i*delta,i=0,ngrid-1)/)
! Assign the degree of the polynomial and the knots.
pointer_bkpt => bkpt
knotsx=d_spline_knots(ndegree, pointer_bkpt)
knotsy=knotsx
! Define the constraints for the fitted surface.
C(1)=surface_constraints(point=(/zero,zero/),type='==',value=one)
C(2)=surface_constraints(derivative=(/1,0/),&
point=(/zero,zero/),type='==',value=zero)
C(3)=surface_constraints(derivative=(/0,1/),&
point=(/zero,zero/),type='==',value=zero)
! Fit the data and obtain the coefficients.
coeff = surface_fitting(spline_data, knotsx, knotsy,&
CONSTRAINTS=C)
! Evaluate the residual = spline - function
! at a grid of points inside the square.
delta=two/(nvalues+1)
x=(/(i*delta,i=1,nvalues)/); y=x
values=exp(-spread(x**2,1,nvalues)-spread(y**2,2,nvalues))
values=surface_values((/0,0/), x, y, knotsx, knotsy, coeff)-&
values
f_00 = surface_values((/0,0/), zero, zero, knotsx, knotsy, coeff)
f_x00= surface_values((/1,0/), zero, zero, knotsx, knotsy, coeff)
f_y00= surface_values((/0,1/), zero, zero, knotsx, knotsy, coeff)
! Compute the R.M.S. error:
sizev=norm(pack(values, (values == values)))/nvalues
PASS = sizev <= TOLERANCE
PASS = abs (f_00 - one) <= sqrt(epsilon(one)) .and. PASS
PASS = f_x00 <= sqrt(epsilon(one)) .and. PASS
PASS = f_y00 <= sqrt(epsilon(one)) .and. PASS
if (PASS) then
write(*,*) 'Example 3 for SURFACE_FITTING is correct.'
end if
end
Example 3 for SURFACE_FITTING is correct.
Example 4: Constraining a Spline Surface to be non-Negative
The review of interpolating methods by Franke (1982) uses a test data set originally due to James Ferguson. We use this data set of 25 points, with unit uncertainty for each dependent variable. Our algorithm does not interpolate the data values but approximately fits them in the least-squares sense. We reset the regularization parameter values of flatness and thinness, Hanson (1995). Then the surface is fit to the data and evaluated at a grid of points. Although the surface appears smooth and fits the data, the values are negative near one corner. Our scenario for the application assumes that the surface be non-negative at all points of the rectangle containing the independent variable data pairs. Our algorithm for constraining the surface is simple but effective in this case. The data fitting is repeated one more time but with positive constraints at the grid of points where it was previously negative.
USE surface_fitting_int
USE rand_int
USE surface_fitting_int
USE rand_int
USE norm_int
implicit none
! This is Example 4 for SURFACE_FITTING, tensor product
! B-splines approximation, f(x,y). Use the data set from
! Franke, due to Ferguson. Without constraints the function
! becomes negative in a corner. Constrain the surface
! at a grid of values so it is non-negative.
integer :: i, j, q
integer, parameter :: ngrid=9, nord=4, ndegree=nord-1, &
nbkpt=ngrid+2*ndegree, ndata = 25, nvalues=50
real(kind(1d0)), parameter :: zero=0d0, one=1d0
real(kind(1d0)), parameter :: TOLERANCE=1d-3
real(kind(1d0)), target :: spline_data (4, ndata), bkptx(nbkpt), &
bkpty(nbkpt),coeff(ngrid+ndegree-1,ngrid+ndegree-1), &
x(nvalues), y(nvalues), values(nvalues, nvalues), &
delta
real(kind(1d0)), pointer :: pointer_bkpt(:)
type (d_spline_knots) knotsx, knotsy
type (d_surface_constraints), allocatable :: C(:)
real(kind(1e0)) :: data (3*ndata) = & ! This is Ferguson's data:
(/2.0 , 15.0 , 2.5 , 2.49 , 7.647, 3.2,&
2.981 , 0.291, 3.4 , 3.471, -7.062, 3.5,&
3.961 , -14.418, 3.5 , 7.45 , 12.003, 2.5,&
7.35 , 6.012, 3.5 , 7.251, 0.018, 3.0,&
7.151 , -5.973, 2.0 , 7.051, -11.967, 2.5,&
10.901, 9.015, 2.0 , 10.751, 4.536, 1.925,&
10.602, 0.06 , 1.85, 10.453, -4.419, 1.576,&
10.304, -8.895, 1.7 , 14.055, 10.509, 1.5,&
14.194, 6.783, 1.3 , 14.331, 3.054, 1.7,&
14.469, -0.672, 2.1 , 14.607, -4.398, 1.75,&
15.0 , 12.0 , 0.5 , 15.729, 8.067, 0.5,&
16.457, 4.134, 0.7 , 17.185, 0.198, 1.1,&
17.914, -3.735, 1.7/)
spline_data(1:3,:)=reshape(data,(/3,ndata/)); spline_data(4,:)=one
! Define the knots for the tensor product data fitting problem.
! Use the data limits to the knot sequences.
bkptx(1:ndegree) = minval(spline_data(1,:))
bkptx(nbkpt-ndegree+1:nbkpt) = maxval(spline_data(1,:))
delta=(bkptx(nbkpt)-bkptx(ndegree))/(ngrid-1)
bkptx(nord:nbkpt-ndegree)=(/(bkptx(1)+i*delta,i=0,ngrid-1)/)
! Assign the degree of the polynomial and the knots for x.
pointer_bkpt => bkptx
knotsx=d_spline_knots(ndegree, pointer_bkpt)
bkpty(1:ndegree) = minval(spline_data(2,:))
bkpty(nbkpt-ndegree+1:nbkpt) = maxval(spline_data(2,:))
delta=(bkpty(nbkpt)-bkpty(ndegree))/(ngrid-1)
bkpty(nord:nbkpt-ndegree)=(/(bkpty(1)+i*delta,i=0,ngrid-1)/)
! Assign the degree of the polynomial and the knots for y.
pointer_bkpt => bkpty
knotsy=d_spline_knots(ndegree, pointer_bkpt)
! Fit the data and obtain the coefficients.
coeff = surface_fitting(spline_data, knotsx, knotsy)
delta=(bkptx(nbkpt)-bkptx(1))/(nvalues+1)
x=(/(bkptx(1)+i*delta,i=1,nvalues)/)
delta=(bkpty(nbkpt)-bkpty(1))/(nvalues+1)
y=(/(bkpty(1)+i*delta,i=1,nvalues)/)
! Evaluate the function at a rectangular grid.
! Use non-positive values to a constraint.
values=surface_values((/0,0/), x, y, knotsx, knotsy, coeff)
! Count the number of values <= zero. Then constrain the spline
! so that it is >= TOLERANCE at those points where it was <= zero.
q=count(values <= zero)
allocate (C(q))
DO I=1,nvalues
DO J=1,nvalues
IF(values(I,J) <= zero) THEN
C(q)=surface_constraints(point=(/x(i),y(j)/), type='>=',&
value=TOLERANCE)
q=q-1
END IF
END DO
END DO
! Fit the data with constraints and obtain the coefficients.
coeff = surface_fitting(spline_data, knotsx, knotsy,&
CONSTRAINTS=C)
deallocate(C)
! Evaluate the surface at a grid and check, once again, for
! non-positive values. All values should now be positive.
values=surface_values((/0,0/), x, y, knotsx, knotsy, coeff)
if (count(values <= zero) == 0) then
write(*,*) 'Example 4 for SURFACE_FITTING is correct.'
end if
end
Example 4 for SURFACE_FITTING is correct.
