Akiko Sasaki

Julien Trébosc

Hiroki Nagashima

Jean-Paul Amoureux

On the applicability of cosine-modulated pulses for high-resolution solid-state NMR of quadrupolar nuclei with spin > 3/2

In MQMAS-based high-resolution solid-state NMR experiments of half-integer spin quadrupolar nuclei, the high radiofrequency (RF) field requirement for the MQ excitation and conversion steps with two hard-pulses is often a sensitivity limiting factor in many practical applications. Recently, the use of two cosinemodulated (cos) low-power (lp) pulses, lasting one-rotor period each, was successfully introduced for efficient MQ excitation and conversion of spin-3/2 nuclei with a reduced RF amplitude. In this study, we extend our previous investigations of spin-3/2 nuclei to systems with higher spin values and discuss the applicability of coslp-MQ excitation and conversion in MQMAS and MQ-HETCOR experiments under slow and fast spinning conditions. For the numerical simulations and experiments we used a moderate magnetic field of 14.1 T. Two spin-5/2 nuclei ( 85 Rb and 27 Al) are mainly employed with a large variety of CQ values, but we show that the practical set up is also available for higher spin values, such as spin-9/2 with 93 Nb in Cs4Nb11O30. We demonstrate for nuclei with spin value larger than 3/2 a preferential use of coslp-MQ acquisition for low-gamma nuclei and/or large CQ values with a much reduced RF-field with respect to that of hard-pulses used with conventional methods.

I. Introduction

Despite that nuclei with spin value I > ½ account for more than 70% of NMR-active isotopes, many practical applications of these nuclei are severely limited due to the presence of the quadrupolar interactions, which cause spectral broadenings and thus result in a significant loss of sensitivity and resolution. To eliminate this broadening, which is proportional to the quadrupolar coupling constant, CQ = e 2 qQ, magic angle spinning (MAS) 1 is routinely employed in solid-state NMR (ssNMR). For half-integer spin quadrupolar nuclei (I = 3/2, 5/2, 7/2 and 9/2), high-resolution two-dimensional (2D) experiments, such as multiple-quantum MAS (MQMAS) 2 or satellite-transition MAS (STMAS), 3 are performed to obtain truly isotropic spectra. Although MQMAS is more widely used than STMAS as a routine method owing to its ease in experimental setup, these two methods are complementary to each other as their comparison allows to identify the potential presence of dynamics close to the quadrupolar nuclei. 4,5 Moreover, following the development of fast MAS probes with small-diameter rotors, 1 H-detected HETCOR (hetero-nuclear correlation) experiments with quadrupolar nuclei have been proposed, [6][7][8][9][10][11][12][13] which utilize the improved 1 H resolution under fast spinning. 14,15 The acquisition of MQMAS/STMAS and MQ/ST-HETCOR 2D spectra allows a detailed analysis of the local structure around the quadrupolar nuclei of interest, in much the same way that has been employed for I = 1/2 spin systems.

In the MQ-based high-resolution experiments, the intrinsically low efficiency of the MQ excitation/conversion steps is often a sensitivity-limiting factor in practical applications. This is because the coherences on the three-quantum level are 'forbidden' from the quantum mechanics point of view, and therefore the MQ excitation and conversion processes with conventional hard-pulses (hp) require a high radiofrequency (RF) field, 16 which moreover depends on the CQ value.

Such large RFs are not always attainable, especially with low-gamma nuclei, or when large diameter rotors (  3mm) are used for insensitive isotopes. To circumvent this high RF requirement of the conventional hard-pulses, a variety of experimental approaches has been implemented, such as RIACT, 17 FAM, 18,19 DFS, 20 FASTER, 21 HS, 22 RAPT, 23 and QCPMG-MQMAS, 24 although they often require a time-consuming optimization process of multiple parameters.

Recently, two novel variations for MQ excitation/conversion were devised, namely lp-MQMAS 25 and coslp-MQMAS, 26 which utilize two identical low-power long-pulses (lp) of one rotor period (R) each, enabling the MQ excitation/conversion with a reduced RF field. These two approaches stem from the lp-STMAS scheme, 27 previously proposed to manipulate the innermost STs with low RF requirement. In this article and in the previous one, 31 we use a notation of the sequences that slightly differs from the original ones. 25,26 We call them lp-MQMAS and coslp-MQMAS, instead of lpMQMAS and cos-lpMQMAS, to emphasize the fact that they are MQMAS methods, and we add WURST at the beginning to emphasize the fact that this sweeping is a very important part of these sequences.

In coslp-MQMAS, the amplitude of the two identical R-pulses is cosine-modulated, and it has been shown that with respect to lp-MQMAS, which uses two square pulses, the efficiency is about 30% higher and the isotropic linewidths are narrower and hence more resolved owing to its robustness to spinning instability. 26 Since the lp-and coslp-MQMAS sequences begin with a centraltransition (CT) selective 90 pulse, they both benefit from the existing CT signal enhancement schemes. Among those well-established approaches, the WURST (wideband, uniform-rate, smooth-truncation) [28][29][30] scheme was chosen for the ease of implementation in the aforementioned lp-and coslp-MQMAS studies. 25,26 Recently, we have made a thorough investigation on I = 3/2 nuclei of the practical aspects of MQMAS and MQ-HETCOR experiments under slow and fast MAS conditions, using 87 Rb, 71 Ga and 35 Cl isotopes. Our study concluded that WURST-coslp-MQbased experiments are highly recommended for (i) low-gamma nuclei, (ii) large rotor diameters, and (iii) samples with large CQ values. 31 For I > 3/2 spin systems, at the time of writing, only 17 O WURST-coslp-MQMAS results have been reported in the acquisition of 3D 1 H/ 13 C/ 17 O, 32 or 2D 1 H/ 17 O, 33 isotropic HETCOR spectra at B0 = 18.8 T.

In this study, we extend our previous investigations of spin-3/2 nuclei to higher spin systems, mostly to I = 5/2, and discuss the applicability of coslp-MQ excitation/conversion in the context of MQMAS and MQ-HETCOR under slow and fast spinning conditions. All experiments and simulations presented in this article have used a moderate magnetic field of B0 = 14.1 T. Two I = 5/2 isotopes, 85 Rb and 27 Al, have been employed with a large range of CQ values. With the aid of numerical simulations, the robustness of WURST sensitivity enhancements has been verified for I = 5/2 nuclei. Using 85 Rb in RbNO3 (CQ = 3-4 MHz) 34,35 and Rb2SO4 (CQ = 5 and 11 MHz), 36 a preferential use of WURST-coslp-MQMAS over the conventional approaches is demonstrated for low-gamma nuclei and either large CQ values or rotor diameters. Additionally, using 27 Al in AlPO4-berlinite (CQ = 4 MHz), 35,37 Al(acac)3 (Al acetylacetonate: CQ = 3 MHz), 38 Al(lact)3 (Al lactate: CQ = 5 MHz), 15 and ipa-AlPO4-14 (CQ = 1.8-5.6 MHz), 5,39,40 we briefly comment on the z-filter and shifted-echo versions of coslp-MQMAS acquisition of I = 5/2 nuclei and also on a possible use of { 27 Al}-1 H WURST-coslp-MQ-HETCOR experiments at R = 62.5 kHz. We then extend our analysis to a higher spin value, using 93 Nb (I = 9/2) in Cs4Nb11O30 (CQ  15 MHz). 35,41 In the final section, we account for the differences of coslp-MQ behaviors between I = 3/2 and 5/2 spin systems, using numerical simulations with 23 Na and 27 Al nuclei.

II. Methods

All experiments were performed using a Bruker Avance NEO spectrometer with a B0 = 14. The MQMAS/STMAS pulse sequences for I > 3/2 spin systems used in this study are summarized in Figs. 1 andS1, for split-t1 shifted-echo (or full-echo) and z-filter acquisitions, respectively, and in Fig. S2 for {I}-1 H MQ/ST-HETCOR experiments. It must be noted that even if the number and type (hp, lp, coslp or CT-selective) of rf-pulses in the sequences do not depend on the spin value, the timings and selected coherence levels used in the pulse programs are different between those used with I = 3/2, 31 or I > 3/2 (here). Examples of TopSpin pulse programs for WURST-coslp-MQMAS and WURST-coslp-MQ-HETCOR experiments are provided in the SI of this article for I > 3/2 and in that of our previous study for I = 3/2. 31 Prior to any STMAS acquisition, an accurate adjustment of the spinning axis to the magic angle was performed using a DQF version of the split-t1 shifted-echo pulse sequence on the sample of interest itself. 42 Chemical shift scales, shown in ppm, were referenced using the sample of interest itself as a secondary reference, and all 2D spectra were referenced according to the unified representation. 43 All simulations were performed using the SIMPSON program. 44 The quadrupolar interaction up to the second-order was taken into account with Q = 0, without scalar or dipolar coupling. The input variables were, the nucleus of interest, the spinning frequency (R), the RF amplitude (1) and offset (off), and the quadrupolar coupling constant (CQ). To calculate the transfer efficiency, coherences were chosen to emulate the pulse sequences without any delay, except for the rotor-synchronized STMAS experiments (t1min = R -p1/2 -p4 -p2/2). In high-resolution methods, the efficiency is defined as the ratio between the MQMAS or STMAS CT final signal and that observed after a CT-selective 90° pulse. The powder averaging parameters (i.e. crystal file and number of -angles) and the maximum time step t over which the Hamiltonian is considered time-dependent were tested for convergence, and a combination of ZCW54, 45- 47 with 10 -angles and t = 0.1 s was sufficient for the given range of CQ values used in the plots. Since both the real (Re) and imaginary (Im) parts of the signal were found to be significant, its magnitude ((Re 2 + Im 2 )) is plotted as the signal intensity, unless stated otherwise. Further simulation details are given in the figure captions.

III. Results and discussion

