Bivariate smoothing spline for scalar data defined over an elliptic domain

Hello,

I am interested in fitting a 3D smoothing spline to scalar, gridded, cartesian data that is non-zero only in a non-rectangular volume (e.g., cylinder, sphere, other volumes) embedded within a cartesian grid. For elements (voxels) outside the volume, the scalar field is zero. I would like to fit a smoothing quintic spline because I need to take multiple spatial derivatives of the scalar field. For this reason, I have been working with the function spaps.

The issue is that spaps fits to the zeros outside the non-rectangular volume, severely affecting the quality of the fit, especially for high tolerance (high smoothing). In 1D this can be addressed by specifying a vector of weights such that zeros outside the interval of interest are ignored (weights are set to zero for those elements and one for elements in the interval of interest). However, when extending this idea to multivariate splines, I noticed that the weights can only be specified as a vector for each spatial dimension, as opposed to a 2D or 3D array of weights, so it becomes difficult to zero the weights outside of the non-rectangular volume of interest.

For the purposes of illustrating the issue, I have inserted code in which I define a hypothetical flow field through a duct with elliptic cross section and have added a bit of Gaussian noise. The velocity component parallel to the duct axis (scalar quantity) is a noisy elliptic paraboloid.

%Bivariate tensor product smoothing spline fit to elliptic paraboloid 2D
%flow field representing laminar flow through a duct with elliptic
%cross-section

%Define vectors for each spatial coordinate
xVec = -1:0.01:1; 
yVec = -0.5:0.01:0.5; 

%Define mesh using meshgrid. spaps documentation recommends ndgrid, but spaps
%seems to work with meshgrid if y is passed as the first index and x as the second index.
[xMesh,yMesh] = meshgrid(xVec,yVec);

%Define 2D flow field for elliptic parabolic flow (magnitude of out of
%plane, z-velocity, as a function of position x and y)
a = 0.5; %Major axis along x direction
b = 0.25; %Minor axis along y direction
maxVel = 1; %Maximum velocity
zVel = -1*((xMesh./a).^2 + (yMesh./b).^2) + maxVel; %Velocity field
ductMask = zVel > 0; %Binary mask 1 in duct 0 outside
zVel(zVel <= 0) = 0; %Zero all negative velocities (outside duct)

%Scaled image plot of z-velocity component in an xy plane
figure
imagesc([min(xVec) max(xVec)],[min(yVec) max(yVec)],zVel)
xlabel("x-position (m)")
ylabel("y-position (m)")
clbrObj = colorbar;
clbrObj.Label.String = 'z-velocity (m/s)';
axis image

%Add gaussian noise to the z-velocity component
rng(39);
noiseFrac = 0.01;
zVelNoisy = zVel + noiseFrac*max(zVel,[],"all")*randn(size(zVel));
zVelNoisy(~ductMask) = 0; %Zero velocities outside duct

I proceed to fit a bivariate quintic spline with a large tolerance (high smoothing) to this scalar field and plot the velocity component along a cutline through the minor axis of the elliptic cross section and one parallel to it but off-center. It is evident that the zero velocity component values outside the duct are affecting the fit inside the duct.

%----Uniform weights----%

%Define tolerances and weights for each spatial coordinate.
%Same tolerance in both directions and weights in both
%directions are 1 everywhere (default)
xTol = 3;
yTol = 3; 
xWeights = ones(size(xVec));
yWeights = ones(size(yVec));

%Fit bivariate quintic smoothing spline
zVel_pp = fn2fm(spaps({yVec,xVec},zVelNoisy,{yTol xTol},...
    {yWeights,xWeights},3),"pp");
zVelSmooth = fnval(zVel_pp,{yVec,xVec});

%Plot z-velocity component y-direction cutline through the minor axis
figure
hold on
plot(yVec,zVel(:,101),"LineWidth",2)
plot(yVec,zVelSmooth(:,101),"LineWidth",2)
plot(yVec,zVelNoisy(:,101),"ko")
ylim([-0.2 1.2])
legend("Exact","Spline","Data")
xlabel("y-position (m)")
ylabel("z-velocity (m)")

