(Using Blender 4.2.5)
Results and resources
Figure 1: Displacement interpolation on a cylindrical grid (Top) and on a uniform Cartesian grid (Bottom). Ico Spheres are displayed at the control points, except in the center.
Figure 2: (Top) Displacement interpolation on a butterfly grid. (Bottom) Butterfly grid made with a XY Mirror modifier, followed by a Catmull-Clark Subdivision modifier.
Resources:
(Blender 4.2.5+)
Overview
To achieve a smooth surface mesh from scattered and highly varying control points, the proposed approach proceeds in two steps. Firstly Catmull Rom type splines are built in azimuthal and in radial directions to define a smooth surface, but over a distorted grid. Then the vertical displacement is interpolated from this surface to a regular grid.
For the demonstration, cylindrical and Cartesian grids are illustrated Fig.(1). Cartesian grid removes high grid density around the origin, but it requires some additional treatment along the border. Taking advantage that the interpolation can be computed over any user defined grid, a mesh with optimized topology can be provided instead of building the grid using nodes (see Fig.(2)).
GeometryNodes Modifier
Nomenclature
Let $(X,Y,Z)$ be the Cartesian coordinates, and $(\theta,r,Z)$ be the cylindrical ones. They are related by the following expressions: $$ \begin{array}{lcl} \left\{ \begin{array}{rcl} X & = & r \cos{(\theta)} \\ Y & = & r \sin{(\theta)} \end{array} \right. & \mbox{} & \left\{ \begin{array}{rcl} r & = & \sqrt{X^2+Y^2} \\ \theta & = & \arctan{(Y/X)} \end{array} \right. \end{array} $$
Let $R$ be the largest radius. The interpolation surface is restricted to $(\theta,r) \in \left[ 0 : \pi \right] \times \left[ -R : R \right]$ to enforce the continuity and the smoothness at $r=0$, instead of the classical interval $\left[ -\pi : \pi \right] \times \left[ 0 : R \right]$. Both intervals are spanning the same Cartesian space because: $$ \left\{ \begin{array}{rclcl} X & = & r \cos{(\theta)} & = & -r \cos{(\theta +\pi)} \\ Y & = & r \sin{(\theta)} & = & -r \sin{(\theta +\pi)} \end{array} \right. $$ So $(\theta,r) \in \left[ -\pi : 0 \right] \times \left[ 0 : R \right]$ is spanning the same Cartesian space as $(\theta,r) \in \left[ 0 : \pi \right] \times \left[ -R : 0 \right]$.
Main graph
Figure 3: Graph building and displaying the interpolated surface.
After the 5 meshes defining the control points are joined in a single geometry, the process is split in 6 group nodes:
Transfer: Group making azimuthal profiles of control points, from user defined radial profiles.Center: Group averaging the control points at the center, to smooth discontinuities.Azimuthal: Group making splines in azimuthal direction.Radial: Group making splines in radial direction.Mapping: Group building the interpolation surface of the vertical displacement.Surface: Group making a sampling grid and interpolating the displacement.
Transfer group node
Figure 4: Graph making azimuthal profiles of control points, from user defined radial profiles.
Dark blue nodes: Extracting the control points
Control points with negative $Z$ coordinate are deleted, to model only the upper surface. It is assumed that each user defined radial profile makes a mesh island, that the number of remaining vertices per mesh island is the same, and that mesh islands are evenly distributed in azimuthal direction (starting at $\theta=0$).
Dark green nodes: Making azimuthal mesh lines
The starting geometry is a Grid such that its $X \equiv \theta$ coordinate is the angular position along an azimuthal mesh line, while its $Y \equiv r$ coordinate is the radial position. Consequently, the Size X parameter is set to $2\pi$, for $\theta$ to be in $\left[ -\pi : \pi \right]$, and the Vertices X parameter is set to the count of mesh islands plus 1, to duplicate the first/last user defined radial profile. The Vertices Y parameter is set to the number of control points per user defined radial profile, recovered by dividing the total number of vertices by the count of mesh islands, and the Size Y parameter is set to 0 because the $r \equiv Y$ coordinate is offset afterwards. Eventually, edges Not aligned with the vector $(1,0,0)$ (i.e. in radial direction) are deleted to keep only edges aligned with $\theta \equiv X$ direction (i.e. in azimuthal direction).
Dark red nodes: Duplicating the first/last control point per azimuthal mesh line
Grid points with Index greater than or equal to the Point Count are copies of the first radial profile. So the index of the point to copy is computed with a Truncated Modulo math node.
Dark pink nodes: Transferring control points coordinates
The Position of the point to copy is recovered through a Sample Index node. Its $\theta$ coordinate is ignored, and its $(r,Z)$ coordinates are combined to offset the grid points.
Center group node
Figure 5: Graph averaging the control points at the center.
The input Geometry is a set of mesh lines, defining the vertical displacement along profiles at constant radius. The control points coordinates are such that $\theta \equiv X \in \left[ -\pi : \pi \right]$ and $r \equiv Y \in \left[ 0 : R \right]$. The $Z$ displacement along the azimuthal profile with $r=0$ is not constant, leading to discontinuities at the origin. Consequently the mean displacement along this profile is substituted to the user defined values.
It is assumed that the Mesh Island with Island Index equal to 0 is the azimuthal profile with $r=0$. This mesh island is selected taking advantage that the type casting of the integer 0 into a boolean value is False, while the type casting of positive integers into a boolean value is True. This boolean value is reversed using a Not boolean math node before it is input into an Attribute Statistic node computing the Mean of $r \equiv Z$, and also into the Set Position node overwriting the previous position.
Azimuthal group node
Figure 6.1: Graph making splines in azimuthal direction.
The input Geometry is a set of mesh lines, defining the vertical displacement along profiles at constant radius. The control points coordinates are such that $\theta \equiv X \in \left[ -\pi : \pi \right]$ and $r \equiv Y \in \left[ 0 : R \right]$.
Dark pink nodes: Extending mesh lines in azimuthal direction
To enforce the continuity and the smoothness of the interpolating spline at $\theta \equiv X =\pm\pi$, 3 copies of each profile are merged together, to mimic periodic conditions for the central copy. The leftmost copy (i.e. with Duplicate Index of 0) is translated in $\theta \equiv X$ direction by $-\pi$, the central copy (i.e. with Duplicate Index of 1) is translated by $\pi$, and the rightmost copy (i.e. with Duplicate Index of 2) is translated by $3\pi$. So $\theta \equiv X \in \left[ -2\pi : 4\pi \right]$ after this step, with $\left[ 0 : 2\pi \right]$ the interval to accurately model.
Dark red nodes: Making splines in azimuthal direction
The resulting mesh lines are converted in Catmull Rom type splines using a Mesh to Curve/Set Spline Type combo. The Resolution is adjusted to achieve a smooth evolution (see Fig.(6.2-Top)).
Figure 6.2: (Top) Illustration of Catmull Rom type spline and control points, displayed in the $(\theta \equiv X,Z)$ plane. (Bottom) 3D view of the set of azimuthal splines.
Radial group node
Figure 7.1: Graph making splines in radial direction.
The input Geometry is a set of splines, collected as a single curve, modelling the smooth vertical displacement along profiles at constant radius (see Fig.(6.2)). The control points coordinates are such that $\theta \equiv X \in \left[ -2\pi : 4\pi \right]$ and $r \equiv Y \in \left[ 0 : R \right]$.
Dark green nodes: Making sampling mesh lines
The starting geometry is a Grid such that its $X \equiv f$ coordinate is the Factor along an azimuthal spline, while its $Y \equiv i$ coordinate is the index of such a spline (starting at 0). Consequently, the Size X parameter is set to 1, for the Factor $f$ to be in $\left[ 0 : 1 \right]$, and the Vertices X parameter is set to 1 plus 6 times the number of edges spanning $\pi$. The factor 6 is to discretize the $\mathbf{6}\pi$ wide interval in $\theta \equiv X$ direction of the input Geometry. The Vertices Y parameter is set to the number of azimuthal splines, recovered through a Domain Size node, and the Size Y parameter is set to this count less 1 for the spacing in $Y \equiv i$ direction to be 1. The resulting grid is offset with a Transform Geometry node by half its size such that its lower left corner is at $(f,i) \equiv (X,Y) = (0,0)$. Eventually, edges aligned with the vector $(1,0,0)$ (i.e. at constant spline index $i$) are deleted to keep only edges aligned with $Y$ direction (i.e. at constant $f$ value).
Dark pink nodes: Sampling azimuthal splines
$(\theta,r,Z)$ coordinates are interpolated along azimuthal splines from the $(f,i) \equiv (X,Y)$ doublet using a Sample Curve node set in Factor mode. Because of the slope in Z direction, a uniform distribution of curvilinear abscissa along an azimuthal spline is not matching a uniform distribution in $\theta$ direction. Consequently, the achieved mesh lines are not at constant $\theta$ value.
Dark blue nodes: Extending the modelling space in radial direction
The half space with $r \equiv Y \le 0$ (see the Nomenclature subsection) is filled by offsetting the $\theta \equiv X$ coordinate by $-\pi$ and by multiplying the $r \equiv Y$ coordinate by $-1$ using a Transform Geometry node. Both sets of mesh lines are joined then merged to make continuous mesh lines extending in radial (i.e. $Y$) direction from $-R$ to $+R$.
Dark red nodes: Making splines in radial direction
The resulting mesh lines are converted in Catmull Rom type splines using a Mesh to Curve/Set Spline Type combo. The Resolution is adjusted to achieve a smooth evolution (see Fig.(7.2-Top)).
Dark grey nodes: Restricting the modelling space in azimuthal direction
Curves not fully spanning the $r \equiv Y \ \left[ -R : R \right]$ interval are deleted. By construction, these are such that for $r = 0$, $\theta \lt -\pi$ or $\theta \gt 2\pi$. The value of $\theta(r=0)$ is approximated as $\theta(f=0.5)$. This $\theta \equiv X$ value is computed by a Sample Curve node for which the Curve Index is provided by a Curve of Point node. The achieved set of splines is illustrated Fig.(7.2-Bottom).
Figure 7.2: (Top) Catmull Rom type spline and control points close to $\theta \equiv X =0$, displayed in the $(r \equiv Y,Z)$ plane. (Bottom) 3D view of the set of radial splines, with an ellipse extruded along each spline, for "Count Theta" set to 8.
Mapping group node
Figure 8.1: Graph building the interpolation surface of the displacement in cylindrical coordinates.
The input Geometry is a set of splines, collected as a single curve, modelling the smooth vertical displacement along radial profiles (see Fig.(7.2)). The control points coordinates are such that $\theta \equiv X \in \left[ -\pi : 2\pi \right]$ and $r \equiv Y \in \left[ -R : R \right]$.
Dark green nodes: Making the sampling grid
The starting geometry is a Grid such that its $Y \equiv f$ coordinate is the Factor along a radial spline, while its $X \equiv i$ coordinate is the index of such a spline (starting at 0). Consequently, the Size Y parameter is set to 1, for the Factor $f$ to be in $\left[ 0 : 1 \right]$, and the Vertices Y parameter is set to 1 plus 2 times the number of edges spanning $R$. The factor 2 is to discretize the $\mathbf{2}R$ wide interval in $r \equiv Y$ direction of the input Geometry. The Vertices X parameter is set to the number of radial splines, recovered through a Domain Size node, and the Size X parameter is set to this count less 1 for the spacing in $X \equiv i$ direction to be 1. The resulting grid is offset with a Transform Geometry node by half its size such that its lower left corner is at $(i,f) \equiv (X,Y) = (0,0)$.
Dark pink nodes: Sampling radial splines
$(\theta,r,Z)$ coordinates are interpolated along radial splines from the $(i,f) \equiv (X,Y)$ doublet using a Sample Curve node set in Factor mode. Because of the slope in Z direction, a uniform distribution of curvilinear abscissa along a radial spline is not matching a uniform distribution in $r$ direction. Consequently, the achieved mesh lines are not at constant $r$ value (see Fig.(8.2)).
Figure 8.2: (Top) 3D view of the interpolation surface. (Bottom) Top view of the interpolation surface, for "Count Theta" and "Count R" set to 8.
Surface group node
Figure 9: Graph making a sampling grid and interpolating the displacement.
The input Geometry is the interpolation surface with $(\theta,r) \equiv (X,Y) \in \left[ -\pi : 2\pi \right] \times \left[ -R : R \right]$ (see Fig.(8.2)). The largest radius $R$ is recovered as the highest value in $Y$ direction through a Bounding Box node.
Dark blue nodes: Making the Cartesian sampling grid
This grid size is made 10% larger than the original data set diameter, thus the factor 2.1 multiplying $R$. The Grid node is returning a mesh in the $XY$ plane with $Z=0$. So the conversion from $(X,Y)$ to $(\theta,r)$, stored as $(X,Y)$, is straightforward (see the Nomenclature subsection).
Dark red nodes: Making the cylindrical sampling grid
Setting the Size X parameter to $2\pi$ and the Size Y parameter to $R$, the Grid node is returning a mesh with $(\theta,r) \equiv (X,Y) \in \left[ -\pi : \pi \right] \times \left[ -\frac{1}{2}R : \frac{1}{2}R \right]$. Consequently, it is offset by $\frac{1}{2}R$ in $Y$ direction using a Transform Geometry node such that $r \equiv Y \in \left[ 0 : R \right]$.
Dark green nodes: Shaping the coordinates to interpolate
The sampling half-space $(\theta,r) \equiv (X,Y) \in \left[ -\pi : 0 \right[ \times \left[ 0 : R \right]$ is offset by $\pi$ in $\theta \equiv X$ direction (see the Add math node), and inverted in $r \equiv Y$ direction (see the Multiply math node), to match the sampling half-space $(\theta,r) \equiv (X,Y) \in \left[ 0 : \pi \right[ \times \left[ -R : 0 \right]$ (see the Nomenclature subsection).
Dark pink nodes: Interpolating
The transfer of the displacement in $Z$ direction from the input interpolation surface to the output sampling mesh is based on a Sample Nearest Surface node. To restrict the interpolation in the $XY$ plane with $Z=0$, the $Z$ coordinate of the interpolation surface is nullified (see the Multiply vector math node) before it is input as Source Mesh. Consequently, this coordinate must be preserved through a Capture Attribute node beforehand to be transferred.
Dark grey nodes: Computing Cartesian coordinates
The conversion from $(X,Y,Z) \equiv (\theta,r,Z)$ to $(X,Y,Z)$ is straightforward (see the Nomenclature subsection).