We first briefly remind of a few essential practical points behind these WURST-lp/coslp-MQ-based experiments established in the previous studies of I = 3/2 spin systems. [25][26][27]31 (i) For MQ excitation/conversion, the two R-pulses need to be identical and symmetrically applied to produce a coherent evolution of the two inner STs, 25,26 to serve as inversion pulses between the 1Q-CT and 3Q coherences. (ii) In lp-MQMAS, the R-pulses are applied far off-resonance from the CT signals with a large offset (irr) on the ST manifold, 25 whereas in coslp-MQMAS, the R-pulses are applied near to the CT signals (irr), ideally on the center-band position of the inner STs, and are cosine-modulated (cos). 26 (iii) The optimum RF field of the R-pulses depends on the spinning frequency and the spin and CQ values, and its amplitude in coslp-MQMAS, 1,opt  [CQR/(2I(2I -1))] 1/2 , (1) is approximately 2 higher compared to lp-MQMAS. 26,31 In the following paragraphs, these technically essential points are similarly illustrated using I = 5/2 spin systems. 85 Rb is a I = 3/2 low-gamma nucleus (0  58 MHz at 14.1 T) with a natural abundance of 72.2 % and along with the 87 Rb other isotope (I = 5/2), 85 Rb has been employed in ssNMR method development owing to its favourable relaxation properties (e.g. T1  a few hundred ms). 36 The sensitivity-limiting factor in conventional MQMAS/STMAS signal acquisition of such low-gamma nuclei is the small maximum RF field attainable for the hard-pulses by the conventional ssNMR setup. Here, using 85 RbNO3 and 85 Rb2SO4, we demonstrate the sensitivity advantage of WURST-coslp-MQMAS acquisition for low-gamma nuclei within the maximum RF field achievable (43 kHz) by our 3.2 mm MAS probes. Three 85 Rb species exist in RbNO3 and two in Rb2SO4, with (cs (ppm), CQ (MHz), Q) = (-27.4, 3.3, 0.2), (-28.9, 4.2, 0.9), (-31.6, 3.4, 0.6) and (42, 5, 0.9), (16, 11.0, 0.1). 34,35 Fig. 2 shows a comparison of the isotropic projections of 2D 85 Rb DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra of RbNO3 and Rb2SO4 with z-filter and full-echo acquisition. We observe that the WURST-coslp-MQMAS sensitivity is 3-4 times higher than the conventional hp-MQMAS and as good as the STMAS equivalent, especially with full-echo acquisition. Moreover, the best resolution is always observed with WURST-coslp-MQMAS, as best seen for the two close resonances of RbNO3. For crystalline compounds, the full-echo acquisition shows a sensitivity advantage over the z-filter equivalent, especially for samples with large CQ values and/or small losses (few or no 1 H nuclei). 48 Thanks to these advantages of sensitivity and resolution, WURST-coslp-MQMAS allows with full-echo acquisition the observation in Rb2SO4 of the species with CQ = 11 MHz at ca. -140 ppm. This was impossible with all other sequences. It should be emphasized that the WURST sensitivity enhancement gives rise to a factor of 2.5-3.0 increase in the lp/coslp-MQMAS intensity (not shown). This is verified via numerical simulations (Fig. S3), which summarize at R = 21 kHz the simulated 85 Rb WURST signal gain for RbNO3 and Rb2SO4, with respect to the WURST RF field (1,WURST) and offset frequency (off). It can be seen that WURST is robust with respect to RF field and offset, given the range of 1,WURST = 10-20 kHz and off > 200 kHz. It must be noted that on the contrary, WURST sweeps do not bring any signal enhancement with hp-MQMAS.

III

Similarly, Fig. 3 shows the simulated 85 Rb signal intensity profiles of the same four MQMAS/STMAS sequences with respect to CQ and 1 RF-field of the MQ/ST excitation/conversion pulses. In consistence with the experimental results, the conventional hp-MQMAS is globally inefficient within the practically attainable 85 Rb RF range limited to 1 < 43 kHz for  = 3.2 mm, whereas WURSTcoslp-MQMAS performs even better than the DQF-STMAS counterpart, owing to the intrinsically low RF requirement of the Rpulses.

Here, it is worth reminding that, upon setting up the R-pulses of the coslp-MQMAS experiment, the RF field (1) and the cosine modulation (cos) and carrier (irr) frequencies need being optimized on the sample, especially with large CQ values. 26,31 The first two of these points are illustrated in Fig. S4, using the experimental 85 Rb WURST-coslp-MQMAS signal intensity profiles of RbNO3 and Rb2SO4. For these two compounds, cos = 100-200 kHz was found to be optimum and, since CQ is larger in Rb2SO4 (CQ = 5 and 11 MHz) than RbNO3 (CQ = 3.3-4.2 MHz), the optimum RF is larger in Rb2SO4 (1  30 kHz) than RbNO3 (1  20 kHz) (see Eq.1). For the third point, the offset positioning is more adverse in Rb2SO4 than RbNO3 because the inner ST center-band lies far away from the CT for large CQ sites. This is illustrated in Fig. S5, by a comparison of experimental 85 Rb WURST-coslp-MQMAS signal intensity profiles of RbNO3 and Rb2SO4, with respect to the irradiation frequency (irr). In RbNO3, changing the offset position of Rpulses by 8 kHz results in no apparent change in intensity, whereas in Rb2SO4, irr  8 kHz needs to be applied for the best sensitivity.

Overall, our experiments and numerical simulations have verified that the WURST-coslp-MQMAS sensitivity is 3-4 times higher than the conventional hp-MQMAS for 85 Rb, and even better than the STMAS equivalent without the precise adjustment of spinning axis to the magic angle. Hence, WURST-coslp-MQMAS is highly recommended for low-gamma nuclei, for which the large RF requirement of conventional hp-MQMAS cannot be satisfied and/or with large CQ samples for which the MQ excitation/conversion efficiency is expected to be intrinsically low.

III.2. 27 Al NMR of AlPO4-Berlinite (CQ = 4.0 MHz) at R = 16 kHz

Having demonstrated the sensitivity advantage of WURST-coslp-MQMAS with the low-gamma 85 Rb isotope, we envisaged to extend the use of cosine-modulation to moderate-gamma nuclei, 27 Al, aiming to extrapolate the results to { 27 Al}-1 H MQ-HETCOR experiments under fast MAS conditions.

First, we employed a single 27 Al site compound, AlPO4-Berlinite (CQ = 4.0 MHz), 35,37 at R = 16 kHz to check the 27 Al WURSTcoslp-MQMAS sensitivity in comparison with other conventional methods. The isotropic projections of 2D 27 Al DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra of AlPO4-berlinite are summarized, with z-filter and full-echo acquisitions, in Fig. S6.

As often verified for samples with small-moderate CQ values, the full-echo approach is here less favourable than the z-filter counterpart, due to the required long echo delay (echo = 5 ms). With this high Larmor frequency (156.4 MHz), a sufficiently high RF-field (1 = 61 kHz) is achievable with  = 3.2 mm rotors, and hence STMAS is best for sensitivity, on condition that the spinning axis is accurately set to the magic angle prior to the signal acquisition. However, although not as good as STMAS, the WURST-coslp-MQMAS sensitivity was found better than the conventional hp-MQMAS, with a very much reduced RF requirement (1  18 instead of 61 kHz). We then chose a 1:1 molar mixture of Al(acac)3 + Al(lact)3, each containing one single 27 Al site, with CQ = 3.0 and 5.0 MHz, respectively. 15,38 By observing the 27 Al isotropic resonances, we reckon that Al(acac)3 exhibits a crystalline nature with a sharp, well-defined peak, while Al(lact)3 has a more distributed character with a broad, featureless peak. This is reflected in Fig. 4, which shows the isotropic projections of this mixture recorded with 2D 27 Al DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS sequences at R = 16 kHz, with z-filter and full-echo acquisitions.

These spectra are fully in agreement with those observed with AlPO4-Berlinite: (i) the z-filter acquisition is better than full-echo for these protonated samples with moderate CQs, (ii) DQF-STMAS shows the highest sensitivity, but with 1 = 61 kHz and a perfectly adjusted magic-angle, and (iii) for Al(acac)3 WURST-coslp-MQMAS with only 1 = 15 kHz is better than the conventional hp-MQMAS with 61 kHz. For Al(lact)3, the WURST-coslp-MQMAS sensitivity gain is slightly smaller due to the distribution of surrounding.

III.4. 27 Al NMR of ipa-AlPO4-14 (CQ = 1.8-5.4 MHz) at R = 62.5 kHz

Aiming for { 27 Al}-1 H WURST-coslp-MQ-HETCOR detection, we then switched to a small rotor diameter,  = 1.3 mm, and performed the equivalent set of experiments with R = 62.5 kHz, using ipa-AlPO4-14. This compound contains four 27 Al sites showing isotropic peaks from left to right: (cs (ppm), CQ (MHz), Q) = (42.2, 1.8, 0.6), (45.0, 4.1, 0.8), (27.0, 5.4, 0.9), (-0.2, 2.6, 0.9). It has been known that the STMAS signals of ipa-AlPO4-14 are motionally broadened due to the presence of s dynamics around Al nuclei, 5 which makes high-resolution MQMAS and STMAS acquisitions complementary to each other. For coslp-MQMAS, it is worth remembering that the optimum RF-field of the R-pulses depends not only on the CQ value, but also on the spinning frequency (Eq.1). 26,31 This is verified in our 27 Al experimental observations, shown in Figs. 4,5, and S6, where the optimum RF-field of the Rpulses increased from 1  18 to 30-66 kHz, as the spinning frequency increased from R = 16 to 62.5 kHz.

Upon WURST-coslp-MQMAS acquisition of ipa-AlPO4-14, we observed two additional complexities.

First, for such samples containing multiple 27 Al sites with a wide range of CQ values, the coslp-MQMAS variables need to be chosen as a compromise over the optimum ranges for all sites, which consequently led to only a marginal increase in the sensitivity overall. As a result, the signal intensities observed with WURST coslp-MQMAS are smaller than with DQF-STMAS, and they are similar to those recorded with hp-MQMAS (Fig. 5). However, the required RF-field was much smaller with WURST-coslp-MQMAS (30 or 66 kHz) than with DQF-STMAS and hp-MQMAS (121 kHz), and no stringent magic angle set up was required.

Second, when phasing a 2D full-echo hp-MQMAS spectrum, one can usually use the same phase correction parameters as for the analogous 1D spectrum from which the 2D dataset is created. However, when phasing the 2D full-echo WURST-coslp-MQMAS spectra, we needed to perform the 2D phase correction from scratch (i.e. discarding the phase parameters from the 1D dataset) to find the parameters that simultaneously produce a properly-phased 2D spectrum along the MAS and isotropic dimensions. We suspect that this 2D phasing complexity originates from the phase evolution during the R-pulses. This phase issue did not arise in AlPO4-berlinite with a single Al site or Al(acac)3 + Al(lact)3 with close CQ values, but was noticeable with ipa-AlPO4-14 containing multiple 27 Al sites with a wide range of CQ values. Based on our hands-on experience, we suggest the following step-by-step phasing protocol when encountering such phasing issues upon full-echo acquisition: 1) reset all phases. 2) Calculate an approximate firstorder correction in F2, such as 180*(echo + 3R)/(dwell-time). Here, the additional evolution on 3Q during the long R pulses must be considered. 3) Do the 2D Fourier transform. 4) Phase the F1 dimension to have two neighbouring peaks well phased. 5) Correct the first-order term in F2 by a small amount. 6) Repeat steps 3 to 5 until all peaks are well phased. It must be noted that this way of phasing is very close to that described previously for 2D lp-MQMAS spectra. 25 Fortunately, this phase issue is less critical in the 2D { 27 Al}-1 H MQ/ST-HETCOR spectra, because the signals are narrow along the two dimensions with high spinning speeds. Fig. 6 compares the 1D { 27 Al}-1 H MQ/ST-HETCOR 1 H spectra of ipa-AlPO4-14 at R = 62.5 kHz, recorded with ST and ST-SPAM, 15 hp-MQ and hp-MQ-SPAM, 14,15 WURST-lp-MQ and WURST-coslp-MQ approaches. The SPAM (Soft-Pulse Added-Mixing) scheme has been used here, because it increases the efficiency of the second hard-pulse and only differs from the conventional hp-MQMAS and STMAS approaches by a change in the phases. [48][49][50] Due to the fact that the coslp-MQ variables result from a compromise over the wide CQ range, we observed no sensitivity advantage in { 27 Al}-1 H WURST-coslp-MQ-HETCOR acquisition of ipa-AlPO4-14. It should be remembered that, as in STMAS and STMAS-SPAM, the isotropic dimension of ST-HETCOR spectra of ipa-AlPO4-14 (Fig. 5) is broadened by s dynamics around Al nuclei, 5 and thus the required accuracy of magic angle setting is then not as stringent as expected for other crystalline samples. Therefore, in this particular case, we recommend the use of the two MQ/ST-SPAM-HETCOR acquisitions, for the ease of setup and the best sensitivity, in spite of the higher RF-field than with WURST-coslp-MQ-HETCOR (121 v.s. 54 kHz),