%Plot z-velocity components for a y-direction cutline
%parallel to the minor axis
figure
hold on
plot(yVec,zVel(:,63),"LineWidth",2)
plot(yVec,zVelSmooth(:,63),"LineWidth",2)
plot(yVec,zVelNoisy(:,63),"ko")
ylim([-0.1 0.5])
legend("Exact","Spline","Data")
xlabel("y-position (m)")
ylabel("z-velocity (m)")

I then pass weight vectors that set the fitting weights to zero outside a rectangle bounding the ellipse and to one within the rectangle. This greatly improves the fit for the cutline through the minor axis, but not through the parallel cutline because some zero velocity values outside the ellipse are still weighted for this cutline.

%----Non-uniform weights----%

%Define tolerances and weights for each spatial coordinate.
%Specify weight vectors which are zero outside rectangle containing duct
xWeights = zeros(size(xVec));
yWeights = zeros(size(yVec));
[row,col] = find(ductMask);
xWeights(min(col):max(col)) = 1;
yWeights(min(row):max(row)) = 1;

%Fit bivariate quintic smoothing spline with non-uniform weights
zVelWeighted_pp = fn2fm(spaps({yVec,xVec},zVelNoisy,{yTol xTol},...
    {yWeights,xWeights},3),"pp");
zVelWeightedSmooth = fnval(zVelWeighted_pp,{yVec,xVec});

%Plot z-velocity component y-direction cutline through the minor axis.
%Unweighted (weight = 0) data points are not shown.
figure
hold on
plot(yVec,zVel(:,101),"LineWidth",2)
plot(yVec,zVelWeightedSmooth(:,101),"LineWidth",2)
plot(yVec(min(row):max(row)),zVelNoisy(min(row):max(row),101),"ko")
ylim([-0.2 1.2])
legend("Exact","Spline","Data")
xlabel("y-position (m)")
ylabel("z-velocity (m)")

%Plot z-velocity components for a y-direction cutline
%parallel to the minor axis
%Unweighted (weight = 0) data points are not shown.
figure
hold on
plot(yVec,zVel(:,63),"LineWidth",2)
plot(yVec,zVelWeightedSmooth(:,63),"LineWidth",2)
plot(yVec(min(row):max(row)),zVelNoisy((min(row):max(row)),63),"ko")
ylim([-0.1 0.5])
legend("Exact","Spline","Data")
xlabel("y-position (m)")
ylabel("z-velocity (m)")

Can a matrix of weights be passed to spaps to unweight elements outside the ellipse? If not, is there another way to ensure that a smoothing spline is only fit within the elliptical domain (or other non-rectangular domain)?

Answer 1

[Abhinaya Kennedy]

Hi Julian,

You are looking to fit a 3D smoothing spline to data within a non-rectangular volume while ignoring or minimizing the influence of the surrounding zero-valued data, and if “spaps” can be used with a 3D weight matrix.

One approach is to create a custom weight matrix that reflects the shape of your domain (e.g. Ellipse). This matrix would then need to be flattened into a vector to serve as weights within “spaps”. Alternatively, you could parameterize your domain, and fit the spline within this more regular parametric space.

Another option is to divide the domain into rectangular subdomains, fit splines to each of these areas, and ensure that they are continuous at the boundaries. You might consider other fitting techniques like radial basis functions (RBFs) or kriging, which are better suited to irregular domains.

You could also preprocess your data to substitute zero values outside the domain with NaNs or a similar placeholder, allowing the use of a fitting method that can handle missing data. Lastly, developing a custom fitting function that creates a weight matrix specifically for your domain and iteratively fits a spline could be effective.

You can find additional information on “spaps” here: https://www.mathworks.com/help/curvefit/spaps.html?searchHighlight=spaps&s_tid=srchtitle_support_results_1_spaps

Hope this helps!