In this section, we develop a general method for smoothing functions on surfaces, which are represented by triangular meshes. The proposed algorithm consists of several steps: optimization of the SF, extraction of the surface mesh information, a smoothing of the operator structuring, a smoothing process and loop parameter calculation, and finally an evaluation of the smoothing effectiveness based on the curvature change in the contour lines of the SF.
Figure 1 illustrates the process of these steps.
2.1. Previous Work
The SFM is one of the most widely used methods in gradient coil design, and its current-carrying surface is mainly a developable surface. Applying the SFM to undevelopable surfaces via the extrinsic method has also been discussed in detail in Ren’s paper [
13]. This work was aimed at coil designs on undevelopable surfaces via the intrinsic method, but it is also compatible with developable surfaces. The SFM on non-developable surfaces is an important prerequisite for the implementation of this work. Therefore, we provide here a brief explanation of the SFM, as well as the issues related to this article.
The SFM simplifies the complex calculations in direct problems by converting the current density vector into a scalar SF. The conversion relationship is represented by the following formula:
where
represents the surface current density;
is the tangential gradient operator on the current-carrying surface;
is a scalar stream function that is defined on the surface; and
represents the unit normal vector of the surface.
The core of the SFM is to find a suitable scalar function
, such that that the corresponding surface current density structure can produce a magnetic field that meets design requirements. This can be expressed as an optimization function:
Equation (2) shows a standard optimization problem, where
is the main objective and
is a regularization term used to solve ill-posed problems [
10,
14]. And
is a function of
, representing the actual magnetic field generated by the coils in the imaging region;
is a set of constants representing the target gradient magnetic field within the imaging region;
represents the auxiliary objective; and
is the regularization parameter, which is a key factor determining the smoothness of the coil.
We discretize the current-carrying surface in Equation (2) and choose the coil power consumption as the auxiliary objective. Thus,
and
, where
is the sensitivity matrix of magnetic field
with respect to
, and
is the power consumption matrix, which uses coil energy consumption as the auxiliary objective. The discrete expression of the above equation is as follows:
We can obtain the solution of Equation (3) by setting the first-order derivative
of the SF to zero. In this step, the choice of the optimal regularization parameter plays a crucial role in determining the weight balance between the two optimization objectives. Solutions obtained with smaller values of
may achieve lower main objective values, but they may also result in coil structures with a poor smoothness. On the other hand, larger values of
tend to make the coils smoother but introduce more errors in magnetic field linearity:
Therefore, in each instance of SFM, we need to balance the solution by adjusting the value of so as to balance magnetic field accuracy and coil power consumption. Thus, a smooth coil configuration can be obtained.
In order to choose an appropriate regularization parameter, it is often necessary to sample it multiple times and obtain an L-curve, which shows the relationship between the main and auxiliary objectives for the different values of the regularization parameter. In this way, we can select the value that satisfies the design objective from around the inflection point of the L-curve. Although this method can achieve the purpose of coil optimization design, it results in iteration calculations that uses multiples, which greatly increases the time and labor cost of the coil design. In contrast, smoothing methods can make the SF configuration evolve towards smoothness through linear calculations. Using this method, it only takes a short time to evolve a poorly smooth initial SF configuration to a smooth SF configuration.
2.2. Smooth Function on a Surface
The smoothness of the SF contour has always been an unavoidable research topic in coil design. The Tikhonov regularization method is a solution to this topic [
10,
14,
15]. Another method for smoothing SF contours is to use a smoothing operator, where the simplest form is the filtering method, which averages the stream function values of each node by applying certain weights [
16,
17]. The weights are usually assigned according to the distance information between the node and its neighbors in a well-defined neighborhood. This method changes the interpolation of the equivalent points within the unit by filtering the SF values of the nodes; thus, the SF contour within the design surface tends to be smooth. Before conducting the work in this paper, we attempted to use distance weighting for function filtering on non-developable surfaces, but we were unable to obtain a smooth contour solution. This is because the filtering matrix composed of distance weighting only includes the relative position information of the points on the surface and ignores the curvature degree of the surface itself.
Therefore, a more general smoothing operator is preferred that could include the geometric information of the current-carrying surface, such that the smoothing of the SF within the surface can be conducted inside only a manifold. An important property of SF is that it is embedded in the current-carrying surface. In the process of smoothing evolution, the smoothing operator needs to preserve the embedding property and use the geometric characteristics of the current-carrying surface to smooth the SF. In this way, the structure of the SF contour lines along the surface scale can be simplified [
18].
We let
be a two-dimensional manifold embedded in
;
be the coordinates in a two-dimensional manifold; and let
be the SF defined on the surface Γ [
19]. We then need to evolve the implicit function
such that its contour line,
, evolves toward a geodesic curve, where
is the SF value at a certain point
on the surface. Such problems have been described when using the heat diffusion equation that balances changes in the concentration in space [
20]:
where
represents the time of diffusion. Essentially, this involves using a filtering method with a second-order PDE. In the non-developable surface designed in this work, we need to extend the Laplace operator onto two-dimensional manifolds, where the Laplacian operator
needs to be replaced by the Laplace–Beltrami operator
[
21]. Then, the diffusion equation for functions on the surface
is
Equation (6) needs to be discretized in the time domain as follows:
where
represents the time step. After a simple derivation, the following form can be obtained:
Thus, we obtain the smoothing operator:
This means that we only need to solve a linear equation to achieve the smoothing of the surface functions. Here,
is the identity matrix,
is the time step, and
denotes the function value on the two-dimensional manifold
at time
in the time domain. The Laplace–Beltrami operator in the discretized surface representation can be calculated using the following equation [
21]:
where
and
are the angles of the triangle formed by the vertices of the mesh edge
and
is the actual neighborhood area under the control of node
, as shown in
Figure 2. For ease of illustration, a planar triangular mesh is shown in
Figure 2. In practice, the L–B operator is usually applied in two-dimensional manifolds, that is,
for a convex surface.
To speed up the smoothing process, this work normalizes the operator by using Desbrun’s method [
17], which allows for the use of larger time steps in explicit integration methods, as shown in Equation (10):
2.3. Smoothing Coefficients Based on the Objective Function Control
In the previous section, we normalized the time step for each iteration. However, in practical applications, there is often an abrupt decrease in the accuracy of the magnetic field during the initial stages of the smoothing iteration. Maintaining a constant smoothing step size can easily lose the documented changes occurring during the iteration process, resulting in a rapid loss of accuracy in the coil’s magnetic field. Therefore, the main objective function of the coil design process in this work is introduced to assist in the selection of the smoothing step size:
where
represents the main objective value for the nth iteration step. Before each iteration in the program, an estimate of the smoothed main objective value
is calculated. Then, the program determines the value of the iteration step
based on the size of the objective value in reverse.
That is, the program adjusts the size of the time step for each iteration in real time during the smoothing process. If the increase in the objective function is too large, a smaller step size coefficient (no less than 0.01) is assigned. If the increase in the objective function is relatively small, a larger step size coefficient (no greater than 1) is assigned to accelerate the iteration. Certainly, the objective function may be decreased during the smoothing process. In this case, we assign the maximum step size coefficient, which is 1.
2.4. Controlling the Spacing between Contour Lines Based on the Tangential Gradient
In most cases of gradient coil design and manufacturing, we want to maintain a large number of wire turns to reduce the current value of the wires and maintain the stability of the magnetic field formed by the coil. However, in the case of a high number of wire turns, some areas will have wire layouts that are too dense. This can easily cause an increase in the inductance of the coil. Due to the frequent switching of gradient coils during the operation of MRI equipment, coils generating the gradient magnetic fields should have low inductance [
22,
23]. The method of reducing the density of the coil is an empirical way to reduce inductance. As mentioned earlier, the contour lines in the SFM are the coil structure. Therefore, reducing the sampling interval of the contour lines can easily reduce the density of the coil. However, this will also decrease the coil density in already sparse areas. Of course, we do not want to excessively sacrifice the accuracy and stability of the magnetic field to reduce the inductance. Fortunately, during the diffusion process of the function, the coil layout naturally evolves toward uniformity. All we need to do is monitor the changes in wire spacing in real time during the design process. Thus, given the number of coil turns, we can select a smooth result that satisfies the design conditions while maximizing the wire spacing between the wires. This makes it an important indicator for assisting in the selection of a smooth result. In this work, the gradient operator
is used to perform first-order gradient calculations on the SF. The theoretical minimum wire spacing that can be achieved under the current result can be estimated by combining the numerical value of the SF with the number of discrete wires as follows:
where the
is the number of discrete coil turns and
is the tangential gradient vector of the SF along the surface. By using the tangential gradient operator, we can monitor the changes in wire spacing during the smoothing process of the SF on any surface. This allows us obtention of the optimal number of discrete coil turns that meet the design goals.
2.5. Curvature Changes in Implicit Contour Expression
The goal of this work is to smooth the SF contour lines. Thus, this section describes a quantitative index for smoothing. The smoothness of a curve in space at any point can be expressed by its curvature. A higher curvature at a point indicates that the curve is more bent at that position, while a lower curvature at a point indicates a smoother curve. This section calculates the implicit curvature value of the contour at each node on the current-carrying surface. Thus, the maximum value of the implicit curvature is used as the index for the level of function smoothness.
In the SFM, the contour lines of the SF at any point inside the coil design surface
can be implicitly expressed by the function
, where
can take on any value of the SF within the design domain and where
. Thus, the curvature value of the curve can be calculated by using a second-order differential of its implicit function expression:
During the smoothing process, the maximum curvature value
on the current-carrying design domain gradually decreases and converges. Based on this, the convergence condition for smoothing can be set to
where
represents the maximum curvature value of the SF contour configuration in the i-th iteration.