To summarize, for high-gamma nuclei (e.g. 27 Al) with moderately large CQ values (e.g. > 4 MHz), 27 Al WURST-coslp-MQMAS can perform slightly better than the conventional hp-MQMAS, but with a much reduced RF requirement. For { 27 Al}-1 H MQ-HETCOR under fast MAS, MQ/ST-SPAM-HETCORs may be the methods of choice for the ease of setup, especially when multiple sites with a large variety of CQ values are present. In a last step, we have extended the R-pulses high-resolution methods to higher spin-value nuclei, and we have chosen 93 Nb (I = 9/2) of Cs4Nb11O30, which contains a single species with CQ  15 MHz and Q  0.6. 35, 41 We have used a fast spinning speed of R = 62.5 kHz, to minimize the numerous spinning sidebands of this broad second-order MAS spectrum. In Fig. 7, we show the isotropic projections of the 93 Nb 2D DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra, with z-filter and full-echo acquisitions. The best sensitivity is observed with DQF-STMAS, but this experiment requires a perfect setting of the magic angle, a very stable spinning speed and a high RF-field (1 = 100 kHz). The second best sensitivity is observed with WURST-coslp-MQMAS, especially with the z-filter acquisition. This robust method is more efficient than the hp-MQMAS sequence and only requires a small RF-field of 14 kHz, instead of 100 kHz for the hp-MQMAS version.

III.6. Comparison between 23 Na (I = 3/2) and 27 Al (I = 5/2) WURST-coslp-MQMAS at R = 20 kHz

In the previous sections, we have demonstrated that for nuclei with spin value higher than 3/2, WURST-coslp-MQMAS always requires a small RF-field and is hence simultaneously favorable from the sensitivity point of view in the cases of (i) low-gamma nuclei, (ii) large rotor diameters used for insensitive nuclei, and (iii) large CQ sites. However, its sensitivity is less well performing with MQMAS and MQ-HETCOR acquisitions of high-gamma nuclei. Here, we remind of our previous I = 3/2 investigations, 31 which revealed a general large sensitivity advantage of WURST-coslp-MQMAS based acquisitions, not only for low-gamma isotopes (e.g. 35 Cl) or large CQ sites (e.g. 71 Ga), but also for high-gamma 87 Rb and 23 Na MQMAS and MQ-HETCOR experiments. In this last section, we briefly account the difference in WURST-coslp-MQMAS behaviors between I = 3/2 and 5/2 spin systems, using numerical simulations of 23 Na and 27 Al nuclei. We chose this pair of nuclei as their Larmor frequencies are close to each other (158.7 and 156.4 MHz at 14.1 T, respectively).

In Fig. S7 we show the simulated 23 Na and 27 Al WURST gains at R = 20 kHz, with respect to the WURST RF amplitude (1,WURST) and offset (off). We used CQ values of 2 and 6.6 MHz for 23 Na and 27 Al, respectively, to account for the factor defining the quadrupolar interaction with respect to Q = 3CQ/(2I(2I -1)). We observe that the combination of WURST with off = 300 and 1 = 20 kHz can be safely employed for both nuclei. We also note that the WURST gain increases from ca. 2.0 to 2.7 when going from a spin-3/2 to a spin-5/2 nucleus. This rise is related to the increased number of STs; the higher the spin value the larger the WURST gain.

In Fig. 8, we summarize the simulated 23 Na (I = 3/2) and 27 Al (I = 5/2) signal intensities of hp-MQMAS and WURST-coslp-MQMAS at R = 20 kHz, with respect to CQ and the RF-field (1) of the excitation/conversion pulses.

First, upon comparison of the conventional hp-MQMAS results, one may note that I = 5/2 requires a lower minimum RF-field than I = 3/2, ca. 40 instead of 70 kHz, for a relatively effective 3Q excitation/conversion. This point partly contributes to the favorable comparison of the WURST-coslp-MQMAS performance of 3/2 spin over the 5/2 equivalent, as the sensitivity is generally compared with respect to the conventional hp-MQMAS.

Second, the overall WURST-coslp-MQMAS efficiency is globally twice higher for I = 3/2 than for 5/2. The efficiency of the coslp-MQMAS transfer is related to the coherent inversions of the two inner STs by the R-pulses, such that the initial 3Q coherences created by the first R-pulse are refocused by the second one. 25,26 For I = 3/2 spin systems, only these inner STs exist, and thus the coherences are well-confined, and this refocusing is consequently highly efficient. The signal is observable with an RF-field of ca. 20 kHz, and is maximum at about 70 kHz, with intensities of ca. 0.6 and 0.9, respectively (Fig. 8b). For I > 3/2 nuclei, on the contrary, due to the presence of other STs, a coherence leakage to these STs is inevitably observed during the selective inversion of the inner-most STs. The signal starts with an RF-field of ca. 10 kHz, and is maximum with an intensity of ca. 0.5 at about 20 kHz. For higher RF-fields, the leakage process largely decreases the sensitivity (Fig. 8d).

IV. Conclusions

Following the recent development of WURST-coslp-MQMAS, which enables an efficient MQ excitation/conversion of I = 3/2 nuclei with a reduced RF-field, we extend our previous discussions to higher spin values in the context of high-resolution MQMAS and MQ-HETCOR experiments.

As general rules, we have shown that for all spin values the WURST CT enhancement is effective leading to a signal gain of ca. 2-3, and that the WURST-coslp-MQMAS acquisition (i) does not have to be rotor-synchronized, contrary to STMAS, to expand the spectral width, as example when using large rotor diameters, (ii) requires a small RF-amplitude, contrary to STMAS and hp-MQMAS methods, and (iii) enhances the resolution with respect to other conventional methods.

The sensitivity advantage for spin-5/2 nuclei of WURST-coslp-MQMAS acquisition has been experimentally demonstrated for lowgamma nuclei and/or large CQ values using 85 RbNO3 and 85 Rb2SO4.

For high-gamma nuclei with spin larger than 3/2 (e.g. 27 Al) and with moderate CQ values, the use of WURST-coslp-MQMAS and {I}-1 H WURST-coslp-MQ-HETCOR experiments under fast MAS conditions has been demonstrated. However, the sensitivity advantage related to the two 1Q-CT  3Q coslp transfers is then largely decreased, especially in the case of samples with very different CQ values, and other approaches based on composite pulses, eg. hp + SPAM, are more efficient.

The WURST-coslp-MQMAS high resolution method has been demonstrated with higher spin values, i.e. I = 9/2 with 93 Nb isotope, and the method has also shown a small RF requirement.

The difference in WURST-coslp-MQ behaviors between I = 3/2 and I  5/2 spins has been elucidated using simulations on 23 Na and 27 Al nuclei. For I = 3/2 spin systems, only inner STs exist and the transfer of coherences is well-confined and hence efficient. For I  5/2 nuclei, on the contrary, due to the presence of other STs, a coherence leakage to these extra STs occurs during the selective inversion of the inner-most STs, and the efficiency is decreased. Nevertheless, the sensitivity advantage of this method, which uses weak RF-fields, offers a great promise for low-gamma nuclei, (ii) large rotor diameter rotors used with insensitive nuclei, or (iii) samples with large CQ values; for which conventional approaches fail to produce sufficient MQ signals with hard-pulses. We predict that WURST-coslp-MQMAS will expand the range of nuclei and CQ values that can be further investigated using high-resolution quadrupolar NMR.

. 1 .

 1 85 Rb NMR of RbNO3 (CQ = 3.3-4.2 MHz) and Rb2SO4 (CQ = 5 & 11 MHz) at R = 21 kHz

III. 3 .

 3 27 Al NMR of Al(acac)3 + Al(lact)3 (CQ = 3, 5 MHz) at R = 16 kHz

III. 5 .

 5 93 Nb (I = 9/2) NMR of Cs4Nb11O30 (CQ = 15 MHz) at R = 62.5 kHz

Fig. 1 .

 1 Fig.1. Pulse sequences and coherence transfer pathways for split-t1 shifted-echo (or full-echo) acquisition: (a) DQF-STMAS, (b) hp-MQMAS, (c) WURST-lp-MQMAS and (d) WURST-coslp-MQMAS experiments, where {RMQ, RST} = {19/12, 7/24} for I = 5/2. It must be noted that MQMAS experiments shown in (b-d), can be performed unsynchronized (t1  nR) to increase the indirect spectral width.

Fig. 2 .

 2 Fig.2. Comparison of the isotropic projections of 2D 85 Rb DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra of (a,b) RbNO3 and (c,d) Rb2SO4 with RD = 0.25 s, and (a,c) z-filter (Z) or (b,d) full-echo (FE) acquisition (echo = 4 ms), respectively. On each spectrum, the maximum intensity is indicated relative to that observed with DQF-STMAS. For the WURST coslp-MQMAS spectrum of Rb2SO4 acquired with full-echo, the second isotropic resonance at ca. -140 ppm is indicated with an arrow. WURST80: WURST = 1 ms, sweep = R = 21, 1,WURST = 13, off = 300 kHz. DQF-STMAS and hp-MQMAS: the RF-field for the two hard-pulses was fixed at its maximum safe value (43 kHz) for the probe: 1 {p1, p2, p3, p4} = {43, 43, 10, 10} kHz. RbNO3: Texp  0.8 h each. NS {MQ/ST} = {24/32 (Z), 48/64 (FE)}, t1 {MQ/ST} = {50/100} s, Nt1 = {384 (Z), 192 (FE)}. WURST-lp-MQMAS: 1 {Z, FE} = {23, 19}, irr {Z, FE} = {350, 400} kHz. WURST-coslp-MQMAS: 1 {Z, FE} = {19, 19}, irr = 0, cos = 150 kHz. Rb2SO4: Texp = 2.0 h each. NS {MQ/ST} = {192/192 (Z), 384/384 (FE)}, t1 {MQ/ST} = {45/90} s, Nt1 = {150 (Z), 75 (FE)}. WURST-lp-MQMAS: 1 {Z, FE} = {30, 23}, irr {Z, FE} = {400, 300} kHz. WURST-coslp-MQMAS: 1 {Z, FE} = {30, 30}, irr = 6, cos = 150 kHz.

Fig. 3 .

 3 Fig.3. Simulated 85 Rb signal intensity of DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS, with respect to CQ and RF-field (1) of the MQ/ST excitation and conversion pulses. DQF-STMAS/hp-MQMAS: 1 (p3 & p4) = 10 kHz, {p1, p2, p3, p4} = {1.8, 1.2, 8, 16}/{4.0, 1.8, 8, 16} s. WURST-lp-MQMAS: irr = +340 kHz. WURST-coslp-MQMAS: cos = 200 kHz. WURST80: WURST = 1 ms, sweep = R = 21, 1,WURST = 15, off = 300 kHz.

Fig. 4 .

 4 Fig.4. Al(acac)3 + Al(lact)3: comparison of isotropic projections of 2D 27 Al DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra with (a) z-filter and (b) full-echo (echo = 10 ms) acquisitions. RD = 0.5 s, Texp = 1.3-2.0 h each. On each spectrum the intensities are given relative to those with DQF-STMAS. hp-MQMAS/DQF-STMAS: NS = {48/32 (Z), 96/64 (FE)},t1 = {62.5/125} s, Nt1 = {300 (Z), 150 (FE)}, 1 {p1, p2, p3, p4} = {61, 61, 5, 5} kHz. WURST-lp-MQMAS: 1 {Z, FE} = 15, irr {Z, FE} = {150, 200} kHz. WURST-coslp-MQMAS: 1 = 18, irr = 0, cos = 200 kHz. WURST80: WURST = 1 ms, sweep = R = 16, 1,WURST = 10, off = 250 kHz.

