寻求任意区域内三次样条曲面插值方案（附椭圆抛物面示例）

Hey there! I totally get the frustration of only finding cubic spline surface interpolation examples limited to the [0,1]×[0,1] unit square when you need to work with arbitrary domains—whether that's a grid of 16 points, a quadrilateral with 4 corner points and full derivative info, or your elliptic paraboloid test case. Let’s break down how to solve this:

1. Core Trick: Map Your Arbitrary Domain to the Unit Square #

Nearly all standard cubic spline implementations are built for the unit square because it simplifies the math, but you can easily adapt them to any rectangular (or even quadrilateral) domain using affine coordinate mapping. Here’s how:

Suppose your target domain is ( x \in [x_{\text{min}}, x_{\text{max}}] ) and ( y \in [y_{\text{min}}, y_{\text{max}}] ). Map each (x,y) point to a (u,v) point in [0,1]×[0,1] with:

u = (x - x_min) / (x_max - x_min)
v = (y - y_min) / (y_max - y_min)

Once you’ve interpolated values in the (u,v) space, you don’t need to map back—just use the (u,v) to (x,y) relationship to target any point in your original domain.

2. Handling Different Point Scenarios #

Scenario 1: Regular Grid Points (e.g., 3x3, 4x4/16 points) #

For structured grids like your elliptic paraboloid test case (( x \in {100,101,102} ), ( y \in {100,101,102} )), use a tensor-product cubic spline:

• First, map all your (x,y,f) points to (u,v,f) using the coordinate transform above.
• Use a standard tensor-product spline implementation (like scipy.interpolate.RectBivariateSpline in Python) on the (u,v,f) grid. This function handles cubic spline interpolation for 2D regular grids natively.
• When you need to interpolate at a new (x_target, y_target), convert it to (u_target, v_target) first, then query the spline.

Example snippet for your test case:

import numpy as np
from scipy.interpolate import RectBivariateSpline

# Your elliptic paraboloid points
x = np.array([100, 101, 102])
y = np.array([100, 101, 102])
f = np.array([
    [1805.56, 1827.67, 1850.88],  # Fill in your actual f(x,y) values
    [1827.67, 1849.78, 1872.99],
    [1850.88, 1872.99, 1896.20]
])

# Map to unit square
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()
u = (x - x_min) / (x_max - x_min)
v = (y - y_min) / (y_max - y_min)

# Create cubic spline
spline = RectBivariateSpline(u, v, f, kx=3, ky=3)  # kx/ky=3 specifies cubic order

# Interpolate at a new point (e.g., x=100.5, y=100.5)
x_target = 100.5
y_target = 100.5
u_target = (x_target - x_min) / (x_max - x_min)
v_target = (y_target - y_min) / (y_max - y_min)
f_target = spline(u_target, v_target)[0][0]

Scenario 2: 4 Corner Points with Full Derivative Info #

If you have 4 corner points (e.g., a quadrilateral) and know ( f, f_x, f_y, f_{xy} ) for each, use bicubic Hermite interpolation. This method is designed exactly for this case—no need for a full grid.

Steps:

1. Map your 4 corner (x,y) points to the unit square's corners (0,0), (1,0), (0,1), (1,1) using affine mapping.
2. For each mapped corner (u,v), convert the partial derivatives ( f_x, f_y, f_{xy} ) to the (u,v) coordinate system using the Jacobian of the transform:
   • ( f_u = f_x \cdot (x_{\text{max}} - x_{\text{min}}) )
   • ( f_v = f_y \cdot (y_{\text{max}} - y_{\text{min}}) )
   • ( f_{uv} = f_{xy} \cdot (x_{\text{max}} - x_{\text{min}}) \cdot (y_{\text{max}} - y_{\text{min}}) )
3. Construct the bicubic Hermite surface using the mapped (u,v) values and transformed derivatives. The surface is defined by 16 basis functions, each weighted by the corner values/derivatives.

You can implement this manually with the bicubic Hermite basis formulas, or use specialized interpolation libraries that support Hermite splines for 2D.

3. Notes for Non-Rectangular Arbitrary Domains #

If your domain isn’t a rectangle (e.g., an irregular polygon), you’ll need to use splines over triangulated domains (like thin-plate splines or radial basis functions) instead of tensor-product methods. But for most practical cases (like your rectangular elliptic paraboloid example), the affine mapping + standard spline approach works perfectly.

内容的提问来源于stack exchange，提问作者user9582211