Fig. 5 .

 5 Fig.5. ipa-AlPO4-14: comparison of isotropic projections of 2D 27 Al DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra with (a) z-filter and (b) full-echo (echo = 6 ms) acquisitions. RD = 0.25 s, Texp = 3.0 h each. On each spectrum, the intensity at iso  35 ppm is indicated relative to that observed with hp-MQMAS. The DQF-STMAS spectra are not included in the comparison due to the motional broadening. hp-MQMAS/DQF-STMAS: NS = {72/192 (Z), 144/384 (FE)}, t1 = {16/32} s, Nt1 = {600/200 (Z), 300/100 (FE)}, 1 {p1, p2, p3, p4} = {121, 121, 10, 10} kHz. WURST-lp-MQMAS: 1 {Z, FE} = {45, 30}, irr {Z, FE} = {400, 300} kHz. WURST-coslp-MQMAS: 1 {Z, FE} = {66, 30}, irr = 0, cos {Z, FE} = {400, 250} kHz. WURST80: WURST = 1 ms, sweep = R = 62.5, 1,WURST = 13, off = 250 kHz.

Fig. 6 .

 6 Fig.6. ipa-AlPO4-14: comparison of 1D { 27 Al}-1 H MQ/ST-HETCOR spectra with (a) ST and ST-SPAM, (b) hp-MQ and hp-MQ-SPAM, (c) WURST-lp-MQ and WURST-coslp-MQ approaches. NS = 96, RD = 0.5 s, Texp = 48 s each. On each spectrum the maximum intensities are given relative to that with ST-SPAM. MQ, MQ-SPAM, ST, ST-SPAM: 1 {p1, p2, p3, p4} = {121, 121, 10, 10} kHz. WURST-lp-MQ: 1 = 54, irr = 400 kHz. WURSTcoslp-MQ: 1 = 54, irr = 0, cos = 400 kHz. WURST80: WURST = 1 ms, R = sweep = 62.5, 1,WURST = 13, off = 250 kHz.

Fig. 7 .

 7 Fig.7. Cs4Nb11O30. Comparison of isotropic projections of 2D 93 Nb DQF-STMAS and hp-, WURST-lp-, WURST-coslp-MQMAS spectra with (a) z-filter and (b) full-echo (echo = 0.3 ms) acquisitions. RD = 0.5 s, Texp = 1.0 h each. On each spectrum the intensity is given relative to that observed with DQF-STMAS. hp-MQMAS/DQF-STMAS: NS = 96 (Z), 192 (FE), t1 = {16/32} s, Nt1 = 70 (Z), 35 (FE), 1 {p1, p2, p3, p4} = {100, 100, 20, 20} kHz. WURST-lp-MQMAS: 1 = 11, irr = 230 kHz. WURSTcoslp-MQMAS: 1 = 14, irr = 0, cos = 200 kHz. WURST80: WURST = 1 ms, sweep = R = 62.5, 1,WURST = 20, off = 400 kHz.

Fig. 8 .

 8 Fig.8. Simulated (a,b) 23 Na (I = 3/2) and (c,d) 27 Al (I = 5/2) signal intensity of hp-MQMAS and WURST-coslp-MQMAS, with respect to CQ and the RF-field (1) of the excitation/conversion pulses. hp-MQMAS: {p1, p2, p3, p4} s = {3.5, 3.5, 12.5, 25} for 23 Na, {4.0, 1.8, 8, 16} for 27 Al, with 1 (p3, p4) = 10 kHz. WURST-coslp-MQMAS: cos = 250 kHz. WURST80: WURST = 1 ms, sweep = R = 1,WURST = 20, off = 300 kHz.

  1 T wide-bore magnet at a Larmor frequency of 0 = 58 ( 85 Rb), 156.4 ( 27 Al) and 146.9 ( 93 Nb) MHz, equipped with HX MAS probes using either  = 3.2 mm rotors at R = 16 or 21 kHz or  = 1.3 mm at R = 62.5 kHz. The spinning stability was maintained by a MAS III unit within 10 Hz. Maximum RF fields of 1 = 43 ( 85 Rb) and 61 ( 27 Al) kHz were attained with  = 3.2 mm, and 1 = 130 ( 27 Al) and 100 ( 93 Nb) kHz with  = 1.3 mm. All powder samples were packed as purchased or as synthesized.

Acknowledgements

H.N. acknowledges JSPS Grant-in-Aid for Early-Career Scientists (JP20K15319) and JST, PRESTO Grant Number JPMJPR2276. Authors thank Z. Gan and I. Hung for their fruitful scientific discussions.

Conflicts of interest

There are no conflicts to declare.


domain: info math

text:
Carlotta Giannelli

Sofia Imperatore

Angelos Mantzaflaris

Dominik Mokriš

Leveraging moving parameterization and adaptive THB-splines for CAD surface reconstruction of aircraft engine components

Keywords:

Reconstruction of highly accurate CAD models from point clouds is both paramount and challenging in industries such as aviation. Due to the acquisition process, this kind of data can be scattered and affected by noise, yet the reconstructed geometric models are required to be compact and smooth, while simultaneously capturing key geometric features of the engine parts. In this paper, we present an iterative moving parameterization approach, which consists of alternating steps of surface fitting, parameter correction, and adaptive refinement using truncated hierarchical B-splines (THB-splines). We revisit two existing surface fitting methods, a global least squares approximation and a hierarchical quasi-interpolation scheme, both based on THB-splines. At each step of the adaptive loop, we update the parameter locations by solving a non-linear optimization problem to infer footpoints of the point cloud on the current fitted surface. We compare the behavior of different optimization settings for the critical task of distance minimization, by also relating the effectiveness of the correction step to the quality of the initial parameterization. In addition, we apply the proposed approach in the reconstruction of aircraft engine components from scanned point data. It turns out that the use of moving parameterization instead of fixed parameter values, when suitably combined with the adaptive spline loop, can significantly improve the resulting surfaces, thus outperforming state-of-the-art hierarchical spline model reconstruction schemes.

Introduction

Constructing continuous and flexible geometric models that are easy to shape and manipulate is a fundamental requirement in many applications. For example, they can serve for visualization and direct measurements within a design phase or a manufacturing process, as well as for simulations, i.e. performing numerical computations directly on the model or its components. In Computer Aided Design (CAD) and Computer Aided Engineering (CAE) software, spline representations are constructed in terms of B-splines or their non-uniform rational extension (NURBS). This choice relies on the properties of these functions, namely non-negativity, locality, and partition of unity, which results in free-form shaped and flexible geometric models that are easy to manipulate. However, real-world data sets do not usually consist of such models, since the data acquisition procedures, e.g. via laser scanners, diagnostic devices, or photogrammetry schemes, yield raw-data points. Furthermore, depending on the application and the collection method, these data can result in structured point grids, meshes, or scattered point clouds. Therefore, the design of continuous free-from geometric models relies on different data fitting procedures. In this paper, we propose an accurate adaptive approximation framework, applied to industrial scattered data acquired through optical scans of aircraft turbine blade components. We use adaptive spline models combined with a moving parameterization to ensure high flexibility, while simultaneously handling scattered parametric data configurations within the adaptive loop.

The first fundamental step of any parametric fitting method consists of assigning a parameter value to each data point. While this parameterization process plays a crucial role it is still an open research topic and several solutions have been proposed, that are not optimal in a universal way. Parameterization methods for triangulated surfaces were proposed in [START_REF] Floater | Parametrization and smooth approximation of surface triangulations[END_REF][START_REF] Floater | Mean value coordinates[END_REF] and later extended to point clouds [START_REF] Floater | Meshless parameterization and surface reconstruction[END_REF] and periodic surfaces [START_REF] Groiss | 27 variants of Tutte's theorem for plane near-triangulations and an application to periodic spline surface fitting[END_REF]. These methods rely on a barycentric mapping induced from the local neighborhoods of the points. Recently, new data-driven approaches, such as neural networks, have also emerged as viable options for tackling this problem. These include [START_REF] Laube | Deep Learning Parametrization for B-Spline Curve Approximation[END_REF][START_REF] Scholz | Parameterization for polynomial curve approximation via residual deep neural networks[END_REF] for univariate spline and polynomial approximation respectively, while a more general framework was proposed in [START_REF] De | Parameterization learning with convolutional neural networks for gridded data fitting[END_REF][START_REF] Giannelli | Learning meshless parameterization with graph convolutional neural networks[END_REF] for multivariate spline approximation. We refer the reader to the related survey articles [START_REF] Floater | Surface parameterization: a tutorial and survey[END_REF][START_REF] Zhu | A review of 3D point clouds parameterization methods[END_REF] for more details on the topic of parameterization.

Often the parameterization and fitting problems are treated separately, that is, first a fixed parameterization is constructed and is kept fixed throughout the fitting algorithm. Hoschek introduced in [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF] the idea that the parameterization must also be adjusted to the fitted surface: the so-called intrinsic parametrization (or parameter correction method) is (iteratively) computed starting from a certain initial approximation. The fitted points are projected onto the surface and the projected footpoints take the place of the previous parameter values. The re-fitted surface can then significantly improve the reconstructed geometric model. The core of the method relies on efficient footpoint projection, that is reduced to a non-linear minimization problem. The iterative approach of [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF] for finding the footpoint (i.e., the closest point on the surface) has been revisited in [START_REF] Saux | An improved hoschek intrinsic parametrization[END_REF] where a Newton-like method is proposed for the problem. As in all non-linear minimization schemes, the quality of the initial point, the nature of the objective function, and the different stepping strategies play a crucial role in the computation. Among the different contributions that incorporate parameter correction to improve spline approximation and fitting, we just mention the case of structured adaptive data fitting with T-splines [WZ13, SFF * 19] and of adaptive template mapping [START_REF] Sajavi | Template Mapping Using Adaptive Splines and Optimization of the Parameterization[END_REF].

As far as surface fitting is concerned, current CAD/CAE technologies rely mostly on tensor-product constructions, which pose several challenges for approximating scattered data point clouds, as well as dealing with noise or outliers or other local properties. To this direction, the authors in [START_REF] Raffo | Weighted quasi-interpolant spline approximations: Properties and applications[END_REF] propose a weighted quasi-interpolant, where the definition of suitable weights aims at recovering the local information of the raw data points. A penalized global least squares model is exploited in [MJM23, LYM * 23]. In particular, locality is tackled by automatically adapting the penalization power according to the data density. Regarding shape representations, the need to overcome the intrinsic limits of tensorproduct spline constructions resulted in the development of new adaptive spline spaces that allow local refinement. Among others, these include T-splines [START_REF] Sederberg | Tsplines and T-NURCCs[END_REF], locally refined (LR) Bsplines [START_REF] Dokken | Polynomial splines over locally refined box-partitions[END_REF], hierarchical B-splines (HB) [START_REF] Kraft R | Adaptive and Linearly Independent Multilevel B-Splines[END_REF] and their truncated formulation as truncated hierarchical B-splines [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF] (THB). The corresponding family of fitting methods typically accept as input a scattered point cloud, together with a suitable parameterization, and provide as output a locally refined spline geometric model, characterized by a certain prescribed accuracy. In particular, we mention the THB-spline methods developed in [KGZ * 14, BGS17, BGGS18] and [BGG * 22].

In the present work we follow up on the idea that, as the refinement proceeds in an adaptive setting, not only the geometric model but also the parameter values of each data point should be optimized. In fact, when adaptively designing a parametric geometric model, there is no evidence that the initial (possible) optimality of the parameterization is maintained when changing the approximation space. We revisit two existing surface fitting schemes, a global least squares approximation and a hierarchical quasi-interpolation scheme, both based on THB-splines. At each step of the adaptive loop, we update the parameter locations by solving a non-linear optimization problem to infer footpoints of the point cloud on the currently fitted surface. We compare the behaviour of different optimizers for the critical task of distance minimization, by also relat-ing the effectiveness of the correction step to the quality of the initial parameterization. In addition, we apply the proposed approach in the reconstruction of aircraft engine components from scanned point data. It turns out that the use of a moving parameterization instead of fixed parameter values, when suitably combined with the adaptive spline loop, can significantly improve the resulting surfaces, thus outperforming state-of-the-art hierarchical spline model reconstruction schemes.

The paper is organized as follows. Section 2 provides a brief overview of (TH)B-splines constructions. After presenting the scattered data fitting problem, Section 3 introduces adaptive THBspline approximations with moving parameters. In particular, two THB-spline adaptive fitting algorithms for scattered data [KGZ * 14, BGG * 22] are revisited with this new perspective. Section 4 shows the benefits of employing the proposed method on a selection of geometric models representing aircraft engine components.

B-spline constructions

In this section, we introduce the main concepts of polynomial Bsplines and their hierarchical extension.

B-splines

Let Ω = [a, b] ⊂ R be a real interval, d ∈ N a polynomial degree, k := d + 1 the corresponding order, and let t t t := [t 0 , . . . ,t n+k ] a vector of non-decreasing real values, with t k-1 ≡ a and t n+1 ≡ b, called knot-vector. Thereby, n + 1 univariate B-splines of degree d over [a, b] can be recursively generated [START_REF] De Boor | A Practical Guide to Splines[END_REF]. To this purpose, let β j,k : R → R indicate the j-th B-spline of order k ≥ 1, for each j = 0, . . . , n. For u ∈ R, if k = 1 ,

β j,1 (u) = 1 if u ∈ [t j ,t j+1 ) 0 otherwise and β j,1 (t n+1 ) = 1. If k ≥ 2, β j,k (u) = ω j,k (u)β j,k-1 (u) + 1 -ω j+1,k (u) β j+1,k-1 (u),
where ω j,k : R → R is a piecewise linear polynomial of the form

ω j,k (u) = u-t j t j+k-1 -t j , if u < t j+k-1 0 otherwise.
Note that each B-spline β j,k for j = 0, . . . , n is a piecewise polynomial function of degree d, which has maximum local smoothness on each subinterval of the partition t t t. On the other hand, the regularity at each unique knot t j ∈ t t t is dµ j , where 1 ≤ µ j ≤ d is the number of times the knot value t j appears in t t t. Univariate B-splines are characterized by three key properties: (i) nonnegativity,

β j,k (u) ≥ 0 for all u ∈ R, (ii) local support, β j,k (u) = 0 if u ̸ ∈ [t j ,t j+k ), and (iii) partition of unity, ∑ n j=0 β j,k (u) = 1 if u ∈ [t k-1 ,t n+1 ] ≡ [a, b].
The spline space of order k and knots t t t is defined as V := span{β 0,k , . . . , β n,k }.

Univariate B-splines can be easily extended to higher dimensions by considering a tensor-product construction. In particular, let Ω be a hypercube of R D , d d d = (d 1 , . . . , d D ) the polynomial multidegree, k k k = (k 1 , . . . , k D ) the corresponding polynomial multi-order and G the tensor-product mesh, defined by suitable knot vectors in each parametric direction t t t 1 , . . . ,t t t D . A tensor-product B-spline β j : R D → R is then defined as the product of D univariate Bsplines, for each j ∈ Γ k k k := { j = ( j 1 , . . . , j D ) | j h = 1, . . . , n h , h = 1, . . . , D}. Because of the tensor construction, many of the simple algebraic properties of univariate B-splines hold. In particular, nonnegativity, locality, and partition of unity. The tensor-product spline space of multi-order k k k with respect to the tensor-product grid G is the space

V := span{β j | j ∈ Γ k k k }.

Truncated hierarchical B-splines

When dealing with multivariate settings, the tensor-product constructions allow only a global distribution of the degrees of freedom for the problem at hand, whereas local spline refinement enables the possibility of achieving flexible and accurate models by strongly reducing the total number of degrees of freedom when compared with standard tensor-product B-spline representations. This type of mesh design is intrinsically supported by hierarchical splines [START_REF] Forsey | Hierarchical b-spline refinement[END_REF], where local refinement is achieved by introducing tensor-product B-splines on multiple hierarchical levels. A selection mechanism to properly identify a hierarchical B-spline basis was originally introduced in [START_REF] Kraft R | Adaptive and Linearly Independent Multilevel B-Splines[END_REF].

Let Ω be a hypercube of R D and let d d d be a polynomial multi-degree, while k k k is the corresponding multi-order. We consider a nested sequence of L tensor-product B-spline spaces

V 0 ⊂ . . . ⊂ V L-1 defined on Ω ⊂ R D , so that for each level ℓ = 0, . . . , L -1, V ℓ = span β ℓ j | j ∈ Γ ℓ k k k , where Γ ℓ k k k
is the set of indices for the tensor-product B-spline basis of level ℓ. In addition, each V ℓ is associated with a tensor-product mesh G ℓ and their (nonempty) quadrilateral elements q are commonly addressed as mesh cells of level ℓ. We also consider a nested sequence of closed domains Ω ≡ Ω 0 ⊃ . . . ⊃ Ω L = ∅, so that each Ω ℓ is the union of a subset of cells of the tensor product mesh G ℓ . The hierarchical mesh is defined as

G := q ∈ G ℓ | ℓ = 0, . . . , L -1 , where each G ℓ := q ∈ G ℓ | q ⊂ Ω ℓ \ Ω ℓ+1 is called the set of ac- tive cells of level ℓ.
The hierarchical spline construction consists in replacing any Bspline of level ℓ with support completely contained in Ω ℓ+1 by Bsplines at successively refined levels. More precisely, the hierarchical B-spline basis is defined as

H k k k := β ℓ j | j ∈ A ℓ k k k , ℓ = 0, . . . , L -1
where

A ℓ k k k := j ∈ Γ ℓ k k k | supp β ℓ j ⊆ Ω ℓ ∧ supp β ℓ j ̸ ⊆ Ω ℓ+1
is the set of indices of active functions and supp β ℓ j denotes the intersection of the support of β ℓ j with Ω 0 . The corresponding hierarchical space is defined as span {H k k k }. This simple definition, however, leads to different overlaps of coarse and fine B-spline supports, and partition of unity of the basis is lost. Truncated hierarchical B-splines (THB-splines) were introduced in [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF] to reduce the interaction between hierarchical functions at different levels and recover the partition of unity property. For any ℓ = 0, . . . , L -2, let s ∈ V ℓ ⊂ V ℓ+1 a tensor-product spline represented in terms of the basis of the refined space V ℓ+1 as

s(u u u) = ∑ j∈Γ ℓ+1 k k k c ℓ+1 j (s)β ℓ+1 j (u u u), for u u u ∈ Ω, for suitable coefficients c ℓ+1 j (s), for each j ∈ Γ ℓ+1 k k k . The truncation of s ∈ V ℓ at level ℓ + 1 is defined as trunc ℓ+1 (s) := ∑ j∈Γ ℓ+1 k k k supp(β ℓ j )̸ ⊆Ω ℓ+1 c ℓ+1 j (s) β ℓ+1 j
and the cumulative truncation with respect to all finer levels is

Trunc ℓ+1 (s) := trunc L-1 trunc L-2 . . . trunc ℓ+1 (s) . . . ,
with Trunc L (s) ≡ s, for s ∈ V L-1 . Finally, the THB-spline basis for the hierarchical space can be defined as

T k k k := τ ℓ j = Trunc ℓ+1 β ℓ j | j ∈ A ℓ k k k , ℓ = 0, . . . , L -1 ,
where the B-spline β ℓ j is the mother B-spline of the truncated τ ℓ j . THB-splines are non-negative, have local support, form a partition of unity, and span the same space of the hierarchical B-spline basis, span

{H k k k } ≡ span {T k k k } [GJS12]
. The effectiveness of employing THB-splines both for geometric design and isogeometric analysis has been shown in [GJK * 16], among others.

Because of their properties, THB-splines are a desirable tool for building flexible geometric models. They have been here presented in their more general form, where their construction can be developed in any parametric and physical dimension, namely for any D, N ∈ N ≥1 . In particular, a (TH)B-spline object is a linear combination of (TH)B-splines, defined as

s s s(u u u) = L-1 ∑ ℓ=0 ∑ j∈A ℓ k k k c c c j τ ℓ j (u u u), u u u ∈ Ω, (1)
with c c c j ∈ R N for j ∈ A ℓ k k k , ℓ = 0, . . . , L -1. For D = 1 a planar (N = 2) or spatial (N = 3) (TH)B-spline curve can be specified, whereas for D = 2 and N = 3 a (TH)B-spline surface can be assembled. Throughout the rest of the paper, by moving to industrial applications, we will always assume D = 2 and N = 3.

Adaptive THB-splines fitting with moving parameters

Let P = {p p p i ∈ R 3 | i = 1, . . . , m} be a (noisy) scattered point cloud of m data belonging to the physical space, the problem of surface fitting consists in finding a surface s s s : Ω ⊆ R 2 → R 3 , which approximates each point of P, hence s s s i ≈ p p p i for each i = 1, . . . , m, where s s s i is a point on the surface associated to the observation p p p i for each i = 1, . . . , m. This results in an exact match s s s i = p p p i , if we require the approximating model s s s to interpolate the data, for each i = 1, . . . , m.

Due to the THB-spline properties and their effectiveness as geometric design tools in CAD/CAE frameworks, we require s s s to be a THB-spline surface of a certain fixed bi-degree d d d and bi-order k k k as defined in (1). A crucial problem in any parametric surface reconstruction scheme consists in the data parameterization, namely assigning to each point p p p i ∈ R 3 a corresponding parameter value u u u i ∈ Ω ⊂ R 2 for i = 1, . . . , m in order to define the parametric set

U := {u u u i ∈ Ω ⊂ R 2 |u u u i = (u i , v i ) , i = 1, . . . , m}.
Since the parameters encode intrinsic characteristics of the surface representation, estimating a good point cloud parameterization is a fundamental and delicate issue. State-of-the-art THB-splines fitting methods address the parameterization problem, the construction of the hierarchical space, and the definition of the control net separately. In particular, given P, a suitable parameterization U is computed and, on such a fixed parameterization, both the hierarchical space and the approximating spline model are defined by addressing alternately the knot placement and the computation of the control net within an adaptive loop [KGZ * 14, BGS17, BGG * 22]. In this section instead, we propose a strategy to properly embed the parameterization within the adaptive fitting schemes.

Initial spline configurations

To approximate the industrial scattered point clouds examined in this paper, the heuristics developed in [START_REF] Floater | Meshless parameterization and surface reconstruction[END_REF] are employed as initial parametrization. These methods consist of two steps: given a point cloud P, firstly its boundary points P B ⊂ P are individuated and parameterized with univariate techniques, see e.g. [Lee89, PT95, Lim99, SA04, BÖKK20], so that their parametric values lay anticlockwise on ∂Ω. Subsequently, a graph-connectivity is determined, by defining for each interior point p p p i ∈ P \ P B a certain neighbourhood, e.g. by determining the k nearest points to p p p i for some appropriate k ∈ N ≥1 , among others. Thereby, the interior point parameters are computed as weighted convex combinations of the parameters of their neighbours. Note that, specific choices of neighbourhoods and weights considered in the convex combinations characterize the final parameterization and lead to different results, see [START_REF] Floater | Meshless parameterization and surface reconstruction[END_REF] for more details. In principle, however, any parameterization method can be employed as a starting parameterization to be further improved within the adaptive scheme.

In addition to the parameterization, the design of spline models is characterized by further unknowns, i.e., the knot lines placement, which defines the initial knot configuration, and, consequently, the control points of the geometric model s s s. In a more general context, the weights associated with the data observations, a smoothness functional and the real factor controlling its influence are all additional potential unknowns. However, their analysis goes beyond the scope of this paper. For interpolation problems, the amount of degrees of freedom is known a priori since it is constrained to the number of data observations. On the other hand, the input data P hereby considered are industrial real-world data, namely noisy and unevenly distributed. Hence, any interpolation scheme would potentially result in a costly model with as many degrees of freedom as the number of input points, being nevertheless numerically inaccurate. Therefore, we seek a geometric model s s s, which approximate the data P within a certain tolerance ε ∈ R >0 , in the sense that dist (s s s i , p p p i ) ≤ ε for each i = 1, . . . , m, where dist (•, •) is a suitable distance metric. In particular, given the data P, their parameterization U and a tensor-product mesh G, we exploit two methods for the computations of the control points: a penalized global least squares and a quasi-interpolation scheme. The final amount of degrees of freedom needed to obtain the desired accuracy is usually derived through an iterative process which performs adaptive spline refinement.

Hierarchical spline approximation

Once an initial parameter and mesh configuration is chosen, any adaptive approximation procedure is characterised by four main steps which are successively repeated: (1) computation of the approximation on the current mesh, (2) error estimation, (3) marking and (4) refinement strategies to suitably identify the adaptive mesh to be used in the next iteration of the adaptive loop. In particular, we revisit the adaptive global least squares method (LS) proposed in [KGZ * 14] and the adaptive hierarchical quasi-interpolation (QI) approach proposed in [BGG * 22] to define adaptive THB-spline surface fitting schemes with parameter correction.

Within the surface fitting framework, for a fixed THB-spline space, the first step of the adaptive loop consists in computing the control points c c c ℓ j for each j ∈ A ℓ k k k and ℓ = 0, . . . , L -1 to define the geometric model in (1). As concerns LS, this is achieved by solving the penalized least squares problem

min c c c ℓ j , j∈A ℓ k k k , ℓ=0,...,L-1 1 2 m ∑ i=1 ∥s s s (u u u i ) -p p p i ∥ 2 2 + λJ (s s s) ,
where the penalization term J is the thin-plate energy functional, whose influence is controlled by a weight λ ≥ 0., i.e. for u u u = (u, v) ∈ Ω,

J (s s s) = Ω ∂ 2 s s s ∂u∂u 2 2 + 2 ∂ 2 s s s ∂u∂v 2 2 + ∂ 2 s s s ∂v∂v 2 2 dudv. ( 2
)
As concerns the QI, each control point c c c ℓ j depends on a (local) subset of the input point cloud P j ⊂ P and its corresponding parameters U j ⊂ U, i.e. c c c ℓ j = c c c ℓ j P j , U j for each j ∈ A ℓ k k k and ℓ = 0, . . . , L -1 and therefore, the resolution of a global linear system is avoided. Note that by exploiting THB-spline constructions, it is possible to define hierarchical quasi-interpolation schemes without any additional efforts with respect to their tensorproduct formulations [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF]. More specifically, the computation of the QI control net consists of implementing a two-stage algorithm by combining local penalized least squares spline approximations (first stage) with the assembly of the hierarchical quasiinterpolant (second stage). In particular, for each j ∈ A ℓ k k k and ℓ = 0, . . . , L -1, let β ℓ j be the mother tensor-product B-spline of the truncated τ ℓ j , i.e., τ ℓ j = Trunc ℓ+1 (β ℓ j ), and Ω j ⊂ Ω a suitably chosen local subdomain which has a non-empty intersection with β ℓ j , namely Ω j ∩ supp(β ℓ j ) ̸ = ∅. Subsequently, individuate the indices I j = i | u u u i ∈ U ∩ Ω j ⊂ {1, . . . , m}, so that n min ≤ I j ≪m and the related data

P j = p p p i | i ∈ I j , U j = u u u i | i ∈ I j .
(3)

For more details on the choice of n min , I j , P j and U j see [BGG * 22].

Finally, denote with

β ℓ r | r ∈ Λ ℓ, j k k k ⊂ Γ ℓ k k k
the subset of tensor-product B-splines which do not vanish on Ω j (by definition β ℓ j belongs to it). We then approximate P j with

s s s j (u u u) = ∑ r∈Λ ℓ, j k k k α α α ℓ r β ℓ r (u u u) (4)
by computing the following penalized tensor-product B-spline local approximation, min

α α α ℓ r ,r∈Λ ℓ, j k k k ∑ i∈Ij s s s j (u u u i ) -p p p i 2 2 + λJ(s s s j ),
where the penalization term has been chosen again as the thin-plate energy in (2). The computation of all the local approximations (4) for each j ∈ A ℓ k k k and ℓ = 0, . . . , L -1 terminates the first stage. The second stage consists of assembling the global approximation defined in (1). In particular, due to the THB-spline properties, each coefficient c c c ℓ j = c c c ℓ j (P j , U j ) in (1) directly corresponds to the coefficient α α α ℓ j associated with β ℓ j in (4).

Subsequently, the second step of the adaptive fitting loop consists of evaluating the THB-spline approximant on the parameter sites u u u i ∈ Ω related to the data points p p p i to compute a suitable error indicator. In particular, we choose the point-wise error distance ∥s s s(u u u i )p p p i ∥ 2 for each i = 1, . . . , m, among others. The error indicator indicates the region of the domain Ω where additional degrees of freedom are needed to meet the prescribed surface accuracy, by individuating the parametric sites u u u i ∈ U where it exceeds a certain input threshold, i.e. ∥s s s(u u u i )p p p i ∥ 2 ≥ ε.

As concerns the LS scheme, for the making strategy we identify the cells of the current hierarchical level ℓ which contain the parameters u u u i identified by the error indicator and mark them for refinement, together with two surrounding rings of cells in the hierarchical mesh. For more details about the feasible configurations of hierarchical meshes with THB-splines for geometric design, see [GJK * 16]. In the QI scheme instead, the marking process is naturally implemented on the basis functions. In particular, the basis function of level ℓ which are active on the parameter sites characterized by ∥s s s(u u u i )p p p i ∥ 2 ≥ ε are marked for refinemnt and replaced by basis functions of the successive hierarchical level ℓ + 1.

Finally, the refinement strategy of the two THB-spline fitting schemes has to be considered. In LS, the marked cells are dyadically split, whereas in QI the refinement strategy is more sophisticated, to better handle scattered data configurations and exploit local adaptivity. In particular, in some situations, the parameter values corresponding to the local data set U j can be concentrated in a small part of the support of a marked function, thus splitting its support may affect the quality of the final approximation. The support of each marked basis function is then analyzed as follows: if it contains a minimum number of data points, its dyadic refinement is considered, otherwise the local refinement of the corresponding area is prevented, see [BGG * 22] for the details.

At this stage, the hierarchical space has been updated and both in [KGZ * 14] and [BGG * 22] a new iteration of the adaptive loops begins. We now enrich the hierarchical approximation scheme by updating the current parameterization, after the refinement strategy. This additional step in the adaptive loop enables the definition of adaptive THB-spline fitting with moving parameters.

Moving parameters

In the context of surface fitting, the parameters update within the described adaptive methodology is performed by adding a parameter correction (PC) routine [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF], at the end of the standard adaptive loop. This method consists of locating the points on the geometric model which are the closest to the data points, in terms of Euclidean distance. Given a point cloud P, its parameterization U and a surface s s s : Ω ⊂ R 2 → R 3 , the surface closest point problem consists in solving the following minimization problem, min

(ui,vi) 1 2 ∥s s s (u i , v i ) -p p p i ∥ 2 2 , for each i = 1, . . . , m.
This two-dimensional nonlinear problem can be explicitly formulated as

(s s s (u i , v i ) -p p p i ) T ∇s s s (u i , v i ) = 0, for each i = 1, . . . , m, (5)
where ∇s s s indicates the gradient of the surface s s s, and solved by employing a suitable optimizer. In view of (5), the vector connecting the data point p p p i to the surface point s s s(u i , v i ) has to be orthogonal to the tangent plane of the surface and s s s (u i , v i ) is then usually called the foot-point of p p p i over s s s for each i = 1, . . . , m. The updated parameterization ū u u i = ( ūi , vi ) is defined as the solution of (5), for i = 1, . . . , m. Note that after performing one step of PC, the geometric model can be updated by fitting again the surface s s s to the points P with the corrected parameters. Consequently, a new projection and the corresponding correction can take place. Usually, a few steps of PC are sufficient to improve the accuracy of the approximation.

The solution of the optimization problem leads to computing the optimal intrinsic parameterization of a given spline model. Within this framework, one step of PC consists of solving (5) for s s s belonging to the refined hierarchical space and subsequently projecting the points of P on the refined THB-spline surface s s s and finally considering their foot-point as the new parameters U. After a certain number of PC steps, the new parameters and the refined adaptive space are employed to implement another iteration of the adaptive loop. Note that the parameter correction is applied globally to the entire approximant. In particular, for the QI scheme, for each local problem the parameters U j individuated in (3) are kept fixed.

At the beginning of the next iteration, the corrected geometry is updated by solving the approximation scheme again and the adaptive procedure is iterated. As for standard adaptive fitting strategies, this loop is performed until the maximum point-wise error is within an input tolerance ε or a maximum number of hierarchical level L is reached. Note that PC steps should naturally be embedded in any adaptive scheme, i.e. not limited to fitting problems, since even if the starting parameterization benefits from optimality within the initial space, there is no evidence that optimality is maintained when moving to wider approximation spaces, related to their adaptive extensions.

Interaction of foot-point projection with adaptive spline configurations

Finding successfully the foot point projection of a point on a surface is a challenging open problem, see e.g. [START_REF] Ko | Orthogonal projection of points in CAD/CAM applications: an overview[END_REF] and references therein. As far as parameter correction is concerned, the foot-point projection in [START_REF] Hoschek | Intrinsic parametrization for approximation[END_REF] is performed iteratively, by linearizing the problem in (5). Subsequently, in [START_REF] Saux | An improved hoschek intrinsic parametrization[END_REF], the iterative procedure has been replaced by the application of a Newton-like approach. Due to the locality of Newton-like gradient methods, to properly correct the parameters, it is paramount to start from a good initial parameterization guess as well as with a reasonably good geometric model.

Figure 1 shows the consequences of incorporating the PC routine at a too-early stage of adaptive fitting schemes. More specifically, we fit a point cloud P of an industrial tensile part, about 2.5 • 10 -2 m long, and consisting of 9281 scattered data, with the LS scheme presented in [KGZ * 14]. An initial tensor-product Bspline approximation of bi-degree d d d = (2, 2) is built on a coarse 4 × 4 tensor-product mesh and illustrated in Figure 1 (a). This initial coarse mesh is then adaptively refined, using a criterion of error threshold with ε = 5e-5m, leading to a final THB-spline model with 1136 degrees of freedom (DOFs) that registers a maximum error (MAX) of 8.333e-5m. The THB-spline approximation obtained by introducing 1 PC step at each iteration of the adaptive loop is shown in Figure 1 (b) and is characterized by 969 DOFs, MAX error of 3.059e-4m. Consequently, we may note that using the inaccurate approximation shown in Figure 1 (a) to compute the foot-point projections and correct the parameters, leads to a final fitting result that is not suitable for further processing due to self intersections seen in Figure 1 (b). We then consider as initial spline configuration the tensor-product B-spline approximation built on a tensor-product mesh with two additional dyadic refinements with respect to Figure 1 (a), i.e., with a 16 × 16 mesh, as shown in Figure 1 (c). These settings lead to stable foot point projection, due to the quality of the initial approximation and the LS scheme with moving parameters results in a model with 607 DOFs and MAX error 8.467e-5m, shown in Figure 1 (d).

The choice of the optimizer and its settings are fundamental for the successful computation of the foot-point projections. Here, we consider the LS adaptive fitting method with moving parameters and compare the results with respect to the use of different optimizers and optimization settings. In particular, we start with the configuration of Figure 1 (c), i.e. bi-degree d d d = (2, 2) and a 16 × 16 tensor-product mesh, and apply 1 PC step at each adaptive iteration. For a better comparison, we also fix the total number of adaptive refinement steps. The results are shown in Tables 1. As a benchmark, we consider the commercial library Parasolid [Sie]. We can see that taking a too-big minimum step s in the gradient descent method (GD) leads to a final result that can be worse than not moving the parameters at all, with respect to to the number of DOFs of the final models and its accuracy in terms of MAX and MSE errors. However, if carefully fine-tuned, the results obtained by Parasolid in terms of DOFs, MAX and MSE, can be reached as in this case for s = 1e-12. In the following numerical examples, we employed the hybrid limited memory Broyden-Fletcher-Goldfarb-Shanno method (HLBFGS) [Liu] with minimum step length s = 1e-9, since it has always shown to provide extremely similar results to Parasolid in all the considered examples. The results obtained employing HLBFGS for the above-mentioned configuration are also shown in Table 1.

Finally, we remark that including the parameter correction rou- tine on a proper initial guess and choosing suitable settings for the optimizer leads to outperforming the state-of-the-art adaptive fitting schemes as we prove in the numerical examples of Section 4.

Numerical examples

Throughout this section, we use two fitting methods, namely the global adaptive least-squares approximation with THB-splines (LS) and the quasi-interpolation based scheme using local spline approximation (QI), both described in the previous section, to assess the benefits of the moving parameterization approach. We compare the results with and without moving parameters on scanned data provided by MTU Aero Engines. All experimentations were performed using the open-source library G+Smo [JLM * 14, Man20], available at https://github. com/gismo.

Example: Tensile

In this section, we revisit the second configuration of Example 1 in [BGG * 22], where the point cloud of the tensile part considered in the previous section, is initially approximated by a bi-quadratic spline on a tensor-product basis defined on a 4 × 16 mesh with tolerance ε = 5e -5m. We compare the performance of the QI scheme without parameter correction (0 PC) and with one or more rounds of PC after the adaptive refinement step. The HLBFGS optimizer with minimum step-size s = 1e-9 is used for the foot point projection. In particular, the resulting geometries when applying 0, 1 or 2 rounds of PC after each refinement step are shown in Figure 2. We can observe that all the fitting surfaces are of good quality. Nevertheless, in this case suitably embedding the parameter correction routine in the adaptive loop led to gaining the same precision in terms of MAX and MSE errors, while using fewer degrees of freedom to reconstruct the final geometric model, as summarized in the first three lines of Table 2. In particular, applying 3 steps of PC after each hierarchical mesh refinement produces an approximation with 15.20% fewer DOFs with respect to the one achieved with the standard adaptive QI scheme, while simultaneously reducing the MAX error by 9.89% and the MSE by 17.20%. Note that additional rounds (from 4 to 6 in Table 2) of parameter correction, yield diminishing improvements, since the DOFs are still reduced, but there is no gain in terms of MAX and MSE. This suggests that from 1 to 3 PC steps are in general a proper choice when including the parameter correction within adaptive spline fitting schemes.

Example: Blade

In this example, we revisit the second example presented in [BGG * 22]. However, in view of the considerations of Section 3.4, the initial tensor-product basis of bidegree d d d = (3, 3) is now more refined, i.e. we start with a 8 × 8 tensor product mesh. The results are still comparable, as the lowest level basis is entirely refined in the configuration considered in [BGG * 22]. We then fit the set of 27191 measured data points from a blade geometry that has a length of about 5 • 10 -2 m. We compare LS and QI, by always considering as smoothing coefficient λ = 1e -8. In both cases, we stop when at least 95% of data points are within the tolerance ε = 2e -5m.

Figure 3 shows the final THB-spline models obtained with LS (top) and QI (bottom) when 0, 1, or 2 PC steps are considered in the adaptive THB-spline fitting schemes. In particular, we can observe that the QI geometries are characterized by a smaller amount of oscillations along the sharp features. To have more insights into this phenomenon, we present the reflection lines on the different THB-spline models in Figure 4. By analyzing the top and bottom row of this figure from left to right, we see that the PC improves the surface quality by reducing the oscillations near the top left corner. Moreover, by comparing the LS (top) and QI (bottom) results in each column, we can see that the oscillations just under the feature in the top right area of the blade are reduced in the THB-spline models obtained with QI scheme.

The quantitative analysis for this example is conducted in Table 3. In particular, both for LS and QI schemes, we report the final percentage of points within the tolerance, the MAX and MSE errors, as well as the number of DOFs, when 0, 1, or 2 parameter corrections are considered. Note that the stopping criterion (95% pts < 2e -5m) is met in all cases, with a significant reduction (∼70%) of DOFs both for LS and QI schemes when 1 or 2 PC steps are applied. This is due to the fact that the prescribed precision is achieved two iterations earlier, preventing the introduction of two additional hierarchical levels. This can be explicitly seen in Figure 5 (LS) and 6 (QI) where the number of points below the prescribed tolerance versus the amount of DOFs at each adaptive iteration is reported. In this example, either one or two PC steps, after each refinement procedure, are a good choice, since the results are very similar, iteration by iteration, both for LS and QI.

Conclusions

In this work, we enhance the power of THB-spline approximations with an update of the point parameters in each adaptive fitting step, by suitably moving the parameterization within the iterative loop, both for least squares and quasi-interpolation schemes. Our study reveals that using the parameter correction step can improve the fitting results while also reducing the number of degrees of freedom required to achieve a certain accuracy. We observe that this correction can lead to earlier termination of the adaptive process, thus providing more compact models with less refinement depth. However, if used on a coarse fitting, the effect can be negative, since it can lead to an invalid parameterization and/or self-intersections in the model. It is then important to have an initial fit that is close enough to the point cloud. The robustness of the footpoint projection is also paramount for the overall process and needs to be accurate enough to gain quality and degrees of freedom.

Figure 1 :

 1 Figure 1: Adaptive THB-spline LS with moving parameters: effect of the initial tensor-product B-spline approximation. The THB-spline model (b) constructed from the initial configuration (a) defined on a 4 × 4 tensor-product mesh is affected by self-intersection due to the too coarse initial mesh. The THB-spline model (d) constructed from the initial configuration (c) defined on a 16 × 16 tensor-product mesh is successfully computed.

Figure 2 :

 2 Figure 2: Adaptive THB-spline QI with moving parameters in Example 4.1. From left to right: models computed with 0, 1, and 2 PC steps.

Figure 3 :

 3 Figure 3: THB-spline models with control nets in Example 4.2. Top row: LS, bottom row: QI. From left to right: 0, 1, and 2 PC steps.

Figure 4 :

 4 Figure 4: THB-spline models with reflection lines in Example 4.2. Top row: LS, bottom row: QI. From left to right: 0, 1, and 2 PC steps.

Figure 5 :Figure 6 :

 56 Figure 5: Percentage of points within the prescribed tolerance during the THB-spline LS scheme in Example 4.2.

Table 1 :

 1 Influence of optimizers used for the foot-point computation loop starting from configuration of Figure 1 (c).

	method	pts < ε	MAX (m)	MSE (m 2 )	DOFs
	no PC	92.770%	1.96820e-4	7.29588e-10	751
	GD, s = 1e -9	91.908%	2.92276e-4	1.06473e-9	748
	GD, s = 1e -10	96.875%	1.62567e-4	4.42720e-10	690
	GD, s = 1e -11	99.289%	8.47310e-5	2.40847e-10	607
	GD, s = 1e -12	99.310%	8.46760e-5	2.37395e-10	607
	HLBFGS	99.310%	8.46695e-5	2.37373e-10	607
	Parasolid	99.310%	8.46716e-5	2.37358e-10	607

Table 2 :

 2 Analysis in terms of points within the prescribed tolerance (%pts < ε), MAX or MSE, and final number of DOFs for the adaptive THB-spline QI scheme in Example 4.1. The percentages of the DOFs reduction with respect to QI without PC are also reported.

Table 3 :

 3 Analysis in terms of points within the tolerance (%pts < ε), MAX or MSE, and number of DOFs for the adaptive THB-spline LS and QI schemes in Example 4.2. The percentages of the DOFs reduction with respect to LS and QI without PC are also reported.

		100
		80
	pts. within tol.	40 60
		20				LS, 0 PC
						LS, 1 PC
		0				LS, 2 PC
		0	2,000	4,000	6,000	8,000
				#DOF

© 2023 The Author(s) Eurographics Proceedings © 2023 The Eurographics Association.

Acknowledgements

CG acknowledges the contribution of the National Recovery and Resilience Plan, Mission 4 Component 2 -Investment 1.4 -CN_00000013 "CENTRO NAZIONALE HPC, BIG DATA E QUANTUM COMPUTING", spoke 6. CG and SI are members of the INdAM group GNCS. The INdAM-GNCS support is gratefully acknowledged. CG, SI, AM also acknowledge the PHC GALILEE project 47786N and AM acknowledges H2020 Marie Skłodowska-Curie grant GRAPES No. 860843.


Gauthier Berthomieu
email: gauthier.berthomieu@univ-brest.fr.

François Salmon

Julian Palacino

Mathieu Paquier

Franc ¸ois Salmon

THE INFLUENCE OF THE NUMBER OF MICROPHONES CONSTITUTING A SPHERICAL MICROPHONE ARRAY ON THE PERCEPTION OF SPATIAL ALIASING

Keywords: Ambisonics, Spherical Microphone Array, Spatial Aliasing, diffuse-field equalization

    The influence of the number of microphones constituting a spherical microphone array on the perception of spatial aliasing

INTRODUCTION

Spherical Microphone Arrays (SMA) are commonly used for capturing sound fields in a physically precise way. Re-cent technology advances in this field allow to create spatial audio content that accurately copies the real-life audio environment. Sound fields recorded thanks to SMAs are commonly described in the Spherical Harmonics (SH) domain by the use of the ambisonics format. This formatting provides features such as rotating the sound field according to the listener's head orientation in real time before rendering it binaurally (dynamic binaural rendering).

Spatial resolution

When a SMA captures a sound field, the pressure that applies on the continuous surface elements (the surface of the sphere along which the microphones are distributed) is only captured at the discrete positions of the microphones. As a result, the SH representation of the sound field (which is theoretically composed of an infinite sum of discrete SH of different degrees and modes) is limited to a finite order. The physical accuracy of the represented sound field thus relies on the maximal SH order of the representation (referred to as the ambisonic order), which is limited by the number of microphones constituting the SMA.

Moreover, the restricted number of sensors leads to spatial ambiguities, where spatial information that belongs to higher-order harmonics are aliased into lowerorder harmonics [START_REF] Rafaely | Spatial Aliasing in Spherical Microphone Arrays[END_REF]. As higher orders harmonics contribute mostly to the representation of the higher frequencies (the contribution of harmonics of order n decay steeply for n > kr, where r is the SMA radius and k the wave number), spatial aliasing implies a corruption of the upper-frequency range of the recorded sound field. The aliasing frequency f A beyond which spatial aliasing is not negligible decreases as the distance between microphones increases [START_REF] Rafaely | Spatial Aliasing in Spherical Microphone Arrays[END_REF]. Hence, for a SMA, f A depends on the number of microphones, the radius and the sampling scheme according to which the microphones are arranged on the sphere.

Perceptual evaluation of the artifacts

The errors between the actual sound field recorded with a SMA and the reproduced sound field can lead to audible artifacts that degrade the quality of the listening experience. As an example, the truncation-induced error is known to reduce the localization accuracy as well as the externalization of binaural reproductions [START_REF] Avni | Spatial perception of sound fields recorded by spherical microphone arrays with varying spatial resolution[END_REF].

Some studies aimed at defining minimal SH orders beyond which the effect of encoding errors on the rendered audio content is negligible. Lübeck et al. [START_REF] Lübeck | Binaural reproduction of dummy head and spherical microphone array data-A perceptual study on the minimum required spatial resolution[END_REF] evaluated Points of Subjective Equality (PSE) between binaural auralizations synthetized from impulse responses captured through SMAs of varying ambisonic orders (inducing varying truncation and aliasing error) and a binaural reference synthetized from a SMA of order N = 29 (for which the introduced artifacts can be neglected). Moreover, they evaluated PSEs between the reference and binaurals signals synthetized from Binaural Room Impulse Response (BRIR) SH interpolations, which suffer from truncation error but for which each individual measurement has the maximum spatial resolution and thus do not exhibit artefacts due to spatial aliasing. Following an ABX three-interval/two-alternative forced choice (3I/2AFC) test design, the authors reported that orders N = 17 -20 were required for the auralizations of SMA data to be indistinguishable from the reference, whereas orders N = 9 -13 were required for the auralization of BRIR data. However, these PSEs were obtained using a baseline approach, and several methods have been developped in the last years to improve the binaural renderings of SMA captures, such as the diffuse-field equalization.

Diffuse-field equalization

As the aliasing error increases, the sound field is not correctly represented in the SH domain. The aliased components induce an increase of the high frequency energy, leading to an audible spectral coloration. The diffusefield equalization step aims at correcting the aliased signal spectrum to that of a diffuse sound field beyond a given frequency.

This equalization was performed on the encoding matrix (which allows to convert the microphonic signals into SH signals) in the samme manner as in McCormack et al. [START_REF] Mccormack | Parametric ambisonic encoding of arbitrary microphone arrays[END_REF] for f > f A .

EXPERIMENT

The present experiment aims at evaluating the minimal number of microphones constituting a SMA beyond which the spatial aliasing error is no longer perceptible. An adaptive procedure was carried out, consisting in pairwise comparisons between spatially aliased stimuli and reference signals that do not exhibit spatial aliasing. SMAs of different radii and sampling schemes were investigated.

The aliased stimuli were computed using different virtual SMAs and reproduced through headphones following a dynamic binaural decoding. In order to evaluate the effect of the diffuse-field equalization on the reproduced signals, the experiment was carried out with both equalized and non-equalized encoding matrices.

Sound sources

Spatial Room Impulse Responses (SRIR) were simulated to model the incident sound field reaching the SMA. A small (5.3 m × 4.1 m × 2.3 m -width × depth × height) room was modeled as a shoe-box, in which five sound sources were placed with an azimuth angle of ±40 • , ±15 • and 0 • from the recording position. The reverberation time of the room was equal to 0.24 s. These SRIR were simulated at the order N = 33 so that the the incident sound field was physically accurate.

To create virtual SMA signals, the microphone signals were first computed using different sampling grid orders N g (which defines the number of microphones), sphere radii and sampling schemes depending on the SMA configuration from the SRIRs of ambisonic order N = 33, and encoded in ambisonics at the order Ñ = 7.

The unaliased reference stimuli were created by directly truncating the simulated SRIRs of order N = 33 at the order Ñ = 7. The resulting SRIRs were then convolved with anechoic sound sources to create the stimuli. In summary, aliased ambisonics signals of order Ñ = 7 exhibiting different amount of spatial aliasing induced error (depending on the characteristics of the SMA and the diffuse-field equalization step) were compared to unaliased ambisonic signals of order Ñ = 7. The order Ñ = 7 was chosen since it has been reported to provide a reasonable trade-off between technical effort and perceptual quality [START_REF] Bernschütz | Microphone arrays and sound field decomposition for dynamic binaural recording[END_REF].

In total, 60 virtual SMAs were created by combining 2 radii (6 cm and 8.75 cm), 15 sampling grid orders (7 to 21) and 2 sampling schemes (Gauss-Legendre [START_REF] Bogaert | Iteration-free computation of gausslegendre quadrature nodes and weights[END_REF] and Extremal [START_REF] Sloan | Extremal systems of points and numerical integration on the sphere[END_REF]). For these sampling schemes, a minimum number of sampling points Q = (N g + 1) 2 and Q = 2 * (N g + 1) 2 are required for extracting SH coefficients at a given sampling grid order N g , respectively. The virtual SMA properties and the corresponding aliasing frequency are shown in Table 1. For these 60 different SMAs, two encoding matrices were created: one being diffuse-field equalized and one not being diffused-field equalized.

Two audio contents were considered in the experiment: a 2 s-long pink noise and a 4.5 s-long musical excerpt of Don't mean a thing from Duke Ellington, played by a quartet comprising a violin, a tenor saxophone, a guitar and a double-bass.

Sound sources from the quartet were located at -40 • , -15 • , 15 • and 40 • in the horizontal plane. The audio signals were convolved with the corresponding SRIRs and the resulting signals were summed to synthesize the 7thorder ambisonic scene. The pink noise signal was convolved with SRIRs computed for a sound source located at 0°in the horizontal plane.

Stimuli were presented over open, circumaural headphones (Sennheiser HD650). The ambisonic signals were converted into binaural headphone signals using the Ambi Head HD VST plug-in from Noise Makers [START_REF] Makers | AmbiHeadHD[END_REF] with HRTFs derived from measurements of a Neumann KU 100 dummy head [START_REF] Bernschütz | A spherical far field HRIR/HRTF compilation of the neumann ku 100[END_REF]. This decoder implements the binaural rendering approach of Schörkhuber et al. [START_REF] Schörkhuber | Binaural rendering of ambisonic signals via magnitude least squares[END_REF] and Zaunschirm et al. [START_REF] Zaunschirm | Binaural rendering of ambisonic signals by headrelated impulse response time alignment and a diffuseness constraint[END_REF]. Head-tracking was performed using the T3 head-tracker from Feichter Electronic [START_REF] Electronics | T3 Head-tracker[END_REF] along with the Ambi Head HD VST plug-in.

Procedure

This experiment intended to highlight to what extent the aliasing artifacts were perceptible in different SMA configurations. Here, the aliasing frequency f A was shifted up or down by increasing or reducing the grid order N g of the SMA, respectively. This was done in an adaptive procedure which aimed at revealing a PSE, in terms of grid order, beyond which the induced aliasing artifacts were not detected by the listeners.

PSEs were evaluated independently for each combination of the SMA characteristics and for equalized and non-equalized encoding matrices. In total, the procedure consisted in a repeated-measures design gathering 16 runs (2 radii × 2 grids × 2 equalization conditions × 2 audio contents).

Each run started at a grid order N g = 7 (i.e. with the lowest number of microphones for the corresponding sampling scheme given Table 1), and the grid order was either increased or decreased after each trial following a one-up-one-down staircase method.

Each trial of the run consisted in a threeinterval/three-alternative forced choice (3I/3AFC) test. During each trial, three intervals (A, B and C) were presented to the listener. Each interval could be assigned either the stimulus which grid order varied accordingly to the procedure (the test stimulus) or the unaliased reference (the reference stimulus). During each trial, two intervals were assigned the same stimulus (either the test or reference stimulus) and the third interval (named the oddball) was assigned the other stimulus. During each trial, each of A, B and C was played once in this order. Participants were instructed to designate which of A, B or C was the oddball.

The grid order was increased by one for the next trial of this run if the listener correctly identified the oddball whereas it was decreased by one if the listener gave the wrong answer. Each run ended after 9 reversals (a reversal is defined as a correct answer followed by a wrong answer, or vice versa).

Playback, commands, and data capture were controlled by a software implemented in Max 8.

The experiment consisted of two 30 min-long sessions (one per sampling scheme condition). A training phase allowed subjects to become thoroughly familiar with the test environment and the grading process prior to formal grading.

Pre-results

The pre-results presented in this paper were gathered on 10 expert participants who were related to the Sound Master's degree of the University of Brest.

Fig. 1 shows that grid order PSEs were lowered by the diffuse-field equalization applied to the encoding matrices for each SMA configuration and each recorded audio content. Grids of order as low as N g = 7 -10 could lead to auralizations of the music excerpt that were indistinguishable from unaliased references providing the encoding matrices were equalized. When the encoding matrices were not equalized, SMAs required N g > 10 in each configuration for the aliased stimuli to be indistinguishable from the unaliased references. As Lübeck et al. reported that truncation error might be audible up to the orders N = 9 -13 [START_REF] Lübeck | Binaural reproduction of dummy head and spherical microphone array data-A perceptual study on the minimum required spatial resolution[END_REF], correcting the aliasing-induced error thanks to a diffuse-field equalization might allow both truncation and aliasing induced artifacts to no longer be audible beyond these orders.

Figure 1 .

 1 Figure 1. Grid order PSE of SMAs using an Extremal (Ext) or Gauss-Lengendre (G-L) sampling scheme, recording noise (left) and musical excerpt (right) with (blue) or without (brown) diffuse-field equalization.

Table 1 .

 1 Number of microphones involved in the listening test for each sampling scheme (left) along with the corresponding aliasing frequency for each radius (right).

		Sampling grid	Radius
		Extremal Gauss-L.	6 cm	8.75 cm
	N g	Microphones		Aliasing freq. (Hz)
	7	64	128	6368	4367
	8	81	162	7278	4991
	. . .	. . .	. . .	. . .	. . .
	21	484	968 19106	13101

Turin, Italy • 11 th -15 th September 2023 • Politecnico di Torino